INTRODUCTION
Previous studies in sports science suggest that regular exercise is beneficial for our health, for example, it reduces the chances of developing heart disease [1–4].
Nevertheless, the consequences of performing exhausting endurance activities remain unclear. For example, a study looking at Tour de France participants found a significantly lower rate of cardiovascular events and longer life expectancy compared to agematched, general French men population [5]. Conversely, subjecting the body to the extremes of rigorous training, such as marathon running, might lead to adverse effects [6, 7]. The reasons behind these divergent outcomes may lie in various factors that influence how our bodies respond to intense physical exertion. It has been posited that age and the conditions under which training occurs could account for these differences [8].
To better understand what happens in the human body, the present study analyzed changes in the gene expression of subjects after they had run a marathon. Previous work focusing on gene expression during an ultra-marathon reported a significant impact on the immune and inflammatory systems [9]. Furthermore, enriched pathways were linked to protein synthesis repression, an altered immune system, and infectious disease-related mechanisms. Some studies suggest that strenuous exercise can induce immune system dysfunction and increase the risk of infection while endurance activities could also have a direct impact on fatty acid metabolism [10, 11]. High-intensity training may increase the prevalence and severity of coronary atherosclerosis and modify insulin resistance in humans [12, 13].
Consequently, research into the impact of endurance sports on health, both at a professional or amateur level, has grown recently. Therefore, the aim of this study was to assess the effect of a demanding intense physical effort on gene expression using RNA extracted from whole blood samples, with a particular focus on understanding the mechanisms involved in recovery to mitigate any negative effects. This was done by identifying genes with significantly different expression levels after the subjects had run a marathon; then studying the genes’ biological implications through pathway and gene ontology term enrichment analyses and assessing their significance during the race; and finally defining their expression recovery levels 24 hours after the strenuous exercise.
MATERIALS AND METHODS
Study Design
A total of 10 mL of whole blood was collected from 78 non-elite athletes at three distinct time points surrounding their participation in the 2016 Barcelona Marathon, a 42-kilometer endurance event. After comprehensive quality control (QC) assessment, which considered RNA concentration, final library concentration, and sequencing performance metrics (see Supplementary Table 1), 18 participants were excluded due to poor RNA quality or incomplete extractions. The final dataset comprised 60 individuals (42 men and 18 women), with participant characteristics summarized in Table 1. The male group showed a mean height of 178.9 ± 5.5 cm, weight of 75.6 ± 7.3 kg, and finishing time of 205.8 ± 25.3 minutes, while the female group presented 165.5 ± 6.2 cm, 57.1 ± 7.0 kg, and 226.5 ± 19.0 minutes, respectively. The final dataset included 180 samples (three per participant) encompassing 14,554 genes annotated with Ensembl IDs.
TABLE 1
Participant characteristics by sex. This table summarizes the demographic and performance characteristics of the study participants stratified by sex. Variables include Sex, Age group (years), Height (cm, mean ± SD), Weight (kg, mean ± SD), and Finishing Time (minutes, mean ± SD). Age groups are presented with the number of participants and the corresponding percentage of the total within each sex.
Blood samples were collected at the following time points: immediately prior to the race (START), immediately after the race (FINISH), and 24 hours post-race (24REST). These time points enabled three comparative analyses: C1, which compared FINISH to START levels; C2, which examined FINISH against 24REST levels; and C3, which assessed 24REST relative to START levels (Figure 1).
FIG. 1
Experimental design and analysis workflow. (1) Blood samples were collected from participants at three time points: START, FINISH, and 24 hours after the marathon. (2) RNA was extracted from the samples and subjected to high-throughput sequencing. (3) Differential gene expression analyses were performed in parallel using edgeR and DESeq2. (4) Significantly regulated genes were subsequently analyzed for functional enrichment using Gene Ontology (GO) terms and KEGG pathways. Finally, genes exhibiting similar expression patterns across the time points were clustered to identify groups of co-regulated genes.
The study was conducted in accordance with the Declaration of Helsinki and was approved by the Ethics Committee at Hospital de la Santa Creu i Sant Pau (Barcelona). All participants provided written informed consent. To ensure compliance with data protection regulations, all data was managed under the framework of the General Data Protection Regulation (GDPR), safeguarding participants’ rights regarding the processing and free movement of personal data.
RNA Extraction and Sample Preparation
Whole blood samples were collected in PAXgene Blood RNA Tubes (QIAGEN GmbH, Germany) following the manufacturer’s recommendations. Intracellular RNA was isolated using the QIAsymphony PAXgene Blood RNA Kit (QIAGEN GmbH), and total RNA was stored at -80°C until further use.
Blood samples from the FINISH time point were collected in a pavilion near the marathon finish line, equipped for blood extractions, while START and 24REST samples were obtained at Hospital de la Santa Creu i Sant Pau. All samples were sent to Catalonia’s Blood and Tissue Bank (Banc de Sang i Teixits, BST) for RNA extraction and sequencing. Globin mRNA was depleted using the GLOBINclear Kit (Invitrogen, Thermo Fisher Scientific Inc., Massachusetts, USA). RNA libraries were prepared using the Illumina Stranded mRNA Prep Kit (Illumina, San Diego, CA) according to the manufacturer’s protocol and sequenced on a NextSeq500 system (Illumina).
Gene expression quantification and linear model fitting approaches
Sequencing was performed using a 150-cycle (2 × 74 bp) paired-end read kit. Reads were aligned to the GRCh38 reference genome with the STAR (v.2.7.10a) algorithm using default parameters, and gene quantification was conducted with Salmon (v.1.10). Gene annotation followed Ensembl guidelines.
Gene expression data were processed using two complementary analytical pipelines: the first combining the EdgeR (v4.0.16) and limma (v3.58.1) packages, and the second employing DESeq2 (1.44.0) [14, 15, 16]. For the EdgeR/limma approach, lowly expressed genes with fewer than five counts were excluded to enhance reliability and avoid lowly expressed bias. The countsper-million matrix was normalized using the weighted trimmed mean of M-values (TMM) method implemented in EdgeR, which adjusts for library size differences. Batch effects were subsequently removed using the removeBatchEffect function from limma.
For the DESeq2 pipeline, the default normalization procedure was applied, in which counts are divided by sample-specific size factors determined by the median ratio of gene counts relative to the geometric mean per gene. In this approach, batch effect correction was not performed as a preprocessing step but was instead included as a covariate in the model. Following normalization, a multidimensional scaling (MDS) analysis was performed to examine sample distribution and identify potential covariates or clustering patterns for consideration in subsequent analyses.
The primary objective of the analysis was to compare gene expression levels across three time points using the pairwise contrasts: C1 (FINISH − START), C2 (24REST − FINISH), and C3(START − 24REST). A design matrix was constructed for each analytical pipeline to incorporate covariates that could potentially bias the model. For the limma/edgeR pipeline, the design matrix included group (Vg), sex (Vs), finishing time (Vf), and age (Va), with the intercept omitted to enable direct estimation of relative changes across all time points without dependence on a reference level:
Individual variability was not included as a fixed effect. Instead, we applied the duplicateCorrelation function prior to model fitting to estimate and account for the correlation among repeated measures from the same subject.
For the DESeq2 pipeline, the same covariates were included, along with the batch effect variable (Vbe). In this case, the intercept was retained, as DESeq2 does not require its removal:
Results from both analyses were compared to validate the robustness of DEG detection and confirm reproducibility of the transcriptomic patterns across time points.
Differential gene expression analysis (DEG)
Differential expression analysis was performed for all three pairwise comparisons. In the EdgeR/limma pipeline, the lmFit function from the sva (v.3.52.0) package was used to fit a linear model to each gene. Based on these linear models, the eBayes function was applied to compute moderated t-statistics, moderated F-statistics, and logodds of differential expression through empirical Bayes moderation of the standard errors toward a common value. This approach enabled the assessment of quantitative differences in gene expression levels between experimental groups.
For the DESeq2 pipeline, the DESeq function was used to identify differentially expressed genes. All p-values were adjusted using the Bonferroni multiple comparison correction method, a conservative approach that minimizes false-positive results [17]. This correction was applied consistently across all analyses conducted in this study.
Likelihood Ratio Test (LRT) analysis in DESeq2
In addition to pairwise comparisons, a Likelihood Ratio Test (LRT) was performed to identify genes exhibiting expression changes across the different conditions. The interaction between gene expression variation and the variables sex and marathon finishing time were also evaluated to determine whether specific genes exhibited differential expression patterns across groups as a function of these covariates. For this purpose, two additional LRT models were fitted in DESeq2: one testing the interaction between the different groups and marathon finishing time, and another testing the interaction between groups and gender, while adjusting for the other covariables in both cases.
Furthermore, significant DEGs identified across the different groups were clustered using a k-means algorithm to detect groups of genes with similar expression patterns [18]. The optimal number of clusters was determined using the elbow method [19]. Beyond general clustering, a subset analysis was conducted by intersecting the LRT results with cluster C3 to assess whether genes not recovered 24 hours after the marathon displayed distinct expression dynamics.
Finally, several visualizations, including hierarchical clustering, were generated using the gplots (v3.5.1) package [20].
Functional Enrichment Analysis: Gene Ontology and KEGG Pathways
Gene Ontology (GO) enrichment analysis was performed using the GOstats (v2.68.0) package to identify both overrepresented and underrepresented GO terms [21]. Significant genes obtained from the linear models described in the previous section (i.e., genes with adjusted p-values < 0.05) were used as input for the GO analysis of each comparison group (C1, C2, and C3). A conditional hypergeometric test was applied to evaluate relationships among GO terms. Using a subclass of HyperGParams, the package computed hypergeometric p-values to assess the over- or underrepresentation of each GO term within the specified gene set.
Results from GOstats were filtered based on three main parameters: (i) the total number of genes included in each GO term (size), (ii) the number of observed genes from the input set (count), and (iii) the probability of observing these genes by chance (odds ratio). The filtered results were visualized using the GOplot (v1.0.2) package [22].
For the KEGG pathway analysis, the same methodology described above was applied. The hypergeometric test was performed using the signatureSearch (v1.16.0) package [23], specifically employing the enrichKEGG2 function, which generates KEGG pathway enrichment results from a provided vector of gene identifiers.
RESULTS
Exploring covariates using multidimensional scaling analysis
To identify potential covariates that might influence the analysis, a multidimensional scaling analysis was conducted. This analysis revealed variations among samples based on normalized RNA counts, as illustrated in Figure 2.
FIG. 2
Principal component analysis (PCA) plot showing the clustering of samples based on transcriptomic profiles. Samples are colored by group: red for START, green for 24REST, and black for FINISH. Each point represents an individual sample, with “Δ” denoting male and “○” denoting female. The X-axis (leading logFC dim 1) explains 20% of the variation, while the Y-axis (leading logFC dim 2) explains 6%. Distinct clustering is observed for the FINISH group, while START and 24REST show more overlap.
In dim1 (X-axis), a distinct separation between the START and 24REST samples was evident in comparison to the FINISH samples. Regarding the covariates, sex was identified as a variable influencing inter-variability, given the clear difference in dim2 between females and males. Additionally, intra-variability among samples also requires consideration.
Differential gene expression analysis (DEGs)
Differential gene expression analysis revealed that 9,861 of the initial 14,554 genes were differentially expressed (DE) in C1, 9,730 in C2, and only 279 in C3. Given how the groups were compared, the overexpression in groups C1 and C2 determines that expression was significantly higher after the race (i.e., after strenuous exercise) than when at rest. While the overexpression of gene levels in comparison C3 indicates there was significantly higher level of gene expression 24 hours after the race compared to baseline levels. Finally, a deeper analysis of the DEGs was also performed for each group comparison, with a particular emphasis on the biological function of the most significant genes.
The genes in this comparison with significantly different expression levels were mainly associated with immune cell markers, chemokines, and interleukins (Table 2). Regarding immune cell markers, CD48, CD19, LRRC8C, LRRC8D and CLEC10A genes were found to be downregulated.
TABLE 2
Differentially expressed genes related to the immune system, ROS environment, and mitochondria in the START vs FINISH comparison C1. Two methods were used, DESeq2 and Limma, and p-values were adjusted by Bonferroni in both cases.
As for the chemokines, most of the downregulated genes were expressed by T cell lymphocytes (CCR4, CCR5, CCR6, CCR9). These cells play an important role in immunity and have also been found to be altered in similar studies [24]. Conversely, not all the genes related to cytokines were downregulated. Important specific cytokine receptors for monocytes, such as CCR2 and CCR1, were found to be overexpressed besides CCL5 being downregulated. Moreover, some genes related to the family of suppressors of cytokine signaling proteins were also downregulated (SOCS2, and SOCS7), while others were overexpressed (SOCS3 and SOCS6).
The last family of genes that stood out as DEGs was interleukins (IL). The most downregulated interleukins included IL-32, IL-10RA, ICOS, and IL5RA. Conversely, IL1R1, IL-18, IL-10RB, IL6R, IL-4R, and IL27 were overexpressed in the FINISH group compared with the START group.
Table 2 also includes proteins that regulate oxidative stress, including the glutathione peroxidase family (GPX7) and iron/manganese superoxide dismutase (SOD1), which were also downregulated, while GPX3, another gene from the glutathione peroxidase family and DNMT3A, which is an important gene related to reducing the oxidative environment was overexpressed.
Subsequently, several differentially expressed genes (DEGs) associated with inflammatory markers, as reported in previous studies were identified [25]. Notably, 21 out of 23 previously documented inflammatory markers were found to be significantly differentially expressed (Table 3).
TABLE 3
Differentially expressed genes related to inflammatory markers in the START vs FINISH comparison C1. Two methods were used, DESeq2 and Limma, and p-values were adjusted by Bonferroni in both cases.
The results of the C2 comparison were very similar to those of C1. All the genes reported in C1 were also differentially expressed in C2.
As previously stated, only 279 genes were differentially expressed in C3, representing less than 5% of the differentially expressed (DE) genes identified in the other two comparisons. In C3 (Table 4), genes related to endocytosis, and vesicle biology such as APP, VPS13D, and LRP1 were overexpressed, whereas TMEM179B, EMP3, and COPE were downregulated. Meanwhile, mitochondrial function and metabolism genes, including NDUFS7, PHB2, and ALKBH7, were downregulated. On the other hand, RNA processing, and ribosomal components, such as EIF3A, PRRC2C, and YLPM1 were overexpressed, while RPL28, GET3, and CARHSP1 were downregulated. Finally, transcriptional regulation and chromatin remodeling genes ARID1A and AHNAK were overexpressed, whereas ZNF581 was downregulated.
TABLE 4
The top 10 differentially expressed genes in START vs 24REST comparison (C3). Two methods were used, DESeq2 and Limma, and p-values were adjusted by Bonferroni in both cases.
Likelihood Ratio Test (LRT) and clusterization between groups
Following the likelihood ratio test, 16,401 significant genes were identified and used to determine clusters of similar expression patterns. The elbow method indicated the presence of two major clusters, and k-means clustering was subsequently applied with this specification. A heatmap incorporating hierarchical clustering via the complete linkage method was generated to visualize these clusters. The results revealed two distinct clustering of FINISH samples, whereas START and 24REST samples exhibited more similar expression profiles, indicating substantial transcriptomic alterations immediately after the marathon. Two primary patterns were observed: one group of genes upregulated at FINISH and another downregulated at FINISH, both returning toward baseline levels at 24REST (Figure 3a).
FIG. 3
Hierarchical clustering and k-means analysis of gene expression trajectories.
Gene expression trajectories were explored using k-means clustering, with the optimal number of clusters determined via the elbow method. Figure 3A shows two predominant clusters: a first group with increased expression immediately post-marathon that returned to baseline after 24 hours, and a second group exhibiting decreased expression post-marathon followed by recovery at 24 hours. Figure 3B highlights a smaller subset of genes (279 genes) displaying three distinct temporal patterns, from differential expression genes 24 hours after the race compared to baseline.
To examine finer subclusters, the analysis was repeated for the 279 genes identified in pairwise comparison 3. In this case, the elbow method suggested three clusters, which were visualized (y-axis) in a heatmap ordered by time group along the x-axis (Figure 3b). The first cluster (cluster 3) showed increased expression at FINISH, further rising at 24REST. The second cluster (cluster 4) exhibited elevated expression at FINISH compared with START, followed by a pronounced decrease at 24REST, falling below baseline levels. The third cluster (cluster 5) displayed decreased expression at FINISH relative to START, with partial recovery at 24REST that did not reach baseline levels.
Consistently, visualization of C1 gene expression stratified by sex (Figure 4) reproduced the five cluster patterns observed in the overall dataset.
FIG. 4
Temporal expression patterns of representative gene subsets across immune pathways.
This figure comprises five subpanels (a–e), illustrating the temporal variation in representative genes across key functional categories: immune markers, chemokines/SOCS, interleukins, Th1/Th2 signaling, and genes with no baseline-level expression at 24 hours.
Finally, to evaluate potential sex-specific transcriptional responses, an interaction model including sex and group was tested using DESeq2 while adjusting for age and batch effects. After bonferroni correction, no genes showed statistically significant interaction effects, indicating that the overall temporal expression changes following the marathon were similar in both men and women
Meanwhile, when comparing potential performance-specific transcriptional responses, the interaction between the model including finishing time and group display IRAK4 (ENSG00000198001) as the only single differentially expressed gene, with a bonferroni adjusted p-value of 0.04 but with a very low logFC -0.0013.
Gene ontology term and KEGG pathway enrichment analyses
After performing conditional GO analysis and filtering the results by size and counts, 310 enriched GO terms were found in C1, 318 in C2, and 43 in C3. A complete list of all GO terms is given in Tables S6 for C1 and S7 for C2. C1 and C2 returned similar results, as previously reported in the DEGs analysis. The top enriched GO terms for C2, related to mitochondrial functions (Figure 5a) and T and B cell activation (Figure 5b) were mainly downregulated while virusrelated (Figure 5c) GO terms were upregulated.
FIG. 5
GO enrichment analysis results in C2. A) Enriched GO terms related to mitochondrial functions in FINISH vs 24REST conditions with an adjusted p-value Bonferroni of < 0.05. B) Enriched GO terms related to immune system in FINISH vs 24REST conditions with an adjusted p-value Bonferroni of < 0.05. C) Enriched GO terms related to viruses (c) in FINISH vs 24REST conditions with an adjusted p-value Bonferroni of < 0.05.
Finally, results obtained in C3 (Table S8) show that most of the enriched GO terms were related to energy generation processes, such as the electron transport chain, several GO terms related with the ATP synthesis and the mitochondrial respiratory chain complex I assembly.
On the other hand, the total number of KEGG pathways enriched were 58, 57, and 12 for comparisons C1, C2, and C3, respectively. Tables S9 and S10 show the total list of pathways for comparisons C1 and C2, while the list of enriched pathways for C3 is shown in Table 5.
TABLE 5
Enriched KEGG pathways in the START vs 24REST comparison (C3). All p-values were adjusted by Bonferroni method.
The pathway analysis results closely resembled those of the GO term analysis, strengthening the significance of the findings. Moreover, the results of the C1 and C2 comparisons were also similar to each other but differed from C3. For comparisons C1 and C2, the pathways related to apoptosis, cellular senescence, mitophagy, and necroptosis were enriched. Further top enriched terms were a group of pathways related to lipid metabolism, including fatty acid metabolism, lipid and atherosclerosis, and some sphingolipid signaling pathways. In C2, insulin resistance and vascular endothelial growth factor (VEGF) signaling pathways were enriched too.
Gene expression trajectories were explored using k-means clustering, with the optimal number of clusters determined via the elbow method. Figure 3A shows two predominant clusters: a first group with increased expression immediately post-marathon that returned to baseline after 24 hours, and a second group exhibiting decreased expression post-marathon followed by recovery at 24 hours. Figure 3B highlights a smaller subset of genes (279 genes) displaying three distinct temporal patterns, from differential expression genes 24 hours after the race compared to baseline.
Finally, in C3, those pathways related to oxidative environment maintenance were among the top significant results (Table 5).
DISCUSSION
Previous studies have investigated the effect of strenuous exercise on human gene expression. The present study was carried out in non-elite athletes to measure the effects of an endurance event such as running a marathon and determine if gene expression levels recovered after 24 hours. The first significant finding was that 60% of the genome was differentially expressed immediately after finishing the marathon compared to baseline levels, indicating that athletes experience major transcriptomic alteration due to endurance exercise. Moreover, 279 genes remained differentially expressed 24 hours after the marathon, suggesting that subjects had not recovered their baseline gene expression levels even after 24 hours of rest.
In general, the results reflect that strenuous exercise can induce an inflammatory response, activate an oxidative environment, and downregulate the immune system. Previous studies have reported that endurance exercise induces an inflammatory environment, which indicates that long periods of strenuous exercise can generally lead to higher levels of inflammatory mediators and therefore may increase the risk of injury and chronic inflammation [26]. The results obtained in this study mirrored these findings for 21 out of 23 inflammatory markers that were significantly differentially expressed at C1 (Table 3).
Immune repression with monocyte activation after marathon running
This figure comprises five subpanels (a–e), illustrating the temporal variation in representative genes across key functional categories: immune markers, chemokines/SOCS, interleukins, Th1/Th2 signaling, and genes with no baseline-level expression at 24 hours.
Subsequently, most immune-related genes were downregulated after the marathon. In Condition 1 (C1), markers involved in the B cell life cycle, specifically those associated with B lymphocyte activation (CD48), pre-B cells (CD19), and LRRC8C and LRRC8D, which are essential for B cell maturation and belong to the T cell activation leucine-rich repeat protein family, were downregulated. Similarly, chemokine receptors expressed by T lymphocytes (CCR4, CCR5, CCR6, and CCR9) showed reduced expression (Table 2). Within the suppressor of the cytokine signaling family, SOCS2 and SOCS7 were downregulated, whereas SOCS3 and SOCS6 were upregulated. In addition, CCL5, a chemokine that mediates chemotaxis for T cells, eosinophils, and basophils and promotes proliferation and activation of natural killer cells, was also downregulated [27].
Furthermore, several interleukins were suppressed in C1, including IL5 and IL32. IL32, a proinflammatory cytokine, contributes to the pathogenesis of autoimmune diseases such as rheumatoid arthritis, while also conferring protection against certain respiratory infections, including tuberculosis [28–29]. Moreover, IL5 regulates genes associated with B cell terminal maturation, as the markers commented above [30].
However, not all immune-related genes were downregulated. The IL6 receptor, a cytokine produced by myocytes during muscle contraction, and genes related to IL10 were significantly overexpressed [31]. This overexpression was not consistent across all IL10-related genes. While ICOS, a protein enhancing IL10 synthesis, and IL10RB were downregulated, IL10RA was upregulated in C1.
Among the overexpressed immune genes, most were specific to monocytes. These included CD14, a surface antigen predominantly expressed on monocytes; CD163, identified as a potential biomarker of monocyte activation in ischemic stroke; and CLEC10A, a gene involved in inflammatory and immune responses (Table 2) [32]. Although CLEC10A mRNA expression occurs in intermediate monocytes, its predominant expression is found in dendritic cells. Additionally, key monocyte chemokine receptors, including CCR2 and CCR1, were overexpressed [33]. IL18 and its receptor IL18R1 were overexpressed, consistent with IL18’s role in enhancing monocyte activation [34]. IL1B, a cytokine primarily produced by monocytes and a key early proinflammatory mediator, was also upregulated, along with its antagonist receptor IL1RN and receptors IL1R1 and IL1R2 [35].
Despite the widespread downregulation observed in other immune components, the upregulation of monocyte markers, specific cytokines, and interleukins indicates that monocyte activation persists following strenuous exercise. Previous studies have shown that intense physical exertion activates monocytes, leading to acute inflammation and hypoxemia. These findings align with the present results, underscoring the pivotal role of monocytes in the immune response to marathon running [36].
Following prolonged exercise, several genes regulating oxidative stress, including the glutathione peroxidase family member GPX7 and superoxide dismutase SOD1, were downregulated. This oxidative environment may represent the molecular link between monocyte chemotaxis and inflammatory pain [37]. Further studies are required to clarify the role and significance of monocytes in this context. Conversely, DNMT3A was overexpressed in C1, consistent with reports showing increased expression in oxidative skeletal muscle after endurance exercise [38]. DNMT3A deficiency is associated with excessive production of reactive oxygen species, a key contributor to muscle dysfunction.
Neuroprotective, chromatin remodeling and mitochondrial adaptations 24 Hours after marathon running
In C3, APP and LRP1 were found overexpressed compared with baseline levels. Both genes are associated with Alzheimer’s disease, and their upregulation following exercise has been reported to exert neuroprotective effects [39]. In particular, LRP1 expression typically declines with disease progression, whereas its enhancement through exercise or pharmacological approaches has been linked to improved amyloid-beta clearance and attenuation of neuroinflammatory responses [40].
EIF3A, a subunit of the eukaryotic translation initiation factor 3 complex, was also overexpressed, suggesting increased regulation of protein synthesis and potential involvement in the cellular response to DNA damage [41]. In contrast, genes related to mitochondrial function, including NDUFS7, PHB2, and ALKBH7, were downregulated, indicating a transient reduction in mitochondrial activity during recovery.
Additionally, the overexpression of ARID1A and AHNAK points to transcriptional and chromatin remodeling responses associated with exercise-induced histone acetylation. Such epigenetic modulation promotes chromatin decompaction and activation of exercise-responsive genes, reflecting ongoing molecular adaptation during the postexercise recovery phase [42].
Immune recovery but persistent mitochondrial alterations 24 Hours after marathon running
The enrichment analysis showed that downregulated genes were primarily associated with Gene Ontology (GO) terms and KEGG pathways related to B and T lymphocytes and other immune system components. This immune suppression may contribute to the increased susceptibility to infection following endurance exercise, as previously reported [43, 44]. These findings are consistent with evidence that endurance athletes experience marked physiological stress during prolonged exercise, leading to transient immunodepression and a higher risk of upper respiratory tract infections due to immune weakening [45]. Conversely, in C3, the results indicated that 24 hours after exercise, immune-related pathways had largely returned to baseline levels, suggesting recovery of immune function during this period.
In C2, enrichment of KEGG pathways related to diabetic cardiomyopathy and atherosclerosis was observed. Previous studies have reported a higher prevalence of coronary artery calcification among long-term marathon athletes, particularly males, which may reflect exercise-induced cardiovascular remodeling [46].
At C3, enriched GO terms and KEGG pathways were associated with mitochondrial activity, energy metabolism, and oxidative stress. Specifically, pathways such as ATP synthesis coupled electron transport, ATP metabolic process, mitochondrial electron transport from NADH to ubiquinone, and oxidative phosphorylation were significantly enriched. These findings indicate that non-elite athletes exhibited reduced mitochondrial energy generation capacity 24 hours after the race compared with their baseline condition.
Identification of two clusters of recovered genes and three clusters of differentially expressed genes 24 hours after marathon completion
Following the enrichment analysis, global transcriptomic profiling revealed two primary gene expression patterns across the three conditions. Genes that returned to baseline 24 hours after the marathon were assigned to one of these clusters, reflecting overall recovery dynamics. Examination of the 279 genes differentially expressed in C3 identified three subclusters with distinct temporal behaviors. Cluster 3, which included APP and LRP1, displayed continued elevation at 24 hours, suggesting prolonged post-exercise activation. Cluster 4 showed upregulation immediately after the marathon, followed by a pronounced decline at 24 hours to levels below baseline. Cluster 5 comprised genes downregulated at the finish line that partially recovered at 24 hours but did not reach baseline levels. These patterns illustrate heterogeneous recovery kinetics among specific genes, highlighting persistent transcriptomic alterations in the early post-exercise period.
Although principal component analysis (PCA) demonstrated clear separation between men and women at the global transcriptomic level (Figure 2), reflecting expected baseline biological differences, our formal interaction analysis did not identify statistically significant sex-dependent transcriptional responses to the marathon after multiple-testing correction. This likely reflects, at least in part, the unbalanced sample composition (42 men and 18 women), which limits power to detect smaller sex effects. The inclusion of sex as a covariate in our principal models mitigates potential confounding, but the present dataset may still underrepresent female-specific molecular signatures. Future studies with larger and more balanced cohorts are warranted to better characterize potential sex-dependent regulatory responses to endurance exercise.
In conclusion, marathon running induces widespread transcriptomic changes, with 9,874 of 14,554 detected genes differentially expressed immediately after the race. We observed immune repression with increased monocyte activity, alongside changes in oxidative stress and inflammatory markers. At 24 hours post-marathon, immune gene expression largely returned to baseline, while mitochondrial function, energy metabolism, and chromatin regulation remained altered. Notably, two genes associated with Alzheimer’s disease, APP and LRP1, remained overexpressed at 24 hours, a pattern previously linked to neuroprotective effects associated with exercise. Clustering analysis identified two groups of recovered genes and three groups of genes with sustained changes, reflecting distinct temporal patterns of recovery.
These results provide a detailed view of transcriptional responses to endurance exercise in non-elite athletes. They highlight gene expression patterns that recover rapidly their baseline expression after endurance exercise and those that remain altered after 24 hours, offering data that could guide personalized training, nutritional strategies, and recovery planning. This study contributes to understanding the temporal dynamics of gene expression following prolonged endurance exercise and identifies targets for future longitudinal studies.
Limitations
While this study provides robust transcriptomic insights into marathon-induced responses, some limitations should be noted. Biochemical parameters such as CRP, creatine kinase, and lactate were not collected due to the practical challenges of multiple blood draws during the event. The dataset was adjusted to account for gender imbalance, though some residual heterogeneity may remain given the higher proportion of men in the cohort. Additionally, key findings were not independently validated by RT-qPCR, which could further support the results. Future studies combining biochemical measures, larger cohorts, and independent validation would further strengthen these observations.
