Genome-Wide Association Study Reveals a Novel Association Between MYBPC3 Gene Polymorphism, Endurance Athlete Status, Aerobic Capacity and Steroid Metabolism

Background The genetic predisposition to elite athletic performance has been a controversial subject due to the underpowered studies and the small effect size of identified genetic variants. The aims of this study were to investigate the association of common single-nucleotide polymorphisms (SNPs) with endurance athlete status in a large cohort of elite European athletes using GWAS approach, followed by replication studies in Russian and Japanese elite athletes and functional validation using metabolomics analysis. Results The association of 476,728 SNPs of Illumina DrugCore Gene chip and endurance athlete status was investigated in 796 European international-level athletes (645 males, 151 females) by comparing allelic frequencies between athletes specialized in sports with high (n = 662) and low/moderate (n = 134) aerobic component. Replication of results was performed by comparing the frequencies of the most significant SNPs between 242 and 168 elite Russian high and low/moderate aerobic athletes, respectively, and between 60 elite Japanese endurance athletes and 406 controls. A meta-analysis has identified rs1052373 (GG homozygotes) in Myosin Binding Protein (MYBPC3; implicated in cardiac hypertrophic myopathy) gene to be associated with endurance athlete status (P = 1.43 × 10−8, odd ratio 2.2). Homozygotes carriers of rs1052373 G allele in Russian athletes had significantly greater VO2max than carriers of the AA + AG (P = 0.005). Subsequent metabolomics analysis revealed several amino acids and lipids associated with rs1052373 G allele (1.82 × 10–05) including the testosterone precursor androstenediol (3beta,17beta) disulfate. Conclusions This is the first report of genome-wide significant SNP and related metabolites associated with elite athlete status. Further investigations of the functional relevance of the identified SNPs and metabolites in relation to enhanced athletic performance are warranted.

Background: The genetic predisposition to elite athletic performance has been a controversial subject due to the underpowered studies and the small effect size of identified genetic variants. The aims of this study were to investigate the association of common single-nucleotide polymorphisms (SNPs) with endurance athlete status in a large cohort of elite European athletes using GWAS approach, followed by replication studies in Russian and Japanese elite athletes and functional validation using metabolomics analysis.
Results: The association of 476,728 SNPs of Illumina DrugCore Gene chip and endurance athlete status was investigated in 796 European international-level athletes (645 males, 151 females) by comparing allelic frequencies between athletes specialized in sports with high (n = 662) and low/moderate (n = 134) aerobic component. Replication of results was performed by comparing the frequencies of the most significant SNPs between 242 and 168 elite Russian high and low/moderate aerobic athletes, respectively, and between 60 elite Japanese endurance athletes and 406 controls. A meta-analysis has identified rs1052373 (GG homozygotes) in Myosin Binding Protein (MYBPC3; implicated in cardiac hypertrophic myopathy) gene to be associated with endurance athlete status (P = 1.43 × 10 −8 , odd ratio 2.2). Homozygotes carriers of rs1052373 G allele in Russian athletes had significantly greater VO 2max than carriers of the AA + AG (P = 0.005). Subsequent metabolomics analysis revealed several amino acids and lipids associated with rs1052373 G allele (1.82 × 10 −05 ) including the testosterone precursor androstenediol (3beta,17beta) disulfate.
Conclusions: This is the first report of genome-wide significant SNP and related metabolites associated with elite athlete status. Further investigations of the functional relevance of the identified SNPs and metabolites in relation to enhanced athletic performance are warranted.
Keywords: GWAS, SNP, metabolomics, metabolites, elite athletes, endurance BACKGROUND Elite athletic performance is a multi-factorial trait with input from both genetic and environmental factors. The superior performance of elite athletes has been historically considered an outcome of a special talent shaped by intensive training. The talent is now believed to be a product of additive genetic components predisposing the athlete to endurance, speed, strength, flexibility and coordination trainability under the control of strong environmental cues including exercise and nutrition. In this model, the genetic predisposition together with ability to respond to training are the keys to the superior physical performance of elite athletes (Georgiades et al., 2017).
Sports can be classified according to the type and intensity of the exercise required to perform during competition. The percentage of maximal oxygen uptake (VO 2max ) is a detrimental factor in the categorization of endurance sports, as it reflects the maximal cardiac output, the oxygen transport capacity, and the blood volume (Bergh et al., 2000). Accordingly, sports can be divided into sport events with low, moderate and high aerobic (dynamic) component (Mitchell et al., 2005). Similarly, the percent of maximal voluntary contraction (MVC), which reflects the greatest amount of tension a muscle can generate and hold, is used to classify sports into sporting disciplines with low, moderate and high power component (Mitchell et al., 2005).
Classical twin and family genetic studies have suggested that VO 2max is up to 94% inherited (Bouchard et al., 1998;Peeters et al., 2009). Genome-wide association studies (GWAS) in athletes versus non-athletes have uncovered many new loci in association with VO 2max (Rankinen et al., 2010; Abbreviations: ACP2, acid phosphatase 2, Lysosomal; ADLQ, anti-doping laboratories in Qatar; FDR, false discovery rate; FMSI, Laboratorio Antidoping, Federazione Medico Sportiva Italiana; GVS, genome variation server; GWAS, genome-wide association studies; HESI-II, high resolution/accurate mass spectrometer interfaced with a heated electrospray ionization; MADD, MAP kinase activating death domain; MAF, minor allele frequency; MVC, maximal voluntary contraction; MYBPC3, myosin binding protein C, cardiac; NR1H3, nuclear receptor subfamily 1 group H member 3; OR, odds ratio; Spi-1, Spi-1 proto-oncogene; UPLC, ultra-performance liquid chromatography; VO 2max , maximal oxygen uptake. Bouchard et al., 2011) and elite endurance performance . A more recent review of genetic predisposition to elite athletic endurance has highlighted 100 endurance variants . However, despite some initial evidence suggesting identification of genetic variants in GWAS studies, further studies did not replicate/validate these findings hindered by a small sample size and complex phenotype . One of the first GWAS in athletes using 143 K single-nucleotide polymorphisms (SNPs) and subsequent meta-analysis of 45 promising genetic markers in 1,520 endurance athletes and 2,760 controls has revealed only one statistically significant marker (rs558129 at GALNTL6) associated with endurance status in world class athletes, but not at genome wide level of significance (Rankinen et al., 2016). Therefore, the genetic predisposition to endurance traits remains unclear, largely due to the relatively underpowered elite athletes' cohorts. Recently, a polymorphism in human homeostatic iron regulator protein was found to be associated with elite endurance athlete status and aerobic capacity in Russian athletes (Semenova et al., 2020).
Metabolomics analysis has presented a novel tool to validate genomics data by providing an intermediate phenotype (metabolites) in association with the identified genetic variants (Kastenmuller et al., 2015;Tanaka et al., 2016). Pilot metabolomics studies have revealed differences in the metabolic signature of moderate and high endurance elite athletes, such as steroid biosynthesis, fatty acid metabolism, oxidative stress and energy-related molecular pathways (Al-Khelaifi et al., 2018, 2019a. Recently, a study investigating metabolic GWAS of elite athletes showed novel genetically influenced metabolites associated with athletic performance. These included two novel genetic loci in FOLH1 and VNN1 in association with N-acetylaspartyl-glutamate and linoleoyl ethanolamide, respectively, and one novel locus linking genetic variant in SULT2A1 and androstenediol (3alpha, 17alpha) monosulfate in endurance athletes (Al-Khelaifi et al., 2019b).
In this study, we aimed to investigate the association of multiple SNPs and endurance athlete status in a relatively large cohort of European elite athletes specialized in sports with high and low/moderate aerobic component using GWAS approach and replicate our findings in elite Russian and Japanese athletes. We also aimed to perform functional validation using VO 2max testing and metabolomics analysis by identifying metabolites that are associated with significant endurance-related SNPs.

Genome-Wide Association Study
Athletes from the discovery cohort were classified into different groups of sports following previously published sports classification criteria (Mitchell et al., 2005), as shown in Table 1.
The principle component analysis (PCA) of the genotyping data revealed no influence of sport disciplines ( Figure 1A) or training modality (i.e., sports with low/moderate versus high aerobic component) ( Figure 1B) on genotype distribution. Following quality control data processing, genotyping of 341385 SNPs in 796 European elite athletes revealed several variants associated with endurance athlete status, but none reached GWAS level of significance. Table 2 shows top SNPs (P < 5 × 10 −5 ) with their odd ratios (OR) in relation to elite athletic endurance, location according to function genome variation server (GVS), gene name and minor allele frequency (MAF) in sports with high and low/moderate aerobic component. MAF in non-elite athletes from 1,000 genome project were used as a reference. Figure 1 shows Manhattan (C) and quartile-quartile (QQ) plots (D) of GWAS hits associated with endurance.

Replication of Endurance SNPs in Russian and Japanese Elite Athlete Cohorts
Replication of results was performed by comparing the frequencies of the most significant SNPs (P < 10 −5 ) in 242 elite Russian high and 168 low/moderate aerobic athletes, and in 60 elite Japanese endurance athletes and 406 controls. Out of the 9 top SNPs identified form the GWAS discovery stage, the rs1052373 (MYBPC3) and rs7120118 (NR1H3) showed significant association with endurance in Russian and Japanese (P < 0.05). However, the association was driven by a dominant model since results of this analysis showed over representation for rs1052373 GG and rs7120118 TT genotypes in the high endurance group. A subsequent meta-analysis has confirmed the over representation of the rs1052373 GG and rs7120118 TT genotypes in high endurance sports at genome-wide and Bonferroni levels of significance (1.43 × 10 −8 and 1.66 × 10 −7 , respectively) ( Table 3). The combined analysis showed no evidence of heterogeneity and direction of association was similar in all three cohorts.
The regional association plot for the rs1052373 G allele in MYBPC3 gene revealed a number of SNPs in the same LD block in association with high endurance including the rs7120118 T allele in NR1H3 gene (Figure 2).
To validate the potential functionality of the identified GWAS SNPs, association of the identified two SNPs (rs1052373 G and rs7120118 T alleles) with VO 2max was investigated in a subgroup of the Russian replication cohort in which VO 2max data was available. This included 32 elite Russian long-distance athletes [19 biathletes, 13 cross-country skiers; 17 females, age 23.5 (3.5) years; 15 males, age 21.3 (4.1) years]. The rs1052373 GG carriers had significantly greater VO 2max than carriers of the AA + AG (P = 0.005 adjusted for sex). Similarly, rs7120118 TT carriers showed a trend of higher VO 2max than carriers of the CC + CT (P = 0.053 adjusted for sex).
For further validation of the potential functionality of the identified GWAS SNPs, metabolomics of 750 metabolites was carried out in a subset of the discovery cohort (n = 490) and enriched metabolic pathways associated with the rs1052373 G allele and rs7120118 T alleles were determined ( Table 4). Among the metabolic pathways associated with rs56330321 and rs7120118, various lipids and amino acids were significantly altered by their genotypes. However, only 5alpha-androstan-3alpha,17alpha-diol disulfate reached Bonferroni level of significance (Table 4), exhibiting higher levels in rs1052373 GG and rs7120118 TT carriers compared to AA + AG and CC + TC carriers, respectively (Figure 3).

DISCUSSION
Genetic predisposition into cardiorespiratory fitness and response to exercise training has been previously described (Lortie et al., 1982;Prud'homme et al., 1984;Hamel et al., 1986;Bouchard et al., 1994Bouchard et al., , 1998Bouchard et al., , 1999. Since endurance performance sports are characterized by increased cardiorespiratory capacity, genetic predisposition into elite endurance performance is also expected to be genetically influenced (Guth and Roth, 2013). However, genetic studies of elite athletic endurance showed inconsistent results (Guth and Roth, 2013;Ahmetov and Fedotovskaya, 2015;Pitsiladis et al., 2016;Wang et al., 2016). The aims of this study were to carry out the largest GWAS study of elite European athletes to date using a unique SNP microarray that is enriched with genes involved in different metabolic pathways with direct influence on various physiological pathways characteristic of elite athletes. GWAS results have revealed a number of novel SNPs associated with endurance but none reached the GWAS level of significance. Replication of the top identified SNP associations in two independent cohorts of elite athletes from Russia and Japan has confirmed the association of rs7120118 and rs1052373 with endurance athlete status. Subsequent meta-analysis of the three cohorts has revealed for the first time that both SNPs were associated with endurance athlete status at genome-wide and Bonferroni level of significance, respectively. Functional validation has revealed the association of the two SNPs with increased Vo 2max and levels of the testosterone precursor 5alpha-androstan-3alpha,17alpha-diol disulfate. Distribution of elite athletes in various categories based on sport type-associated peak dynamic (maximal oxygen uptake percentage; VO 2max ) and peak static (maximal voluntary muscle contraction percentage; MVC) components achieved during competition as described previously (Mitchell et al., 2005).
FIGURE 1 | GWAS data quality control. PCA shows no difference in the genotype distribution among sport disciplines (A) or between groups (sports with low/moderate versus high aerobic component) (B) Manhattan (arrow indicates significant SNPs) (C) and Quantile-quantile (no evidence of genomic inflation, lambda GC = 1.006) (D) plots illustrating GWAS results in association with endurance.
The top identified GWAS significant SNP (rs1052373) is located within MYBPC3 gene. MYBPC3 codes for a myosin-associated protein expressed in the cross-bridgebearing zone (C region) of A bands in striated muscle. The phosphorylation of MYBPC3 protein modulates cardiac contraction (Moss et al., 2015). Mutations in MYBPC3 were previously associated with a lower super-relaxed state in patients with hypertrophic cardiomyopathy (HCM) (McNamara et al., 2017). Intense exercise can trigger heart remodeling to compensate for the elevations in blood pressure or volume by increasing muscle mass. Hence, hearts of the endurance athletes typically exhibit an eccentric cardiac hypertrophy with increased cavity dimension and wall thickness (Pelliccia et al., 1991;Hedman et al., 2015), which is influenced by the type of sport performed (Pelliccia, 1996;Pelliccia et al., 1999;Maron and Pelliccia, 2006). As a result, the endurance-trained  heart can deliver a large maximal systolic volume (35% larger than untrained heart) in order to produce a large cardiac output (Ogawa et al., 1992;Pelliccia et al., 1999). Since carriers of the GG allele exhibit a benign phenotype of HCM according to NIH's ClinVar database (Landrum et al., 2018), the mild phenotype may be enhancing exercise-triggered physiological adaptations.
The seemingly dominant effect of rs1052373 GG on increased VO 2max and endurance may support this added advantage although more studies are needed to confirm this finding. These adaptations, however, might be associated with a greater risk of cardiovascular disease. Indeed, we have recently shown that endurance athletes with high cardiovascular demand (higher blood pressure and stroke volume) show metabolic signature consistent with higher risk of cardiovascular disease (Al-Khelaifi et al., 2019a). When investigating the expression quantitative trait loci (eQTLs) associated with rs1052373, a number of genes was identified including SPI1, MYBPC3, MADD, ACP2 and NR1H3 (Ray et al., 1990;Tang and Chu, 2002;Mannan et al., 2004;Wu et al., 2012;Carrier et al., 2015;Theofilopoulos and Arenas, 2015). Interestingly, eQTL (GTEx) showed that rs1052373 polymorphism is associated with expression level of MADD and ACP2 in heart, but not MYBPC3. Since MAP kinase plays an important role of cardiac hypertrophy (Zhang et al., 2003), the association between rs1052373 polymorphism and VO 2max and endurance may also be explained by MADD expression, although this needs further validatoin. Information related to function and associated diseases with these genes are summarized in Supplementary Table S1. The other significant association was between rs7120118 TT carriers and high endurance. Rs7120118 is located in NR1H3 gene that codes for a nuclear receptor regulating macrophage function, lipid homeostasis and inflammation. NR1H3, also known as liver X Receptor Alpha (LXRA), plays an important role in the regulation of cholesterol homeostasis including adrenal steroidogenesis (Repa et al., 2002;Cummins et al., 2006). The association of rs7120118 with high endurance could be reflecting the high linkage disequilibrium (r 2 = 0.89, P < 0.0001) between rs7120118 TT and the potentially functional rs1052373 GG. It could, however, be related to increased synthesis of the testosterone precursor 5alphaandrostan-3alpha,17alpha-diol disulfate since NR1H3 regulates hypothalamo-pituitary-adrenal steroidogenesis (Handa et al., 2011). Indeed, we have previously shown that high-endurance athletes exhibit elevated levels of several sex hormone steroids involved in testosterone synthesis including 5alpha-androstan-3alpha,17alpha-diol disulfate (Al-Khelaifi et al., 2018) with implication on improving performance due to enhanced glucose metabolism and protein synthesis in the muscle (Sato et al., 2008). The functional relevance of these associations remains to be further validated.
Study limitations: The lack of information about participants and the heterogeneity of their sport groups were major limitations of this study. To overcome these limitations and to increase the power of the study, genotyping was compared between athletes who belong to high endurance versus moderate endurance performance sports instead of power versus endurance due to the overlap between the two classes as per Mitchell's categorization (Mitchell et al., 2005). Other limitations included using add-on replication studies (Russian and Japanese cohorts) rather than using a carefully designed replication. However, differences were confirmed in each study separately and the subsequent metaanalysis confirmed the significance of the association of the two SNPs with endurance.

CONCLUSION
This study reports the first GWAS significant SNP (rs1052373) in MYBPC3 in association with endurance athlete status with a direct relevance to cardiac hypertrophy and contraction. The SNP is associated with increased VO 2max and elevated levels of the testosterone precursor androstenediol (3beta,17beta) disulfate, both phenotypes that potentially contribute to the superior performance of endurance athletes. This study also identifies a second SNP (rs7120118) associated with endurance at Bonferroni level of significance in NR1H3. This SNP could be either working independently of rs1052373 through influencing steroidogenesis or could be acting as a marker of rs1052373. Further investigations of the functional relevance of the identified SNPs and associated metabolites in relation to enhanced athletic performance are warranted.

METHODS
The aim of this study is to investigate the genetic predisposition to elite athletic endurance through conducting the largest GWAS in elite athletes to date, followed by functional validation through aerobic capacity testing and metabolomics analysis to shed light on the underlying mechanisms of genetic associations.

Discovery Study
Seven hundred and ninety six consented European internationallevel athletes (645 males, 151 females) from different sports disciplines who participated in national or international sports events and tested negative for doping substances at anti-doping laboratories in Qatar (ADLQ) and Italy (FMSI) were included in this study. No other information of participants was available due to the strict anonymization process undertaken by the anti-doping laboratories. This study was performed in line with the World Medical Association Declaration of Helsinki -Ethical Principles for Medical Research Involving Human Subjects. All protocols were approved by the Institutional Research Board of ADLQ (F2014000009). Athletes were dichotomized into groups with different aerobic (dynamic) and power (static) components (Table 1) based on their sport types as described previously (Mitchell et al., 2005). Table 1 further lists the number of participants based on various analyses as per sport type in each class/group and their genders.
FIGURE 2 | Regional association plot for the region around rs1052373. The colors correspond to different LD thresholds, where LD is computed between the sentinel SNP (lowest P-value, colored in blue) and all SNPs. Shapes of markers correspond to their functionality as described in the legend.
Aerobic Capacity Testing VO 2max in biathletes and cross-country skiers was determined using an incremental test to exhaustion on a treadmill HP Cosmos (Germany). The initial speed was 7 km/h, the increment was 0.1 km/h every 10 s.VO 2max was determined breath by breath using a MetaMax 3B-R2 gas analysis system.VO 2max was recorded as the highest mean value observed over a 30 s period.

Discovery Study
DNA was extracted from leukocytes (venous blood) samples from all participants using DNeasy Blood & Tissue kit (Qiagen) following manufacturer's instructions. The concentration and the quality of DNA were assessed using the Nanodrop (Thermo Fisher) and Qubit Fluorometer (Invitrogen) to ensure sufficient amount and quality of DNA were obtained for genotyping. Illumina Drug Core array-24 BeadChips was chosen for the genotyping of 476,728 SNPs in the 796 European elite athletes collected for Anti-Doping analysis (discovery cohort). This array contains over 240,000 highly-informative genome-wide tag SNPs and a novel ∼200,000 custom marker set designed to support studies of drug target validation and treatment response. The assay required 200 ng of DNA sample as input with a concentration of at least 50 ng/µl. All further procedures were performed according to the instructions of Infinium HD Assay according to manufacturer's instructions. Briefly, 4 µl of obtained DNA was mixed with Illumina amplification reagents and incubated overnight at 37 o C in hybridization oven. On the second day, enzymatic reagents were used to fragment the amplified DNA then precipitated by centrifugation. Subsequently, re-suspended pellet was loaded in the beadchip then incubated overnight at 48 o C in hybridization oven. On third day, beadchips underwent enzymatic base extension and fluorescent staining. Lastly, after coating, the beadchips were imaged using iScan.

Replication Studies
Molecular genetic analysis in Russian cohorts was performed with DNA samples obtained from leukocytes (venous blood). Four ml of venous blood were collected in tubes containing EDTA (Vacuette EDTA tubes, Greiner Bio-One, Austria). Blood samples were transported to the laboratory at 4 • C and DNA was extracted on the same day. DNA extraction and purification were performed using a commercial kit according to the manufacturer's instructions (Technoclon, Russia) and included chemical lysis, selective DNA binding on silica spin columns and ethanol washing. Extracted DNA quality was assessed by agarose gel electrophoresis at this step. HumanOmni1-Quad BeadChips (Illumina Inc, United States) were used for genotyping of 1,140,419 SNPs in athletes and controls. The assay required 200 ng of DNA sample as input with a concentration of at least 50 ng/µl. Exact concentrations of DNA in each sample were measured using a Qubit Fluorometer (Invitrogen, United States). All further procedures were performed according to the instructions of Infinium HD Assay. For the second replication study, total DNA was isolated from saliva or venous blood using Oragene·DNA Collection Kits (DNA genotek, Ontario, Canada) or QIAamp DNA blood Maxi Kit (QIAGEN, Hilden, Germany), respectively. The total DNA content was measured using a NanoDrop 8000 spectrophotometer (Thermo Fisher Scientific, MA, United States). Subsequently, DNA samples were adjusted to a concentration of 50 ng/µL with TE buffer and were stored at 4 • C. Total DNA samples were genotyped for more than 700,000 markers using the Illumina R HumanOmniExpress Beadchip.

Data Extraction and SNP Identification
Raw data was extracted, peak-identified and QC processed using Illumina iScan hardware and software. These systems are built on a web-service platform utilizing Microsoft's NET technologies, which run on high-performance application servers and fiberchannel storage arrays in clusters to provide active failover and load-balancing.

Metabolomics
Screening of serum metabolites was performed in 490 elite athletes (Supplementary Table S2) using protocols established at Metabolon, Durham, NC, United States. The platform utilizes Waters ACQUITY ultra-performance liquid chromatography (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. Detailed protocol and QC measures were previously published (Evans et al., 2009;Al-Khelaifi et al., 2018).

Statistical Analysis
Following genotyping using Illumina's Drug Core SNP array, analysis was performed using Plink v1.9. Quality control measures were applied to the genotype data set to exclude samples with low genotype call rate or excess heterozygosity. Accordingly, SNPs with a genotype call rate <98%, minor allele frequency <1%, or deviating from Hardy-Weinberg equilibrium (P < 10 −6 ) were excluded. After filtering the data with the above criteria, 341,385 SNPs were used in analysis. Population background was determined using principal component analysis (PCA) in comparision to samples from HapMap project and only samples with European ancestry were included in the analysis. The analysis in European and Russian cohorts was performed using linear or logistic regression models. A model incorporating sports grouped by training modalities (i.e., sports with high versus low/moderate aerobic component) was used for the discovery cohort after incorporating gender and PCA components 1, 2, 3 & 4 as covariates in the model. A stringent Bonferroni level of significance of P ≤ 0.05/341385 = 1.46 × 10 −7 was used to define significant associations. To perform the metaanalysis, the Cochrane Review Manager version 5.3 was used. Random and fixed effect models were applied. The heterogeneity degree between the studies was assessed with the I 2 statistics. Associations between SNPs and metabolite levels were computed using lm function in R (version 3.3.1) while correcting for gender, hemolysis and PCA. An additive inheritance model was used (SNPs were coded as 0,1,2 according to their genotype group. Pathway enrichment analyses were carried out using Chi square tests to identify pathways with enriched metabolites ranked by P-value from the linear model since Bonferroni level of significance was not observed.

ETHICS STATEMENT
This study was performed in accordance with the World Medical Association Declaration of Helsinki. All protocols were approved by the Institutional Research Board of anti-doping lab Qatar (F2014000009). The patients/participants provided their written informed consent to participate in this study.