Joint Identification of Genetic Variants for Physical Activity in Korean Population

There has been limited research on genome-wide association with physical activity (PA). This study ascertained genetic associations between PA and 344,893 single nucleotide polymorphism (SNP) markers in 8842 Korean samples. PA data were obtained from a validated questionnaire that included information on PA intensity and duration. Metabolic equivalent of tasks were calculated to estimate the total daily PA level for each individual. In addition to single- and multiple-SNP association tests, a pathway enrichment analysis was performed to identify the biological significance of SNP markers. Although no significant SNP was found at genome-wide significance level via single-SNP association tests, 59 genetic variants mapped to 76 genes were identified via a multiple SNP approach using a bootstrap selection stability measure. Pathway analysis for these 59 variants showed that maturity onset diabetes of the young (MODY) was enriched. Joint identification of SNPs could enable the identification of multiple SNPs with good predictive power for PA and a pathway enriched for PA.


Introduction
Physical activity (PA) is any bodily movement produced by skeletal muscles that requires energy expenditure [1]. It includes exercise as well as work and recreational activities which involve bodily movements. PA can be quantified by metabolic equivalent task (MET) intensity [1,2], and plays an important role in the morbidity, mortality, and health care costs of obesity and related chronic diseases [3].
The behavior of PA may be determined by genetic factors [4][5][6]. In particular, between 48% and 71% of the variability in adult exercise behavior can be explained by genetic factors [6]. Indeed, evidence from twin and family studies have suggested that genetic factors contribute to the propensity of being sedentary [5]. A novel approach has been proposed to identify genetic variants that are related to leisure-time exercise behavior, by conducting Genome-Wide Association (GWA) analyses using logistic regression to find genes associated with exercisers and non-exercisers [7]. Recently, a genome-wide study with quantitative PA as a phenotype was performed for the Korean population [8]. This study revealed how to define phenotypes of PA in genetic association studies, together with appropriate statistical methods for their analysis [8].
In reporting genetic variants associated with a trait or disease, most traditional GWA studies adopted a single-marker approach that identifies single genetic factors one by one. However, this method is inefficient in predicting joint effects of multiple genetic variants on the common complex trait [9,10]. A multiple-marker approach is the preferred alternative for their joint identification [11]. However, multiple linear and logistic regression models are often ill-defined in GWA studies when the number of predictor variables is larger than the sample size. In addition, collinearity often occurs between predictor variables due to linkage disequilibrium among single nucleotide polymorphisms (SNPs).
To identify multiple genetic variants for common complex traits or diseases, an elastic-net (EN) regularization method had been proposed by Cho et al. [11], along with some consistency measures based on bootstrap sampling. The EN regularization method was originally introduced for model fitting and variable selection in ill-defined multiple regressions. It has been applied to GWA studies [12,13] and provided a better prediction than those based on ordinary regression models [14,15] when variables are larger than sample sizes and multicollinearity problems may exist.
The present study aimed to find genetic variants influencing PA in the Korean population. Single-SNP association tests were initially conducted to assess genetic associations between daily PA and various SNP markers. We then performed multiple SNP analysis with EN regularization to determine genetic variants associated with PA. The SNPs identified provide novel biological evidence to understand the genetics of PA through pathway enrichment analysis.

Physical Activity Levels
Overall, the average daily PA level was 1332 (SD 871) MET·min for the Korean participants. Men (mean 1367, SD 887 MET·min·day −1 ) appeared to be slightly more active than women (mean 1300, SD 856 MET·min·day −1 ). The mean PA level of the Ansung cohort (1678, SD 1046 MET·min·day −1 ) was higher than that of the Ansan cohort (1038, SD 534 MET·min·day −1 ), suggesting that people in the rural community tended to be more active than their city counterparts. Figure 1 shows the box plots of self-reported PA by age groups (40)(41)(42)(43)(44)(45)(46)(47)(48)(49), 50-54, 55-59, 60-64, 65+). Although the median PA levels of these six age groups were similar, the PA distributions exhibited substantial variations. In particular, the majority of younger participants appeared to sustain low PA levels with the exception of a few outliers. For older participants especially those over 65 years, their PA levels varied considerably between individuals, as evident from the wide interquartile ranges.

Individual Single Nucleotide Polymorphism (SNP)-Based Association Analysis
Single-marker association analysis was performed for individual SNP with sex, age, area, and body mass index as covariates. Table 1 presents the results of the single-SNP association tests. The first six columns give the SNP information and the remaining columns summarize the regression results. Figure 2 further shows the Manhattan plot of 344,893 SNPs, where the y-axis represents the log-transformed p-value and the x-axis represents the chromosomes. The horizontal solid line indicates p-value = 10 −5 . Although these SNPs did not achieve the genome-wide level significance, 41 genetic variants (listed in Table 1) emerged among the 344,893 SNPs to have some evidence of association with PA under p-value < 10 −4 , and they were mapped to 27 genes.   and body mass index included as covariates; c p-values from single-marker association test indicates p-value = 10 −5 . Although these SNPs did not achieve the genome-wide level significance, 41 genetic variants (listed in Table 1) emerged among the 344,893 SNPs to have some evidence of association with PA under p-value < 10 −4 , and they were mapped to 27 genes.

Multiple SNP-Based Association Analysis
The multi-stage procedure was applied to identify multiple causal SNPs. After performing single-SNP association tests in the first stage, we chose top 1000, top 2000, top 3000, and top 4000 SNPs that exhibit the strongest individual associations with PA. During the next stage, 639, 1248, 1760, and 2239 from the respective top SNP groups were jointly identified as PA-related genetic variants by the EN regularization method. The elastic-net allows correlation among predictors in variable selection, so there can be different selection results according to pre-screened datasets [11]. At the validation stage, these jointly identified SNPs were further evaluated using the bootstrap selection stability (BSS) measure. Since the lists of SNPs were different depending on the number of pre-screened SNPs, we focused on the common 457 SNPs that were simultaneously identified from all four groups. These commonly selected variants have higher BSS than the variants which are chosen only one pre-defined dataset [11]. Table 2 lists the final 59 variants selected with BSS ≥0.95 and mapped to 76 known genes. The p-values were calculated using a multiple linear regression model with adjustment for sex, age, area, and body mass index. Finally, pathway analysis found one pathway enriched in these 76 genes: Maturity onset diabetes of the young (MODY). This pathway includes GCK and HES1 genes and has a p-value of 0.076.  The predictive power of the 2239 SNPs identified from the top 4000 SNPs was next investigated. SNPs were ranked in order form smallest to largest and selected using a given BSS cut-off value, 0.95. Then, a multiple regression model with the selected SNPs was fitted to compute the corresponding adjusted R 2 value. Figure 3 compares the adjusted R 2 values versus number of SNPs between the multiple SNP analysis and the single-marker approach. It is clear that the predictive power increases with the number SNP for both approaches, but the multi-stage approach always performs better than the single-marker approach for prediction purpose. This shows that multi-stage approach using a BSS cut-off value provides a better explanation of phenotype than the single marker approach.

Discussion
The present study investigated genetic factors associated with PA for the Korean population, by performing large-scale GWA through single-SNP analysis and multiple SNP analysis via the EN regularization method. Single-SNP association tests are appropriate to determine individual associations between each SNP and the trait or phenotype. However, if the purpose is to predict the phenotype, then the joint identification of genetic factors would be powerful and provide a better prediction of the trait when multiple genetic factors exist for a common complex trait. In the presence of multicollinearity due to linkage disequilibrium among SNPs, EN regularization with BSS offers more accurate identification of multiple SNPs than ordinary multiple regression analysis.
Our single-marker analysis results showed that, although the most significant SNP did not attain the genome-wide significance level (rs7023003, p-value = 4.67 × 10 −6 ), 41 SNPs did exhibit some evidence of association with PA at a less stringent significance level (p-value < 10 −4 ). Among the 27 genes identified, CDKAL1, CDK5 regulatory subunit associated protein 1-like 1, is one of the tailanchored protein family and associated with a type 2 diabetes susceptibility gene responsible for tRNALys modification [16]. CDKAL1 explores the TCR40/Get3 assisted pathway for insertion of its C-terminal transmembrane domain into the endoplasmic reticulum [16]. It has been reported that this gene can influence insulin response and risk of type 2 diabetes [17].
TFPI2, tissue factor pathway inhibitor 2, is found to be a tumor suppressor gene and frequently inactivated through promoter methylation in several kinds of tumors [18]. In particular TFPI2 methylation plays a key role in the diagnosis of colorectal cancer and has been demonstrated to exist in colorectal cancer patients' sera [19].
CCHCR1, also called coiled-coil alpha-helical rod protein 1, is a candidate gene for Psoriasis which is one kinds of chronic inflammatory skin disorder [20]. CCHCR1 gene is found to be expressed in psoriatic lesions compared to normal healthy skin or other hyper proliferative skin disorders [21,22]. CCHCR1 has been demonstrated to be involved in steroidogenesis of the skin [20].
Another gene, RHOBTB1, encodes Rho-related BTB domain-containing protein 1 [23]. The protein encoded by this gene belongs to the Rho family of the small GTPase superfamily. It has a GTPase domain, a proline-rich region, a tandem of 2 BTB (broad complex, tramtrack, and bric-a-brac) domains, and a conserved C-terminal region. The protein plays a role in small GTPase-mediated signal transduction and the organization of the actin filament system. Alternate transcriptional splice variants have been characterized. It is known that the gene is highly expressed in skeletal muscle, stomach, placenta, kidney, and testis.
The single-marker analysis also identified ASTN2. It encodes a protein that is expressed in the brain and may function in neuronal migration, based on functional studies of the related astrotactin 1 gene in human and mouse. A deletion at this locus was shown to be associated with schizophrenia [24]. Multiple transcript variants encoding different proteins have been identified in this locus.
Through the multi-stage approach, subsets of multiple SNPs were jointly identified through the EN regularization method. Among the 457 common SNPs found, 59 SNPs with BSS values exceeding 0.95 were mapped to 76 genes. Of these genes, ADAM12 encodes ADAM metallopeptidase domain 12 and is a disintegrin and metalloproteases family member [25]. ADAM12 is also a multidomain type I transmembrane protein that functions both in normal physiology and in diseases [25]. CCNI gene, cyclin I, is also identified. CCNI is known to be expressed in human forebrain cortex [26]. CCNI is also presented in skeletal muscle, heart, and brain and expressed constantly during cell cycle progression [26]. Moreover, PTPN5 encodes protein tyrosine phosphatase, non-receptor type 5 and involves in regulating the occurrence of abnormal stress responses underlying depression-related disorders [27]. The basal levels of PTPN5 expression in the dorsal hippocampus determine an individual's susceptibility for developing stress-related cognitive and morphological changes [27]. Another gene found, NRXN3 (neurexin 3), is a member of the neuroxines protein family that acts in the vertebrate nervous system as cell adhesion molecules and receptors [28]. The mutations of NRXN3 are related to alcohol and nicotine dependence in patients who suffer from schizophrenia [29,30].
Pathway analysis of the 59 SNPs led to the identification of MODY. Its inherence is responsible for non-insulin-dependent diabetes typically diagnosed among young people, especially under 25 years of age [31,32]. MODY is often referred to as monogenic diabetes, which is thought to be different from type 1 and type 2 diabetes [33,34]. MODY is also known to include the GCK gene and the HES1 gene. The gene GCK, glucokinase, is one of the frequently causing MODY genes and accounts for approximately 35% of cases [33]. HES1, Hairy enhancer of split 1, is one of the highly conversed family of Hairy related basic helix-loop-helix (bHLH) proteins [35]. HES1 is also known to be involved specifically in Notch1 signaling in neural cells and in bone narrow [35].
The main limitation of our study concerns the self-reporting nature of PA and the questions used were not sufficiently detailed. In addition to possible recall error by the participants, the phenotypes of total PA level may not be well defined. PA was classified by five categories based on MET intensities with average MET assigned to each category. The estimated total PA level may dilute MET intensities and lead to a low power for the genetic association study. Furthermore, it is not feasible to compare our results with previous studies in the absence of similar genotyped SNPs. Consequently, we adopted the bootstrap sampling scheme to obtain replication data sets with appropriate selection stability measure to confirm the findings [11].

Subjects
Study subjects were selected from an ongoing population-based cohort, as part of the Korean Genome and Epidemiology Study (KoGES). Participants were recruited from residents in two cities (Ansung and Ansan) in Gyeonggi-do, Korea. We enrolled 10,038 males and females between 2001-2002 for a baseline study, whose demographics have been reported [8]. The Korean Genome and Epidemiology Study was launched in 2007, whereby over 10,000 subjects were recruited from two community-based cohorts: the rural Ansung and urban Ansan cohorts in the Gyeonggi of Korea. The initial samples included 5018 and 5020 participants aged 40 to 69 years from the two cohorts, respectively [36]. Table 3 summarizes the demographic characteristics of the participants. There were more female participants in Ansung than Ansan but the Ansung cohort was on average older than the Ansan cohort, reflecting the differences between rural and urban areas. This study obtained approval from the appropriate institutional review boards of each participating institution, and written informed consent was obtained from all participants.

Physical Activity Information
Information on intensity and duration of daily PA was obtained from each participant using a structured questionnaire that included five components on PA: stable (lying down except sleeping), sitting (e.g., during typing, playing cards, driving, office work, attending a class), low intensity (e.g., walking, doing the laundry, cleaning, leisure time ping pong), medium intensity (e.g., walking, carpentering, regular exercise, badminton, swimming, tennis), high intensity (e.g., sports competition, climbing, running, logging, farming). For each type of PA, its duration was measured in minutes.
Since each question contained multiple PAs with varying MET intensities, the average MET intensity was assigned. The total amount of daily PA was then calculated by summing across the products of the average MET and the corresponding durations [37].

Genotypes
Genomic DNA samples were isolated from peripheral blood drawn from the participants. The majority of genomic DNAs were genotyped on the Affymetrix Genome-Wide Human SNP array 5.0 containing 500,568 SNPs. From the total of 10,038 participants, 10,004 available samples were genotyped [36] by implementing Bayesian robust linear modeling with the Mahalanobis distance (BRLMM) algorithm, and standard quality control procedures were adopted. Samples with high missing genotype call rate (>4%, n = 401), high heterozygosity (>30%, n = 11), gender inconsistencies (n = 41), and those obtained from individuals who had developed any form of cancer (n = 101), were excluded from subsequent analyses, along with related or identical individuals whose computed average pairwise identity-by-state value was higher than that estimated from first-degree relatives of Korean sib-pair samples (>0.80, n = 601). Samples whose genotype-inferred sex disagreed with clinical records were re-tested for sex confirmation using the SNaPshot Multiplex System (Applied Biosystems, Life Technologies, Carlsbad, CA, USA). Markers with high missing gene call rate (>5%), low Minor Allele Frequency (MAF) (<0.01) and significant deviation from Hardy-Weinberg equilibrium (p-value < 1 × 10 −6 ) were excluded, leaving a total of 352,228 markers to be examined among 8842 individuals. For multiple SNP analysis, reduced information was common due to missing values. To get complete genotype data, we imputed missing genotypes using the Fastphase software (University of Washington, Seattle, WA, USA) [38]. As a result, a total of 344,893 SNP markers with chromosomes 1-22 were obtained for our study.

Statistical Analysis
Single-SNP association analysis was first performed for each of the 344,893 SNPs using linear regression analysis with adjustments for sex, age, area, and body mass index. A multi-stage approach [11] was next conducted to identify multiple SNPs associated with PA. At the first stage, through single-SNP analysis, we selected the top 1000, top 2000, top 3000, and top 4000 SNPs which exhibit strong associations with PA for dimensional reduction. At the second stage, multiple-SNP analysis was performed with the EN regularization method, by utilizing a subset of SNPs chosen at the first stage. At the final stage, the bootstrap selection stability (BSS) measure was computed for each SNP, which indicated how consistently a SNP was replicated in bootstrap datasets. SNPs with high BSS values tend to have a higher chance of being replicated. More information about a multi-stage approach is described in Cho et al. [11].
To investigate the biological significance of PA-related genetic variants, we mapped the identified SNPs to an exon/intron and performed pathway enrichment analysis. All pathways related to the identified genes were investigated via the Kyoto Encyclopedia of Genes and Genomes (KEGG) database. Pathways were evaluated by the over-representation statistic and the Expression Analysis Systematic Explorer (EASE, Database for Annotation, Visualization and Integrated Discovery (DAVID), Frederick, MD, USA) score [39]. EASE calculates over-representation using Fisher's exact probability with respect to the total number of genes assayed and annotated within each system to measure the gene-enrichment in annotation terms [39]. p-values less than 0.1 defined as significantly regulated pathway. The single-marker association tests and adopted multi-stage approach were conducted using the PLINK [40] and R [41] software.

Conclusions
In this study, we demonstrated that the joint identification of SNPs could enable the identification of multiple SNPs with good predictive power for PA and a pathway enriched for PA. Previous genetic studies have focused on the relationship between PA and health (or fitness) [37,[42][43][44]. Future studies are recommended to determine pertinent genetic factors that influence health-or fitness-related PA. In view of the large population diversities in GWA studies, a systematic comparison is needed between our Korean results and those derived from other populations.