PDF
Citation
Print
INTRODUCTION
Type 2 diabetes mellitus (T2DM) is diagnosed based on discrete values of fasting plasma glucose (FPG) and glycosylated hemoglobin (HbA1c). Metabolic dysregulation presents as a continuum of abnormal glucose level that is above normal glucose level and include prediabetes and diabetes. There has been a number of studies that investigated the natural course in the development of T2DM [1]. Based on a prospective cohort specifically designed to study the trajectories of diabetes, those who progressed to diabetes had elevated FPG level 10 years before the development of diabetes compared to those who remained normal glucose tolerance [2]. There seems to be a rapid increase in FPG 2 to 6 years before diagnosis of diabetes [2,3]. In addition, FPG level increased as age increases, at a rate of 12.6 to 19.8 mg/dL per age decade [4]. However, there might be individual variation in the rate of longitudinal glucose change. This is also suggested by the fact that only about 25% of prediabetes subjects progress to T2DM in 5 years of follow-up [5]. The genetic contribution on this longitudinal change in FPG level has not been thoroughly investigated.
Genome-wide association studies (GWAS) on T2DM have identified more than 400 distinct genetic loci associated with T2DM [6,7] and nearly 120 loci for FPG and fasting insulin level [8]. However, it should be noted that genetic loci for FPG and T2DM does not show complete overlap [6]. This discrepancy could be due to multiple factors including different range of glycemia studied for FPG and T2DM, interaction by environmental factors, and most of the studies being cross sectional in nature. Although there have been several studies that investigated the genetic association for longitudinal change in FPG [7], there has not been a locus that reached genome-wide significance.
Therefore, the main interest lies in dependence of longitudinal change in FPG on genetic variants, and their interaction with time. To this end, there are two popular methods to deal with repeated measurements: linear mixed model (LMM) and generalized estimating equations (GEE). LMM allows for correlation within subject to vary by a specific pattern which produces models that better fit the data. Additionally, subjects with missing outcome and different numbers of visits can be included as long as the time intervals are correctly specified in the model. Because LMM uses maximum likelihood estimation, it is robust against data which is missing at random [9]. GEE has the same advantages as LMM as it allows correlation structures through working correlation matrix and it can take several types of covariates. However, the main difference is that GEE only provides population average estimates and not individual level information for random effects. Since it requires complete data or missing completely at random, it is less robust than LMM in missing data scenarios [9,10].
In this study, we hypothesized that using LMM will allow us to investigate the impact of genetic variants on FPG and its longitudinal change in an East Asian population. For this end, we used two large-scale cohorts, Korea Genome and Epidemiology Study (KoGES) and Gene-Environment Interaction and phenotype (GENIE) cohort, to perform GWAS using LMM and meta-analyzed the results.
METHODS
Study participants
The KoGES [11] was initiated in 2002, and 6,122 individuals who did not have T2DM at baseline examination were investigated [2]. The participants consisted of 2,847 males and 3,275 females, with age at baseline range of 40 to 69 years. FPG and post-challenge 2-hour plasma glucose were measured every 2 years from 2001 to 2012. The diagnosis of T2DM was defined according to the American Diabetes Association (ADA) criteria: FPG ≥126mg/dL; post-challenge 2-hour plasma glucose ≥200mg/dL after a 75-g oral glucose load; and HbA1c ≥6.5%. The subjects were followed up for 11 years on average, and the minimum and maximum follow-up time were 6 and 12 years, respectively. During the follow-up period, 790 participants developed incident T2DM, and their FPG levels after the first diagnosis of T2DM were not included in the analyses.
The GENIE cohort consisted of 7,999 participants who had visited Seoul National University Hospital Healthcare System Gangnam Center between 2003 and 2015 for health check-up. A total of 4,406 participants (2,604 males and 1,802 females) who did not have T2DM at their baseline visit were considered for inclusion in the present analyses. The number of visits to the Gangnam Center varied among the participants, and those who visited the center more than three times, but less than 13 times were ultimately including in our analysis. The participants were followed up for 6 years on average. The number of visits to the Gangnam Center varied among the participants, and those who visited the center more than three times, but less than 13 times were ultimately including in our analysis. The participants were followed up for 6 years on average. Most of them had annual or biennial visit for health check-up. Although follow-up intervals were not matched to those of KoGES, we used the same year-unit for follow-up duration in both cohorts. T2DM was diagnosed according to the ADA definition as described above. During the follow-up period, 237 participants developed incident T2DM and their FPG values after the first incidence of T2DM were excluded.
Genotyping and imputation
Genome-wide genotyping on the KoGES and GENIE cohorts was performed using the Affymetrix genome-wide human single nucleotide polymorphism (SNP) Array 5.0 and Affymetrix KOR_v1.0, respectively [12]. Variant calling was conducted with K-medoid algorithm approach [13]. We excluded any SNP marker with more than 5% missingness, minor allele frequency <0.05, deviation from Hardy-Weinberg equilibrium P<1.00×10−6, and any subjects with more than 5% missing genotype calls and sex inconsistency. Genotype data was managed with PLINK (whole genome data analysis toolset) and ONETOOL (a tool for family-based big data analysis) [14]. Imputation was conducted with IMPUTE2 (a computer program for phasing observed genotypes and imputing missing genotypes) using the cosmopolitan reference panel from 1000 Genomes Project phase 3. We used weighted genetic score using imputed genotype probabilities by IMPUTE2. Any imputed SNPs with imputation quality scores <0.4 were excluded. After quality control, 3,758,649 and 3,692,736 SNPs were used for association analyses of longitudinal change in FPG in the KoGES cohorts and GENIE cohorts, respectively. Supplementary Table 1 shows the number of variants in each annotation group categorized by ANNOtate VARiation (ANNOVAR) [15]. Intergenic SNPs and intronic SNPs accounted for 55% and 36% of the total, respectively, with 6% contributed by exonic non-coding RNA. Genotype dosage scores obtained from IMPUTE2 were used for association testing [16].
Genome-wide association analyses
Association analyses were performed with LMM using the lme4 package in R v.3.4.1 [17] with two random effects, intercept and slope over time, for each subject and correlations between two random effects were allowed. We considered several models and the best model was selected by Akaike’s information criterion. We tested a number of different models based on the following criterions: (1) whether to allow correlations between two random effects, intercept and slope over time, for each subject; (2) whether to consider the inverse variance of each time point dataset as weights; (3) whether to adjust body mass index (BMI) as covariate. The selected model has two random effects for intercept and slope over time for each subject, with nonzero correlations between them. SNP, sex, baseline age, BMI, elapsed time on FPG, the ten principal component scores, and the SNP-by-time interaction were included as covariates. For the j-th measurement of the i-th subject, our LMM for both the KoGES and GENIE cohorts can be expressed by
Once the LMMs were applied to both the KoGES and GENIE cohorts, their P values were combined with inverse variance-weighted meta-analyses using METAL software (http://csg.sph.umich.edu/abecasis/metal) [18]. The threshold for statistical significance in this model was P<5.00×10−8, which is conventionally considered to reflect genome-wide significance. Selection of an appropriate suggestive significance threshold was P<1.00×10−5. To verify no inflation in the statistical significance, genomic inflation factor (GIF) was calculated and the GIF close to 1 means no genomic inflation in the GWAS [19]. The SNP based heritability value reflects the relative proportion of phenotypic variance explained by all observed common SNPs. SNP based heritabilities of the average and annual increment of FPG were estimated with genome-wide complex trait analysis (GCTA) v1.91.7 [20]. To estimate SNP heritability of FPG, we first fit the linear model (LM) for repeatedly observed FPG of each subject as follows:
where indicates the average of the ith subject’s age at different time points, and βi0 and βi1 indicate the expected FPG when Ageij= and the average annual increment, respectively. and were estimated for all subjects in both cohorts and were analyzed with GCTA. and were used as the outcome variable, and and sex were used as covariates.
Expression quantitative trait loci analysis and differential expression in T2DM
We explored the expression quantitative trait loci (eQTL) for our GWAS result in the tissues associated with FPG or T2DM using publicly available Genotype-Tissue Expression (GTEx) datasets [21]. Pancreas and adipose tissues, which are known to have a significant influence on FPG or T2DM, were mainly targeted. The candidate gene’s expression was compared for T2DM based on Gene Expression Omnibus (GEO) database of GSE25724 [22]. The platform for GSE25724 is GPL96, [HG-U133A] Affymetrix Human Genome U133A Array, which includes seven normal pancreatic islets subjects and six T2DM pancreatic islets subjects. Gene differential analysis was performed using the ‘limma’ R package (R Foundation for Statistical Computing, Vienna, Austria). Furthermore, Combined Annotation Dependent Depletion (CADD) score was calculated to predict the scoring the deleteriousness of single nucleotide variants as well as insertion/deletions variants in the human genome. CADD score of 10 or greater means top 10% probable functional variants [23].
Investigation for known T2DM and glycemic variants
To investigate whether the already well-known variants of T2DM and FPG are also related to longitudinal change of FPG, we utilized publicly available data via the The DIAbetes Genetics Replication And Meta-analysis (DIAGRAM) consortium (http://diagram-consortium.org/), and GWAS summary results [6,8]. We identified variants associated with T2DM and FPG in genome-wide significance using conditional analyses as implemented in GCTA using GWAS results from 898,130 European-descent individuals [6] and trans-ethnic meta-analyses [6,8]. In total, 231 T2DM-related variants and 38 FPG-related variants were targeted and we considered Bonferroni-controlled P=2.20×10−4 and P=1.30×10−3 as significance threshold, respectively. In total, 231 T2DM-related variants and 38 FPG-related variants were targeted and we considered Bonferroni-controlling multiple-testing procedures [24].
Statistical power calculation
Statistical power to detect effect size 0.1 to 0.5 with sample size of 10,528 at the genome-wide significant level (P=5.00×10−8) is calculated using QUANTO software version 1.2.4 (https://quanto.software.informer.com) [25].
Ethical statements
All datasets involving human subject or genotype data must be approved by the Institutional Review Board (IRB) of the Seoul National University Hospital (H-1704-073-845). The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008. Written informed consent was obtained from all patients.
RESULTS
Clinical characteristics of study population
A total of 6,122 subjects in the KoGES cohort were investigated at regular intervals of every 2 years, whereas the 4,406 individuals in the GENIE cohort attended voluntary annual health screening. The baseline characteristics of each cohort are described in Table 1. The proportion of male subjects was 46% and 59% in the KoGES and GENIE cohort, respectively. During the follow-up period, approximately 10% and 5% of participants developed T2DM in the KoGES and GENIE cohort, respectively. The mean age at baseline was 51.5 years in KoGES, which was 5.8 years older than that of the GENIE cohort (45.7 years). The mean BMI was 24.4 kg/m2 in the KoGES and 23.1 kg/m2 in the GENIE. The mean FPG of the KoGES cohort at baseline was 84.1 mg/dL, which is about 9.8 mg/dL lower than that of the GENIE cohort (93.9 mg/dL). The trend of FPG change over time for the two cohorts is plotted in the Supplementary Fig. 1. The average increase of FPG per year was 0.59 mg/dL and 0.20 mg/dL in KoGES and GENIE cohorts, respectively. Information on genotyping platform and number of variants analyzed are presented in Table 1.
Genetic variants associated with repeated measures of FPG
We conducted GWAS meta-analysis using KoGES and GENIE cohorts. First, we investigated genetic variants associated with repeated measures of FPG using LMM as it allows to test both main effect of SNP and SNP×time interaction effect with longitudinal data. A total of four genetic variants were associated with repeated measure of FPG levels with genome-wide significance (P<5.00×10−8) using the LMM. These include rs78529- 720 (P=2.90×10−16) near glucose-6-phosphatase catalytic subunit 2 (G6PC2) and ATP binding cassette subfamily B member 11 (ABCB11), rs895636 (P=2.00×10−8) near six homeobox 3 (SIX3), rs2971670 (P=8.34×10−20) in glucokinase (GCK), and rs12222793 (P=4.30×10−10) near melatonin receptor 1B (MTNR1B) (Table 2). All of these variants have been reported previously to be strongly associated with FPG level [26-30]. Furthermore, the reported results show a trend in the same direction as observed in our study. Compared to single time point approach like LM analysis, multiple observations per person for LMM increase the power to detect a statistically significant genetic loci associated with FPG (Supplementary Table 2) [31]. This is evidenced by the fact that all four variants had lower P in LMM compared to LM. Quantile-quantile (Q-Q) plot and Manhattan plot for the genetic association using LMM of repeated measures of FPG revealed no evidence for inflation of the test statistics (GIF=1.016) (Fig. 1A and B). The regional association plots theses genetic loci are shown in Supplementary Fig. 2.
Genetic variants associated with longitudinal change in FPG
We further investigated genetic variants associated with longitudinal change in FPG using LMM SNP×time interaction term. The statistical model that was used had good control of type 1 error as evidenced by Q-Q plot (GIF=0.999) and Manhattan plot as shown in (Fig. 1C and D). We found one novel SNP, rs11187850, near phospholipase C epsilon 1 (PLCE1) to be associated with longitudinal change in FPG with genomewide significance (P=4.85×10−8) with CADD score of 10.18. In addition, three variants had suggestive evidence for an association (P<1.00×10−5): rs10947494 (P=3.64×10−6) near nudix hydrolase 3 (NUDT3), rs2414772 (P=6.30×10−6) near microRNA 6085 (MIR6085), and rs16959641 (P=2.64×10−6) near U6 snRNA biogenesis phosphodiesterase 1 (USB1) (Table 3). Regional association plots are shown using linkage disequilibrium information from Asian 1,000 Genomes in Fig. 2.
We fitted a simple LM with FPG as the outcome variable and follow-up time as the independent variable, categorized by the three genotypes (AA/AG/GG) of the rs11187850 variant which showed genome-wide significant association. The alternative allele, G, was significantly associated with increased slope of FPG over time in both cohorts as shown in Supplementary Fig. 3. Individuals with AG and GG genotype showed larger FPG increase of 0.183 mg/dL and 0.366 mg/dL per year compared to those with AA genotype, respectively. We further investigated whether SNP rs11187850 was associated with risk of T2DM and increased FPG. Although this variant was nominally associated with increased risk of T2DM (P=2.91×10−6) and increased FPG level (P=4.75×10−4) in LMM it did not pass the genome-wide significance threshold. Interestingly, this variant was an eQTL for NOC3 like DNA replication regulator (NOC3L) in pancreas (P=2.20×10−11), adipose-subcutaneous (P=6.40×10−11), and adipose-visceral (P=2.70×10−6) in GTEx. Decreased expression of NOC3L gene is also associated with islet dysfunction (estimate of the log2-fold-change corresponding to the non-diabetic is −1.217 and P=8.25×10−4) in GEO datasets (GSE25724).
Association of known glycemic variants with longitudinal change in FPG
A limited number of known genetic variants of T2DM, and FPG were nominally associated (P<0.05) with longitudinal change in FPG (Table 4). Among the already known 231 T2DM-related variants, five SNPs (rs2296172 [P=3.42×10−3] in MACF1, rs243021 [P=3.37×10−2] near LOC101927285 and MIR4432H, rs864745 [P=4.90×10−2] near JAZF1, rs10965250 [P=2.93×10−3] near CDKN2B-AS1 and DMRTA1 and rs9552911 [P=4.89×10−2] near SGCG) showed association with the SNP×time interaction effect. Among the already known 38 FPG-related variants, three SNPs (rs6943153 [P=3.44×10−2] near GRB10, rs10811661 [P=5.10×10−3] near CDKN2B-AS1 and DMRTA1 and rs2293941 [P=2.90×10−2] near PDX1-AS1) had an association with longitudinal change of FPG. However, none of the variants were significant when Bonferroni correction was applied. These findings suggest that the variants associated with T2DM or FPG do not have a significant impact on longitudinal change in glucose level over time.
Statistical power and heritability
Our study has statistical power of 80% with a type 1 error rate of 5.00×10−8 to detect genetic variants, having at least effect size of 0.28 (i.e., 0.28 standard deviation unit difference in the change of FPG over time) in the longitudinal change of FPG (Supplementary Fig. 4). The reason for the difficulty in identifying the SNPs affecting the change in FPG over time can be illuminated based on the GCTA analysis results. As shown in Supplementary Table 3, the GCTA estimated a 10% (P=3.52×10−3) marginal influence of the SNP effect, whereas the influence of the SNP-by-time interaction was only 4% (P=7.58×10−2). Based on these results, genetic effects on longitudinal change of FPG are likely small.
DISCUSSION
In this study, the genetic association of longitudinal change in FPG was analyzed in a total of 10,528 Koreans free from T2DM using KoGES, a large-scale community-based prospective cohort, and GENIE, a hospital-based longitudinal routine health checkup cohort. Since the FPG was repeatedly measured for a long period of time and the number of total measurements of each participants was inconsistent, the analysis was conducted using the LMM method, which is known to have better performance than GEE in this situation [9,10]. Our study population was relatively homogeneous and FPG was measured regularly in 2 years interval for KoGES, both of which lead to increased statistical power. We found a novel association between rs11187850 and longitudinal changes in FPG which passed the genome-wide significance threshold. This variant was associated with NOC3L expression level in pancreas and adipose tissue. Known genetic variants of FPG and T2DM were not associated with longitudinal change in FPG. To the best of our knowledge, our study is the first in non-Europeans to provide a comprehensive assessment of genetic risk factors of longitudinal change in FPG using GWAS meta-analysis.
Using LMM for repeated measure of FPG, we identified four known genetic variants associated FPG and one variant significantly associated with SNP×time interaction. The rs11187850 G allele was associated with more rapid deterioration of FPG over time. Interestingly, this variant was associated with metabolic traits of diastolic blood pressure (P=1.37×10−19), systolic blood pressure (P=2.80×10−15), hypertension (P=8.20×10−12), and ascending aorta diameter (P=1.50×10−9) in Europeans [32-35]. Even though it was located in an intron of PLCE1, this variant was actually associated with expression level of NOC3L. It has been shown that the NOC3L gene is involved in adipocyte differentiation and glucose metabolism [36,37]. Taken together, it could be possible that rs11187850 affects longitudinal change in FPG by regulating NOC3L gene expression in metabolically active tissue such as pancreas and adipose tissue. However, replication of our finding and further functional investigations are required. There was no other genetic variant that was associated with longitudinal change in FPG due to limited statistical power. In addition, it could be possible that unlike static FPG level that is more affected by genetic variants, its longitudinal change could be more affected by environmental factor.
One of the questions we tried to address during this study was whether genetic variants associated longitudinal change in FPG also affect risk of T2DM or degree of hyperglycemia and vice versa. Although the rs11187850 variant showed nominal association for increased static level of FPG or risk of T2DM, it did not reach genome-wide significance. In addition, we were not able to observe a compelling evidence that previously validated genetic variants of FPG or T2DM affect longitudinal changes in FPG. Even though there were five variants of T2DM and three variants of FPG that were nominally associated with longitudinal change in FPG, they were not significant after adjusting for multiple comparison. Our findings suggest that the biological mechanisms or the genetic determinants of FPG itself or T2DM may be different from that of the longitudinal change in FPG.
When we compare our results with that of the European study which investigated nine cohorts to identify genetic variants association with change in FPG over time, none of the variants with suggestive association (P<5.00×10−6) were replicated in our study (P>0.05) [38]. Furthermore, the SNP rs11187850 discovered in this study was not replicated in Europeans (P=0.769). There may be several reasons for the discrepancy. First, the statistical method applied in the two studies were different. We used LMM for each repeated measure, but in the European study the longitudinal patterns of FPG were summarized into an individual slope and was analyzed using LM. Our research method has the effect of reducing the variance of the error term by considering more measurements. In addition, we adjusted for baseline BMI which might impact longitudinal change in FPG. Second, there could be ethnic difference in genetic impact and their interaction with environment on longitudinal change in FPG. It has been shown that relative contribution of β-cell dysfunction and insulin resistance could be different between East Asians and Europeans [11]. In addition, it has been shown that there are ancestry specific variants, such as rs2233580 in PAX4, that affect T2DM and related phenotypes only in East Asians [39].
The present study has certain limitations. First, as shown in the GCTA analysis, genome-wide SNP interaction with time explained only a limited portion of the total variance of the longitudinal change in FPG. Considering that only one SNP associated with longitudinal change in FPG was detected, and four SNPs associated with FPG were identified among many previously reported variants, the current sample size may be insufficient to elucidate additional SNP interaction with time. Second, the two cohorts used in the analysis had subtle difference in clinical characteristics, regarding age, BMI, FPG at baseline, follow-up time, and incidence of T2DM as the KoGES was a community-based cohort and the GENIE was a hospital health checkup cohort in an urban area. For the GENIE cohort, FPG measurement interval varied individually and the follow-up duration was shorter compared to the KoGES. Nevertheless, the fact that the direction of the effect size of the suggestive variants were all estimated to be the same in both cohorts means that our results are noteworthy. Third, since the genotyping platforms of the two cohorts had different genomic coverage for common variants (84.8% for KoGES and 95.4% for GENIE), there would have been variants that were not captured in both cohorts. If we were able to use the most updated genotype platform for the two cohorts, more significant variants affecting the longitudinal change of FPG might have been identified.
In conclusion, large-scale longitudinal cohorts in Korean population enabled us to investigate loci influencing FPG change over time through GWAS-meta analysis with LMM. Our study identified one novel genetic association between rs11187850 and longitudinal change in FPG. This variant was an eQTL for NOC3L expression in T2DM-related tissues such as pancreas and adipose. Further studies are required to validate our findings and to understand the details of the biological mechanisms underlying this association. Through these efforts, we hope to derive novel insights on preventing progression of hyperglycemia and development of diabetes.
XML Download