Investigation of gene-gene interactions of clock genes for chronotype in a healthy Korean population

Chronotype is an important moderator of psychiatric illnesses, which seems to be controlled in some part by genetic factors. Clock genes are the most relevant genes for chronotype. In addition to the roles of individual genes, gene-gene interactions of clock genes substantially contribute to chronotype. We investigated genetic associations and gene-gene interactions of the clock genes BHLHB2, CLOCK, CSNK1E, NR1D1, PER1, PER2, PER3, and TIMELESS for chronotype in 1,293 healthy Korean individuals. Regression analysis was conducted to find associations between single nucleotide polymorphism (SNP) and chronotype. For gene-gene interaction analyses, the quantitative multifactor dimensionality reduction (QMDR) method, a nonparametric model-free method for quantitative phenotypes, were performed. No individual SNP or haplotype showed a significant association with chronotype by both regression analysis and single-locus model of QMDR. QMDR analysis identified NR1D1 rs2314339 and TIMELESS rs4630333 as the best SNP pairs among two-locus interaction models associated with chronotype (cross-validation consistency [CVC] = 8/10, p = 0.041). For the three-locus interaction model, the SNP combination of NR1D1 rs2314339, TIMELESS rs4630333, and PER3 rs228669 showed the best results (CVC = 4/10, p < 0.001). However, because the mean differences between genotype combinations were minor, the clinical roles of clock gene interactions are unlikely to be critical.


Introduction
Many biological processes have circadian rhythms with a cycle of approximately 24 hours. Such circadian rhythms appear at the cellular, molecular, tissue, and organ levels in humans. Both the suprachiasmatic nucleus, which is the internal master clock, and various environmental stimuli (e.g., the light-dark cycle) may play roles together in regulation of circadian rhythm. The sleep-wake cycle is one of the most distinct circadian rhythms. Humans have a preferred timing of sleeping and waking, the so-called chronotype. The degree of morning preference shows a continuum from morning type (with an advanced sleep phase) to evening type (with a delayed sleep phase) [1]. The chronotype is known to be determined in a complex manner by age, sex, and various environmental factors, including level of light exposure [2][3][4]. Developmental changes in chronotype occur; these include earlier preference during childhood, later preference in adolescence and early adulthood, and gradually earlier preference with advancing age [2].
Chronotype has been studied extensively in a number of ways. For convenience in many studies, chronotype has been measured subjectively and arbitrarily classified as morning, intermediate, and evening types, according to the Composite Scale of Morningness (CS) [5,6]. Chronotype shows a normal or near-normal distribution in the population; the heritability of chronotype has been shown to range between 21% and 52% [7]. Therefore, chronotype is presumably a polygenic quantitative trait [8,9]. Circadian rhythms are regulated by a set of circadian genes in mammals [9,10]. These circadian genes are also known as clock genes; they include ADCYAP1, ARNTL, BHLHB2, BHLHB3, CLOCK, PER1, PER2, PER3, NR1D1, NPAS2, NR1D1, THRA, CSNK1D, CSNK1E, CRY2, RORA, RORB, TIMELESS, and VIP. Notably, the list continues to grow due to discoveries in this field. Recently, three genome-wide association studies (GWASs) have been performed in people of European ancestry [11][12][13]. The results indicated meaningful overlap of genes that exhibited significant associations with chronotype. All three GWASs supported the association of four genes-PER2, RGS16, FBXL13, and AK5-with chronotype [14]. A meta-analysis of these three GWASs identified 351 loci associated with chronotype, of which 327 loci were novel and 24 loci had been reported in other GWASs [15].
GWASs can identify only common genetic variants, with small or modest effects for complex traits. Chronotype is a complex trait, which may be influenced by polygenes. A recent GWAS found only common genetic variants with small individual genetic effects on chronotypes. Only a subset of the significant variants found in that GWAS were known clock genes or have been suggested to function as clock genes. There were many other genetic variants without known functions as clock genes. Genetic control of chronotype must occur in multiple ways, with individual genotypic associations of polygenes, haplotypic associations of more crucial genes, and gene-gene interactions of multiple biologically related functional genes. Here, we investigated possible roles for genegene interactions of clock genes in chronotype in a Korean population.

Participants
The study population consisted of 1,293 unrelated healthy Korean individuals. These were the same individuals included as controls in our previous study of bipolar disorder [16]. They consisted mostly of college students, nurses, and public officials, who were recruited after a brief psychiatric interview. Potential participants were excluded if they reported a history of a psychotic disorder, mood disorder, anxiety disorder, substance use disorder, brain trauma, or intellectual disability. All participants were informed of the purpose and methods of the study and provided informed consent before enrollment. The Ethics Committee of Eulji General Hospital approved the study protocol (IRB No. 2016-08-009).

Measurement of chronotype
Chronotype was measured using a self-reported questionnaire. The CS is a 13-item questionnaire, which assesses individual differences in the time of day a person prefers to carry out various activities; it classifies people as morning, intermediate, or evening types [5]. Three items are scored on a five-point scale from 1 to 5; the other 10 items are scored on a four-point scale, from 1 to 4. Higher scores indicate morning preference. All participants completed the CS questionnaire.

Genotyping
The clock genes investigated in this study were BHLHB2, CLOCK, CSNK1E, NR1D1, PER1, PER2, PER3, and TIMELESS. These eight genes were analyzed for 19 different tag single nucleotide polymorphisms (SNPs) with minor allele frequencies exceeding 5% in Asian populations. DNA was extracted from blood and SNPs were genotyped using the TaqMan method (Applied Biosystems, Foster City, CA, USA). Table 1 presents a summary of the minor allele frequencies and chromosomal locations of the SNPs.

Statistical analysis
Individual SNPs were examined for Hardy-Weinberg equilibrium; two SNPs violating Hardy-Weinberg equilibrium were removed. Each SNP association with CS score was analyzed by simple regression analysis. Haplotype association with CS was also analyzed by PLINK if more than two SNPs for each gene were included [17].
Gene-gene interactions were analyzed using the quantitative multifactor-dimensionality reduction (QMDR) method, an extension of the multifactor-dimensionality reduction (MDR) algorithm to work with quantitative or continuous phenotypes [18]. The MDR method is one a commonly used method for detection and characterization of high-order gene-gene or gene-environment interactions in case-control studies; this comprises a nonparametric combinatorial approach that reduces the number of dimensions [19]. For each multi-locus genotype combination, QMDR calculates the mean value of phenotype and compares it to the overall mean to determine the genotype combination is high risk or low risk. By pooling all the genotypes into either high-risk or low-risk groups, a new binary attribute is created. The t-test is used to com-pare the means between high and low risk groups using a t-test and t-statistic is used as a training score to choose the best model. In QMDR, the training and testing score are defined by t-test statistic. The training score is used to determine the best K-order interaction model. QMDR use 10-fold cross-validation and cross-validation consistencies (CVCs) of each model chosen are recorded. The best overall QMDR model is selected as that with the maximum testing score and highest cross-validation consistency. To estimate the p-values of the chosen model, empirical null distribution is used [18].
In this study, interactions of up to three loci were tested using 10-fold cross-validation in a search considering all possible SNP combinations. SNP combination with maximum CVC was considered the best model. p-values were determined empirically by 1,000-fold permutations of case and control labels.

Results
The study population consisted of 481 male participants (37.2%) and 812 female participants (62.8%). Mean ages were 27.5 ± 8.3 years for male participants and 23.7 ± 3.5 years for female participants. Mean CS scores were 32.3 ± 6.4 for male participants and 29.7 ± 5.7 for female participants. The classification of chronotype using total CS score revealed 305 evening type participants (23.6%), 911 intermediate type participants (70.5%), and 77 morning type participants (6.0%) in the study population. The distributions of total CS score according to age and chronotype are shown in Table 2.
Two SNPs of PER1 were excluded from further analysis because they were found to deviate from Hardy-Weinberg equilibrium. Ultimately, 17 SNPs of seven genes were analyzed. In regression analyses, no individual SNP showed a significant association with CS score (Table 3). There were no significant haplotype associations with CS score for any of the genes with more than two SNPs in this study (Table 4). On QMDR analyses, no single locus was found to be associated with chronotype, similar to the results of regression analysis. NR1D1 rs2314339 and TIMELESS rs4630333 were significantly associated with chronotype in a two-locus model (CVC = 8/10, p = 0.041). In the three-locus models, NR1D1 rs2314339, TIMELESS rs4630333, and PER3 rs228669 showed the strongest association with chronotype (CVC = 4/10, p < 0.001). A summary of QMDR results with CVC > 1/10 is presented in Table 5.

Discussion
We hypothesized that circadian genes play an important role in chronotype regulation and that there are gene-gene interaction effects on chronotype. We identified the best interaction models for two and three loci, as well as statistical significances of the best interaction models for chronotype in a Korean population, using the QMDR method and corresponding permutation test. However, it    was difficult to conclude that there were clinically significant genegene interaction effects on chronotype based on our results, since the mean differences in total CS score between genotype combinations were minor. Human chronotype is a heritable polygenic trait, and many groups have searched for the genes involved in chronotype. Many clock genes have been considered as strong candidate genes for chronotype because of their biological functions in circadian networks. There have been many single-gene association studies to identify genetic factors for chronotype. The first report of an association between chronotype and clock genes was related to the 3′-untranslated region of the CLOCK gene (rs1801260) [20]. This finding was not replicated in other populations [21][22][23], but was replicated in a Japanese population [24]. Multiple other clock genes have also been studied. The PER3 gene was shown to be associated with delayed sleep phase syndrome. The studies thus far have focused on a variable number of tandem repeats region located in exon 18 of PER3, but the results of association analysis differed among studies. Specifically, a shorter allele was associated with delayed sleep phase syndrome [25]; a longer allele was also reportedly associated with delayed sleep phase syndrome [26] and was reported to have predictive value in response to sleep loss [27][28][29]. PER1 and PER2 were reported to be associated with advanced sleep phase syndrome in a British population [27,28], but not in a Japanese family [30]. Overall, the results were inconsistent and sometimes contradictory. Table 6 presents a summary of genetic association findings between clock genes and chronotype.
There have been several GWASs regarding chronotype in European populations using data from 23andMe and the UK Biobank [11][12][13]15]. There was a great deal of consistency across these studies, with nine genes identified in two of the three GWASs. Among these nine genes, four showed significant hits and were found in all three GWASs. Thus far, GWASs have been conducted only in populations of European ancestry. Emmanuel and von Schantz reported that six of 41 alleles identified earlier by GWASs in European participants (in the genes RGS16, PER2, and AK5, as well as between the genes APH1A and CA14) were absent from some non-European populations. This study underscores the ancestral diversity of circadian genetics and suggests that future studies should be performed in populations of East Asian ancestry [31].
The proteins encoded by various clock genes cooperate physically with each other and act as transcription factors. Combinations of polymorphisms in these genes may affect phenotype. Therefore, combined analysis of the effects of different clock genes may be more accurate and more revealing, compared to analyses of single genes. A few previous studies examined the effects of gene-gene interactions on chronotype. One study in Korean college students reported a significant interaction for CS score between CLOCK gene 3111 C/T and GNB3 825 C/T, according to regression analysis [32]. Later, the same group reported a genetic interaction for eveningness among ARNTL, PER2, and GNB3, according to MDR analysis [33]. Another study by Pedrazzoli et al. reported that a specific combination of polymorphisms in four clock genes was associated with diurnal preferences in a Brazilian population [34]. They chose four polymorphisms in four clock genes: PER2, PER3, CLOCK, and BMAL1. To the best of our knowledge, other than these studies, there have been no further analyses of the epistatic effects among clock genes for chronotype. Therefore, we attempted to identify gene-gene interactions among clock genes for chronotype in this study. We used QMDR to investigate gene-gene interactions to improve statistical power. We avoided subgrouping based on total CS score because grouping into morning, intermediate, or evening types could be arbitrary; moreover, we could address the quantitative score directly in our analysis. Multifactor dimensionality reduction is a common approach for identification of gene-gene interactions in case-control studies [19]. QMDR is an extension of MDR to handle quantitative phenotypes. Instead of comparing the case-control ratio of each multi-locus genotype to a fixed threshold as in MDR, QMDR compares the mean value of each multi-locus genotype to the overall mean [18].  In particular, gene-gene interactions among NR1D1 rs2314339, TIMELESS rs4630333, and PER3 rs228669 were significantly associated with chronotype in QMDR analyses in the present study. These SNPs did not show any associations as individual SNPs or haplotypes. NR1D1 (nuclear receptor subfamily 1, group D, member 1) is the gene encoding REV-ERBα, located on chromosome 17q21.3. REV-ERBα suppresses the transcription of BMAL1 mRNA [35,36], while BMAL1 activates REV-ERBα; this comprises a feedback loop of the mammalian circadian oscillator. NR1D1 has been reported to show an association with chronotype in healthy Korean young adults [37]. Kang et al. [37] reported a significant association with rs12941497 of NR1D1, but no association with rs2314339, which showed significant gene-gene interactions with other SNPs in the present study. The TIMELESS gene was reportedly associated with morningness-eveningness in healthy university students [38]. The PER3 gene showed conflicting results. Several positive associations for chronotype in European populations have been reported. In addition, negative findings were also reported for European populations [39][40][41] and Asian populations [42,43]. The circadian clock system consists mainly of transcription-translation feedback loops in the internal timekeeping clock. Heterodimers of BMAL (brain and muscle Arnt-like protein-1) and CLOCK (circadian locomotor output cycles kaput) activate the transcription of PER (period) and CRY (cryptochrome) genes. CRY and PER suppress transcriptional activity of BMAL1/ CLOCK [54,55]. Because NR1D1 is related to BMAL1 and PER3 was reported to show a strong genetic interaction with BMAL1, our main finding supports this hypothesis. Unfortunately, no biological mechanism has been reported for chronotype that includes positive gene-gene interactions of all three genes identified in this study. Further studies are needed to understand the complicated biological molecular network of circadian clocks in humans.
Consistent and clear phenotyping is critical when performing genetic studies. We used CS score as a proxy for actual human diurnal preference. However, it evaluates both diurnal preference and sleep homeostasis. There are many ways to determine the timing of the circadian system (e.g., core body temperature monitoring, dim light melatonin onset, and actigraphy-or diary-based midpoint of sleep); self-reported diurnal preference is only a proxy method. Most GWASs have assessed chronotype using only a single question, such as "Are you naturally a night person or a morning person?" Non-precise phenotyping can produce unreliable significant findings that are unlikely to be replicated in subsequent studies [56]. Therefore, future phenotyping should include standardized self-reporting, clinical interviews, or objective assessment of sleep-wake periodicity, such as actigraphy. It is important to validate commonly used self-reported items. Comparisons between self-reported items and biological markers of circadian rhythms are needed to determine which questions are most closely associated with endogenous processes [14].
This study had some limitations that must be considered when interpreting our results. First, this study included only ethnically Korean individuals. Therefore, caution is needed when generalizing our results to other populations. As suggested in previous studies, chronotype and genes for chronotype are likely to differ according to ethnicity and/or between populations. Therefore, this study in a Korean population was necessary. Second, our participants were relatively young, with a mean age of 25 years. Because chronotype is affected by age, our results cannot be applied to other age groups. Third, because of resource limitations, we could include only 17 SNPs of seven circadian genes. Therefore, our results represent only a subset of the real-world epistatic interactions among clock genes. There are likely to be more complicated genegene interactions among circadian genes, as well as epistatic interactions between circadian genes and genes of other biological sys-tems, which are directly and indirectly related to the circadian system.
We could not conclude that clock genes play a critical role in determining chronotype, although QMDR suggested significant gene-gene interactions. Further studies are needed to investigate gene-gene interactions of additional clock genes, especially with respect to SNPs that repeatedly show significant associations in multiple GWASs and studies in populations other than those of European ancestry.

Conflicts of Interest
No potential conflict of interest relevant to this article was reported.