Effects of circadian clock genes and health-related behavior on metabolic syndrome in a Taiwanese population: Evidence from association and interaction analysis

Increased risk of developing metabolic syndrome (MetS) has been associated with the circadian clock genes. In this study, we assessed whether 29 circadian clock-related genes (including ADCYAP1, ARNTL, ARNTL2, BHLHE40, CLOCK, CRY1, CRY2, CSNK1D, CSNK1E, GSK3B, HCRTR2, KLF10, NFIL3, NPAS2, NR1D1, NR1D2, PER1, PER2, PER3, REV1, RORA, RORB, RORC, SENP3, SERPINE1, TIMELESS, TIPIN, VIP, and VIPR2) are associated with MetS and its individual components independently and/or through complex interactions in a Taiwanese population. We also analyzed the interactions between environmental factors and these genes in influencing MetS and its individual components. A total of 3,000 Taiwanese subjects from the Taiwan Biobank were assessed in this study. Metabolic traits such as waist circumference, triglyceride, high-density lipoprotein cholesterol, systolic and diastolic blood pressure, and fasting glucose were measured. Our data showed a nominal association of MetS with several single nucleotide polymorphisms (SNPs) in five key circadian clock genes including ARNTL, GSK3B, PER3, RORA, and RORB; but none of these SNPs persisted significantly after performing Bonferroni correction. Moreover, we identified the effect of GSK3B rs2199503 on high fasting glucose (P = 0.0002). Additionally, we found interactions among the ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, and RORB rs972902 SNPs influenced MetS (P < 0.001 ~ P = 0.002). Finally, we investigated the influence of interactions between ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, and RORB rs972902 with environmental factors such as alcohol consumption, smoking status, and physical activity on MetS and its individual components (P < 0.001 ~ P = 0.002). Our study indicates that circadian clock genes such as ARNTL, GSK3B, PER3, RORA, and RORB genes may contribute to the risk of MetS independently as well as through gene-gene and gene-environment interactions.

Introduction Circadian rhythms are naturally recurring cycles that regulate the timing of biological events such as the sleep-wake cycle and energy metabolism [1]. The intracellular molecular machinery underlying circadian rhythms implicates that circadian oscillations are controlled and maintained by a set of core circadian clock genes, including the aryl hydrocarbon receptor nuclear translocator like (ARNTL), clock circadian regulator (CLOCK), cryptochrome circadian clock 1 (CRY1), cryptochrome circadian clock 2 (CRY2), neuronal PAS domain protein 2 (NPAS2), period circadian clock 1 (PER1), period circadian clock 2 (PER2), and period circadian clock 3 (PER3) genes [2]. The circadian clock genes are prominently involved in the regulation of metabolism and energy homeostasis because circadian rhythms in turn impact on metabolic activity and metabolism [3]. Circadian disruptions have been well documented in adverse effects on human health, influencing lipid and glucose homeostasis, inflammation, and cardiovascular functions [4].
Accumulating evidence indicates that the circadian clock genes may contribute to the development of metabolic syndrome (MetS). Animal studies showed that a loss-of-function mutation of the circadian clock genes, such as Clock and Arntl, may result in metabolic abnormalities of lipid and glucose homeostasis, hallmarks of MetS [5,6]. In humans, evidence has been reported for a potential association of MetS with the CLOCK [7] and CRY2 [8] genes. Previously, Scott et al. showed that a three-marker haplotype (rs4864548, rs3736544, and rs1801260) of the CLOCK gene was associated with MetS, although there were no significant associations between any of these three SNPs and MetS individually [7]. Additionally, Kovanen et al. indicated that CRY2 rs75065406 had a nominally significant association with MetS, but this association did not remain significant after correcting for multiple testing [8]. Moreover, circadian clock genes have been shown to link with the individual components of MetS. Woon et al. found that a three-marker haplotype (rs6486121, rs3789327, and rs969485) of the ARNTL gene was associated with hypertension, an individual component of MetS [9]. Further, Englund et al. reported that NPAS2 rs11541353 was associated with hypertension and PER2 #10870 was linked to blood glucose levels [10]. Kovanen et al. also indicated a nominally significant association of elevated blood pressure with the CRY1 haplotype (rs4964513 and rs12821586); however, the significance did not persist after adjusting for multiple testing [8]. In addition, it has been suggested that circadian clock genes such as ARNTL, CRY1, and PER2 are expressed in human adipose tissue, and their gene expressions were related to the individual components of MetS such as waist circumference [11]. Furthermore, several animal studies indicated that metabolic complications such as MetS, dyslipidemia, glucose intolerance, hypoinsulinaemia, and diabetes can result from deletion of the Arntl, Clock, or Per3 genes, suggesting miscommunication between the circadian clock and metabolic pathways may lead to metabolic disorders [12][13][14].

Study population
This study incorporated Taiwanese subjects from the Taiwan Biobank, which gathered information and specimens from participants in recruitment centers across Taiwan in 2013-2015 [22][23][24][25][26][27]. The Taiwan Biobank was mainly funded by the Taiwan government and aimed to facilitate collaboration for public health-related research concerning local common chronic diseases [22]. The study cohort consisted of 3,000 participants. Inclusion criteria were individuals who were aged 30 years or over and were self-reported as being of Taiwanese Han Chinese ancestry [24]. Participants with a history of cancer were excluded [24]. Ethical approval for the study was granted by the Institutional Review Board of the Taiwan Biobank before conducting the study (approval number: 201506095RINC). Each subject signed the approved informed consent form. All experiments were performed in accordance with relevant guidelines and regulations.
Current alcohol drinkers were defined as persons who reported drinking more than 150 ml of alcohol per week during the last six months [23]. Current smokers were defined as persons who reported smoking every day or some days during the last six months [23]. Physical activity was defined by the amount of excise activity for more than three times and more than 30 minutes each time in each week.

Metabolic syndrome
Measurements of metabolic traits (including as waist circumference, triglyceride, high-density lipoprotein cholesterol, systolic and diastolic blood pressure, and fasting glucose) were obtained when participants underwent general health examinations [22][23][24]. To assess MetS, we used the International Diabetes Federation (IDF) definition, which requires that the participant represented by central obesity (defined as waist circumference ! 90 cm in male subjects and ! 80 cm in female subjects) plus the presence of two or more of the following four components: (1) triglycerides ! 150 mg/dl; (2) HDL cholesterol < 40 mg/dl in male subjects and < 50 mg/dl in female subjects; (3) systolic blood pressure ! 130 mmHg or diastolic blood pressure ! 85 mmHg; and (4) fasting plasma glucose ! 100 mg/dl [28]. Blood pressure was based on the average of two measurements.
Genotyping DNA was isolated from blood samples using a QIAamp DNA blood kit following the manufacturer's instructions (Qiagen, Valencia, CA, USA). The quality of the isolated genomic DNA was evaluated using agarose gel electrophoresis, and the quantity was determined by spectrophotometry [29]. SNP genotyping was carried out using the custom Taiwan BioBank chips and run on the Axiom Genome-Wide Array Plate System (Affymetrix, Santa Clara, CA, USA).
To efficiently obtain maximal genetic information from Taiwanese Han Chinese samples, the custom Taiwan BioBank chips were designed by using SNPs on the Axiom Genome-Wide CHB 1 Array (Affymetrix, Inc., Santa Clara, CA, USA) with minor allele frequencies (MAFs) ! 5%, by using SNPs in exons with MAFs > 10% on the Human Exome BeadChip (Illumina, Inc., San Diego, CA, USA), and by using SNPs previously reported in ancestry information panels, cancer studies, and pharmacogenetic studies. [24]. In this study, the SNP panel consisted of variants from the following 29 circadian clock genes:

Statistical analysis
Categorical data were evaluated using the chi-square test. We conducted the Student's t test to compare the difference in the means from two continuous variables. To estimate the association of the investigated SNP with MetS, we conducted a logistic regression analysis to evaluate the odds ratios (ORs) and their 95% confidence intervals (CIs), adjusting for covariates including age and gender [30]. Furthermore, we estimated the association of the investigated SNP with individual components of MetS by using logistic regression analysis, adjusting for age and gender [31]. The genotype frequencies were assessed for Hardy-Weinberg equilibrium using a χ 2 goodness-of-fit test with 1 degree of freedom (i.e. the number of genotypes minus the number of alleles). Multiple testing was adjusted by the Bonferroni correction. To address the issue of multiple testing, we also calculated Q values and false discovery rates (FDRs) using the qvalue function in an R package [32]. The criterion for significance was set at P < 0.05 for all tests. Data are presented as the mean ± standard deviation.
To investigate gene-gene and gene-environment interactions, we employed the generalized multifactor dimensionality reduction (GMDR) method [33]. We tested two-way up to fourway interactions using 10-fold cross-validation. The GMDR software provides some output parameters, including the testing accuracy and empirical P values, to assess each selected interaction. Moreover, we provided age and gender as covariates for gene-gene and gene-environment interaction models in our interaction analyses. Permutation testing obtains empirical P values of prediction accuracy as a benchmark based on 1,000 shuffles. In order to correct for multiple testing, we applied a conservative Bonferroni correction factor for the number of tests employed in the GMDR analysis.
Based on the effect sizes in this study, the power to detect significant associations was evaluated by QUANTO software (http://biostats.usc.edu/Quanto.html). Table 1 describes the demographic and clinical characteristics of the study population, including 533 MetS subjects and 2,467 non-MetS subjects. The MetS prevalence in our cohort was 17.8%. As shown in Table 1, the distribution was well matched for gender, but not for age. Moreover, there was a significant difference in waist circumference, triglyceride, HDL, systolic blood pressure, diastolic blood pressure, and fasting glucose between the MetS and non-MetS subjects (Table 1; all P < 0.0001, respectively). Furthermore, there was a significant difference in current alcohol consumption (P = 0.029) and smoking status (P = 0.001) between the MetS and non-MetS subjects. However, no significant difference was found between participants with and without the MetS in level of physical activity.
Moreover, the OR analysis showed risk genotypes of variants of ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, and RORB rs972902 after adjusting for covariates, indicating an increased MetS risk among the subjects ( Table 2). As demonstrated in Table 2 for the RORA rs8034880 SNP, there was an indication of an increased MetS risk among the MetS and non-MetS subjects after adjustment of covariates such as age and gender for genetic models, including the dominant model (OR = 1.44; 95% CI = 1.19-1.74; P = 0.0002). Similarly, there was an indication of an increased risk of MetS among the subjects after adjustment of covariates for genetic models in the ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, and RORB rs972902 SNPs (Table 2).
Next, Table 3 and S4 Table show the OR analysis of the ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, and RORB rs972902 SNPs with the individual components of MetS including (a) high waist circumference vs. normal waist circumference  Table 3 for the GSK3B rs2199503 (P = 0.0002), RORA rs8034880 (P = 0.02), and RORB rs972902 (P = 0.015) SNPs, there was an indication of an increased risk of high fasting glucose among the subjects after adjustment of covariates for genetic models. The effect of GSK3B rs2199503 on high fasting glucose remained significant after Bonferroni correction (P < 0.05/75 = 0.0007). Additionally, there was a significantly increased risk of high waist circumference (P = 0.0013) among the subjects for PER3 rs10746473; however, this association did not persist after Bonferroni correction (P < 0.05/ 75 = 0.0007).
In addition, the GMDR analysis was used to assess the impacts of combinations between the five SNPs in MetS and its individual components including age and gender as covariates. Table 4 and S5 Table summarize the results obtained from GMDR analysis for two-way genegene interaction models in influencing MetS with covariate adjustment. As shown in Table 4, there was a significant two-way model involving ARNTL rs10832020 and RORA rs8034880 (P < 0.001), GSK3B rs2199503 and RORA rs8034880 (P < 0.001), GSK3B rs2199503 and RORB rs972902 (P = 0.002), PER3 rs10746473 and RORA rs8034880 (P < 0.001), and PER3 rs10746473 and RORB rs972902 (P < 0.001). The effects of these two-way models remained significant after Bonferroni correction (P < 0.05/10 = 0.005), indicating a potential gene-gene interaction between RORA and ARNTL, between RORA and GSK3B, between RORA and PER3, between RORB and GSK3B, and between RORB and PER3 in influencing MetS. Similarly, there were two-way gene-gene interaction models in influencing individual components such as high waist circumference (P = 0.003) or high fasting glucose (P = 0.007) as shown in S6 Table. There were also a three-way model (P = 0.003) in influencing MetS and a four-way  Table). However, the effects for these models did not persist after Bonferroni correction (P < 0.05/18 = 0.0028). Furthermore, Table 5 and S7 Table show the GMDR analysis of gene-environment interaction models in MetS and its individual components using age and gender as covariates. As shown in Table 5 for MetS, there were a significant two-way model involving RORB rs972902 and smoking (P < 0.001). The effect of the two-way model remained significant after Bonferroni correction (P < 0.05/18 = 0.0028), indicating a potential gene-environment interaction Table 3 ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, and RORB rs972902). between RORB and smoking in influencing MetS. Additionally, there was a significant threeway model (P = 0.001) after Bonferroni correction (P < 0.0028), indicating a potential geneenvironment interaction among PER3, RORB, and smoking in influencing MetS. Similarly, there were significant two-way up to four-way gene-environment interaction models (all P < 0.001) in influencing individual components such as low HDL after Bonferroni correction (P < 0.0028). In addition, there were significant three-way and four-way models (P = 0.001 and 0.002, respectively) in influencing high fasting glucose after Bonferroni correction (P < 0.0028). Finally, statistical power analysis revealed that the present study had a 99.9% power to detect associations of ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, or RORB rs972902 with MetS among the MetS and non-MetS subjects.

Discussion
Our association study is the first to date to examine whether the main effects of 881 SNPs in 29 circadian clock-related genes are significantly associated with the risk of MetS and its individual components independently and/or through gene-gene interactions among Taiwanese individuals. We also investigated the relationship between these genes and health-related behaviors to examine whether these genes confer a risk of MetS according to its effect on geneenvironment interactions. A growing body of evidence indicates that genetic association studies using individual genes may overlook the associations that can only be identified if genegene and gene-environment interactions are explored [34]. Therefore, we further identified five key SNPs in five circadian clock genes such as ARNTL, GSK3B, PER3, RORA, and RORB to assess gene-gene and gene-environment interactions although the association with MetS did not remain significant after performing Bonferroni correction (P < 6 x 10 −5 ). Here, we report for the first time that gene-gene interactions of the ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, RORA rs8034880, and RORB rs972902 SNPs may contribute to the etiology of MetS. Moreover, there were gene-environment interactions between ARNTL rs10832020, GSK3B rs2199503, PER3 rs10746473, and RORB rs972902 with health-related behaviors, such as alcohol consumption, smoking status, and physical activity. Finally, our data revealed that the GSK3B rs2199503 SNP was associated with the individual component of MetS such as high fasting glucose.
Regarding the ARNTL gene, we found a positive association of MetS with 4 SNPs in the ARNTL gene, especially the rs10832020 SNP. The ARNTL gene is located on chromosome 11p15 and encodes a basic helix-loop-helix protein to build the ARNTL/CLOCK heterodimeric protein with CLOCK, another basic helix-loop-helix protein [9]. The ARNTL gene has been associated with metabolic regulation [9,13]. Woon et al. identified a three-marker haplotype (rs6486121, rs3789327, and rs969485) of the ARNTL gene as an indicator of hypertension [9]. Moreover, Shimba et al. demonstrated that deletion of the Arntl gene induced MetS and dyslipidemia in an animal model [13]. In the present study, we detected an association of high blood pressure with ARNTL rs3789327 (P = 0.0346), but not with ARNTL rs6486121 and rs969485 (data not shown). In addition, we did not detect an association of MetS with ARNTL rs6486121, rs3789327, and rs969485 in this study.
Further, five potential candidate SNPs for MetS were identified at the GSK3B gene, particularly the rs2199503 SNP. We also suggested an association of GSK3B rs2199503 with high fasting glucose. The GSK3B gene is located on chromosome 3q13.3 and encodes a subfamily of glycogen synthase kinase (GSK), which comprises two isoforms GSK3A and GSK3B. It has been implicated that GSK3 regulates numerous cellular functions such as metabolism, circadian rhythms, neurological degeneration, and immune responses [35]. Chen et al. reported that the expression of phosphatidylinositol-3 kinase (PI3K)-dependent GSK3 stimulates the production of adiponectin and protects against diet-induced MetS in an animal study [36].
On another note, there was an association of MetS with 6 SNPs within the PER3 gene, notably the rs10746473 SNP. Our analysis also indicated that there was an association of PER3 rs10746473 with high waist circumference. The PER3 gene, located on chromosome 1p36.23, encodes components of the circadian rhythms of locomotor activity, metabolism, and feeding behavior [37]. Dallmann and Weaver reported that the Per3 mutant mice led to higher body mass on a high fat diet, suggesting Per3 may play a role in body mass regulation [38]. In addition, Costa et al. demonstrated that Per3 knock-out mice displayed glucose intolerance as well as altered body composition (such as increased adipose and decreased muscle), indicating Per3 as potent mediator of cell fate [14]. In an animal study, Pendergast et al. also showed that the function of Per3 was tissue specific such that the loss of functional PER3 had an effect on circadian period in about half of the analyzed tissues, demonstrating Per3 as an important player in the mammalian circadian system [39].
Moreover, our results are the first to raise the possibility that 43 SNPs in the RORA gene and 6 SNPs in the RORB gene may contribute to the susceptibility for MetS, especially RORA (rs17237367, rs58469372, rs12591650, rs12594188, rs17270446, rs11630062, rs8029848, rs8034880, rs72752802) and RORB rs972902. Our results also suggested an association of high fasting glucose with the RORA rs8029848, RORA rs8034880, and RORB rs972902 SNPs; however, the significance did not persist after adjusting for Bonferroni correction. The proteins encoded by RORA (located on chromosome 15q22.2) and RORB (located on chromosome 9q2) constitute a subfamily of nuclear hormone receptors [40]. The RORA and RORB genes have been implicated in the regulation of a wide variety of physiological processes, including metabolism, circadian rhythm, and inflammatory responses [40,41]. Kang et al. demonstrated that Rora-deficient mice showed a markedly reduced susceptibility to age-and diet-induced MetS, implicating Rora plays a critical role in the regulation of several aspects of MetS [41].
Previously, circadian clock genes have been reported to associate with MetS, such as a three-marker haplotype (rs4864548, rs3736544, and rs1801260) of the CLOCK gene [7] and CRY2 rs75065406 [8]. In addition, circadian clock genes have been shown to link with the individual components of MetS, including an association of NPAS2 rs11541353 with hypertension [10], an association of PER2 #10870 with high fasting glucose [10], and an association of CRY1 haplotype (rs4964513 and rs12821586) with elevated blood pressure [8]. In this study, we did not weigh CLOCK (rs4864548, rs3736544, and rs1801260), CRY1 (rs4964513 and rs12821586), CRY2 rs75065406, NPAS2 rs11541353, and PER2 #10870 due to lack of these SNPs in the custom chip. It is worth mentioning that the minor allele frequencies of CRY2 rs75065406 and NPAS2 rs11541353 are all 0% in Asian populations (such as Japanese and Han Chinese) as shown in the public data from the 1000 Genomes Project (http://www.1000genomes.org). Moreover, it should be noted that CRY1 (rs4964513 and rs12821586) are annotated within the mitochondrial transcription termination factor 2 (MTERF2) gene, which is in the vicinity of the CRY1 gene, as shown in the dbSNP database (https://www.ncbi.nlm.nih.gov/). On the other hand, we identified other potential candidate SNPs for MetS, including 3 SNPs in the CRY1 gene and 6 SNPs in the NPAS2 gene.
Intriguingly, another finding was that we further inferred the epistatic effects between ARNTL, GSK3B, PER3, RORA, and RORB in influencing MetS by using the GMDR approach. To our knowledge, no other study has been conducted to evaluate gene-gene interactions between these genes. Besides the statistical significance, the potential biological mechanism under the interaction models was our concern. The functional relevance of the interactive effects of ARNTL, GSK3B, PER3, RORA, and RORB on MetS remains to be elucidated. In mammals, the molecular clock mechanism is currently viewed as a complex interplay of transcriptional feedback regulatory loops involving various circadian clock genes to generate and maintain circadian rhythms [40]. The core circadian clock gene ARNTL encodes the ARNTL protein, which is a basic helix-loop-helix protein that forms heterodimers with either CLOCK or NPAS2, two other basic helix-loop-helix proteins [2]. The ARNTL/CLOCK heterodimer initiates the transcription of various target genes such as PER1, PER2, PER3, CRY1, and CRY2 via E-box elements in the promoters of target genes [2]. In turn, the resulting PER and CRY proteins inhibit further ARNTL/CLOCK transcriptional activity, thereby allowing a new cycle to start from a level of low transcriptional activity. In order to improve the robustness of the previous loop, an additional feedback loop involves the nuclear receptors NR1D1 (or Reverb alpha)and RORA, which regulate the circadian transcription of the ARNTL gene and provide a good timing of the core clock machinery [40].
Remarkably, the GMDR analysis of gene-environment interactions reflected the interplay among PER3, RORB, and smoking in influencing MetS. Similarly, we tracked down several significant gene-environment interaction models in influencing individual components such as low HDL and high fasting glucose. It has been highlighted that common diseases are known to have a considerable genetic contribution, although only a proportion of complex diseases overall is explained by the established candidate genes, pointing out that environmental effects and gene-environment interactions will be our key concern in the outlook for genomic studies [42].
In terms of the average ages in the MetS and non-MetS groups, subjects with MetS were significantly older than those without MetS in this study (Table 1). Consistent with our results, Hwang et al. reported that the prevalence of MetS in Taiwan increased with age from 5.2% in individuals aged 20-29 to 36.5% in individuals aged 70-79, reaching peak levels in the 7th decade of life for both men and women [43]. Similarly, the prevalence of the MetS increased with age (ranging from about 7% in persons aged 20-29, about 14% in persons aged 30-39, about 20% in persons aged 40-49, to about 44% in persons aged 60-69) in National Health and Nutrition Examination Survey (1988)(1989)(1990)(1991)(1992)(1993)(1994) in the United States [44].
This study has both strengths and limitations. The main weakness was that our observations require much further research to validate whether the findings are replicated in diverse ethnic groups [45][46][47]. This study may also fail to assess other relevant circadian clock genes such as the basic helix-loop-helix family member e41 (BHLHE41) gene due to lack of those genes in the custom chip. In future work, prospective clinical trials with other ethnic populations are needed to examine a thorough evaluation of the association and interactions of the investigated genes with MetS and its individual components [48,49]. On the other hand, a major strength of our study was that we employed health-related behavior data, which served an exclusionary opportunity to facilitate the interplay between the investigated genes and environmental factors.

Conclusions
In conclusion, we conducted an extensive analysis of the association as well as gene-gene and gene-environment interactions of the circadian clock genes with MetS and its individual components in Taiwanese subjects. Our findings indicate that the ARNTL, GSK3B, PER3, RORA, and RORB genes may affect the prevalence of MetS independently and/or through complex gene-gene and gene-environment interactions. Furthermore, the circadian clock genes may be a determinant of MetS component factors. Independent replication studies with larger sample sizes will likely demonstrate further insights into the role of the circadian clock genes found in this study.