figure a

Introduction

CREBRF (encoding cAMP response element binding protein 3 [CREB3] regulatory factor) was recently identified as an obesity gene in a genome-wide association study (GWAS) in Samoans, a Polynesian population [1]. The discovery variant, rs12513649, was strongly associated with BMI. In Samoans, rs12513649 was in virtually perfect concordance with rs373863828 (G>A), at which the A allele (frequency = 25.9%), which is associated with higher BMI and decreased risk of diabetes, codes for an arginine to glutamine amino acid substitution in CREBRF (p.Arg457Gln or R457Q). A study in Polynesians from New Zealand replicated these associations with BMI and diabetes, with frequencies of the A allele ranging from 9.6% to 23.4% [2]. The BMI association with these variants was also replicated in Tongans, another Polynesian population with an A allele frequency of 15% [3]. A small study found frequency of the A allele to be 0–6% in other Oceanic populations (Melanesians and Micronesians) but its association with BMI was not statistically significant [3]. In the present study, we examined the associations of rs12513649 and rs373863828 with BMI and diabetes in Pacific Islanders, predominantly Mariana Islanders and Micronesians, from Guam and Saipan.

Methods

Participants

The study included 2022 participants from the Micronesian extension of the Family Investigation of Nephropathy and Diabetes, designed to identify genetic determinants of diabetes and kidney disease [4]. Participants were ≥18 years old and were recruited in 2008 and 2009 from dialysis centres, referrals of family members and general advertising. Participants were predominantly Chamorros, indigenous inhabitants of the Mariana Islands, but also included Filipinos and several other Micronesian populations. Fasting plasma glucose was measured in participants without known diabetes. Participants were defined as having diabetes if they had fasting plasma glucose ≥7.0 mmol/l, HbA1c ≥6.5% (44 mmol/mol) or a previous diagnosis of the condition [5]. Fasting serum insulin was measured and measures of insulin resistance and beta cell function were calculated according to the homeostatic model (HOMA-IR and HOMA-B) [6]. Participants reported their maximum lifetime body weight and height at this weight and the corresponding BMI was calculated. Characteristics of participants are shown in electronic supplementary material (ESM) Table 1.

Genotypes

IndexCREBRFvariants rs12513649 and rs373863828 were typed individually by KASP on Demand according to manufacturer’s instructions (LGC, Teddington, England).

Additional CREBRF variants

We analysed 54 additional variants within 200 kb of the transcription start/stop sites for CREBRF. We generated the genotypes as part of a GWAS using the Illumina OmniExpress CoreExome array (Illumina, San Diego, CA, USA).

Genomic control variants

The GWAS-derived genotypes were used for a statistical genomic control procedure (see below). These variants included 269,234 SNPs with minor allele frequency >1%, call rates >95%, discordance rates <2.5% in 100 duplicate specimens, and genotypic distributions consistent with the Hardy–Weinberg equilibrium (p > 0.0001). In addition, whole genome sequence data (BGI, Shenzhen, Guangdong, China) were obtained for a subset of 40 full-heritage Chamorros.

Statistics

Missing genotypes were imputed from phased data of all typed variants with Beagle (http://faculty.washington.edu/browning/beagle/beagle.html; version 3.3.2, accessed 30 May 2013) [7]. Statistical analyses were performed with SAS (SAS Institute, Cary, NC, USA). Associations between genotype and continuous variables (BMI, height, weight, HOMA) were assessed by linear regression models in which phenotype was regressed against number of risk alleles. Covariates included age, sex, end-stage renal disease (ESRD) and the first four principal components from the GWAS, which were included to account for potential population stratification (ESM Fig. 1). To account for residual population stratification, a genomic control procedure was applied across all GWAS markers with a separate inflation parameter for each trait [8]. Association with diabetes was analysed by logistic regression with the same covariates (ESRD was included in the analyses presented to account for the study design; results were similar if it was not included.) A fixed-effects meta-analysis for associations with BMI and diabetes was performed to combine results for the present study with those reported in other Oceanic populations [1,2,3]. Effect sizes were weighted by the inverse of their variance, and heterogeneity was assessed using Cochran’s Q statistic and the I2 measure [9].

Results

Allele frequencies

Frequencies of alleles at rs373863828 and rs12513649 and their haplotypes are shown in Table 1. There was strong linkage disequilibrium between the A allele at rs373863828 and the G allele at rs12513649, as indicated by high D′ values. Since allele frequencies differed, allelic concordance was moderate, with an overall r2 of 0.48.

Table 1 Allele and haplotype frequencies at rs12513649 and rs373863828 in CREBRF among Pacific Islanders from Guam and Saipan

Associations with BMI and diabetes

Associations of rs373863828 and rs12513649 with BMI and diabetes are shown in Fig. 1 a–d. The G allele at rs12513649 was associated with higher BMI (β = 1.55 kg/m2 per copy; p = 0.0026) as was the A allele at rs373863828 (β = 1.48 kg/m2; p = 0.033). These obesity-associated alleles were associated with lower prevalence of diabetes: the OR per copy of the G allele at rs12513649 was 0.63 (p = 0.0063), while the OR for the A allele at rs373863828 was 0.49 (p = 0.0022). To further account for potential population stratification, we additionally adjusted for self-reported primary ethnicity with similar results (ESM Table 2, which also shows results according to primary ethnicity). There was a modest association of the A allele at rs373863828 with taller height (p = 0.0097, ESM Fig. 2); however, both rs12513649 and rs373863828 were significantly associated with weight after adjustment for height (ESM Fig. 3). Higher BMI was strongly associated with higher prevalence of diabetes (ESM Fig. 4), and there was no significant difference in this relationship when rs373863828 A allele carriers were compared with non-carriers (p = 0.81 for interaction).

Fig. 1
figure 1

(a) ‘Bean’ plot of the distribution of BMI according to genotypes at rs12513649. Grey horizonal lines represent individual values, with the length of the line corresponding to the number of observations at each level. The symmetrical plots represent the density at each BMI level, estimated using the kernel density estimation function in SAS. Black horizontal bars represent the mean value. β = 1.55 kg/m2 per copy of the G allele; p = 0.0026. (b) ‘Bean’ plot of the distribution of BMI according to genotypes at rs373863828. β = 1.48 kg/m2 per copy of the A allele; p = 0.033. (c) Prevalence of diabetes according to genotype at rs12513649. OR 0.63 per copy of the G allele; p = 0.0063. (d) Prevalence of diabetes according to genotype at rs373863828. OR 0.49 per copy of the A allele; p = 0.0022. (e) Meta-analysis of the association of the A allele at rs373863828 with BMI, including data from the present study (labelled ‘Guam/Saipan’). Data are presented as the regression coefficient (β, kg/m2 per copy of the A allele) with 95% CI. (f) Meta-analysis of the association of the A allele at rs373863828 with diabetes. Data are presented as OR per copy of the A allele with 95% CI. NZ, New Zealand

By design, a large proportion of individuals had ESRD. The analyses were repeated in the 1682 individuals without ESRD to exclude undue influence from ESRD. The associations of rs12513649 and rs373863828 with BMI and diabetes remained significant after this exclusion (ESM Table 3). Similarly, associations with BMI remained significant when analysed in the 1065 individuals who did not have diabetes or ESRD (ESM Table 3). In the 987 non-diabetic individuals who did not have ESRD, HOMA-B was significantly higher in carriers of the A allele at rs373863828 (by 20% per copy; p = 0.031), but HOMA-IR was not (p = 0.12, ESM Fig. 5); the association with HOMA-B was attenuated with adjustment for BMI (p = 0.22, ESM Fig. 5).

Associations with additional CREBRF variants

In the analysis of all 56 GWAS variants mapping within 200 kb of the transcription start/stop site, rs12513649 achieved the strongest association with BMI (ESM Fig. 6, ESM Table 4), and rs373863828 had the strongest association with diabetes (ESM Fig. 7, ESM Table 4). In models where the effects of rs12513649 and rs373863828 were each conditioned on the effects of the other, rs12513649 remained significantly associated with BMI (p = 0.033, data not shown) while rs373863828 did not (p = 0.97, ESM Fig. 6), and rs373863828 was suggestively, but non-significantly, associated with diabetes (p = 0.080, data not shown) while rs12513649 was not (p = 0.45, ESM Fig. 7). In addition, other than rs373863828, missense variants in CREBRF were not detected in whole-genome sequence data from 40 Chamorros.

Meta-analyses

Meta-analyses combining results of associations with rs373863828 in the present study with those from previous publications are shown in Fig. 1 e–f. The A allele was associated with higher BMI (combined β = 1.38 kg/m2 per allele; p = 2.5 × 10−29) and there was no evidence for heterogeneity across studies (I2 = 0%, p = 0.58). The A allele was also associated with lower risk of diabetes, and there was little evidence for heterogeneity (I2 = 31%, p = 0.23). Although none of the studies investigating the association between rs373863828 and diabetes individually achieved genome-wide statistical significance, the combined effect was significant (OR 0.65, p = 1.5 × 10−13).

Discussion

CREBRF variants are associated with increased BMI in Polynesian populations [1,2,3, 10]. The frequency of the A allele at rs373863828 in Polynesians is 9.6–25.9%, whereas this allele is exceedingly rare in most other global populations (e.g. the frequency was 0.004% in the Exome Aggregation Project; https://www.ncbi.nlm.nih.gov/snp/, accessed 3 November 2018). In Polynesians, the A allele was associated with higher BMI but lower risk of diabetes [1, 2]. In the present study, we extend these findings and show that the A allele is present at low frequencies, between 1.1% and 5.4%, in Marianas and Micronesian populations (Chamorros, Chuukese, Carolinians, Palauans), in whom it is also associated with higher BMI and lower risk of diabetes.

The association with BMI was initially identified with rs12513649, which is in nearly perfect concordance with rs373863828 in Polynesian populations [1]. In the current sample, the concordance is less strong, which, theoretically, allows statistical analyses to separate effects of each SNP. However, results of the present conditional analyses are equivocal and larger sample sizes are required to confidently identify independent effects. Functional studies in murine adipocytes found that the A allele at rs373863828 results in increased lipid accumulation and decreased cellular energy expenditure [1]. The functionality of rs12513649 is not clear, but scores derived from the Genome-Wide Annotation of Variants program (GWAVA; https://www.sanger.ac.uk/sanger/StatGen_Gwava, accessed 5 April 2019), which assesses potential functionality from a variety of genomic annotation features, including epigenetic marks and conservation, are 0.27–0.41; this is only modestly suggestive that it is functional [11]. Given this, and the fact that rs12513649 is not associated with BMI in GWAS from East Asian populations [12], in whom the G allele has a frequency of 6.3% (https://www.ncbi.nlm.nih.gov/snp/, accessed 3 November 2018), it seems likely that rs373863828 is the functional variant. The present study did not find significant evidence for additional variants in CREBRF that influence BMI or diabetes risk.

The mechanisms behind the ‘discordant’ associations, in which the allele associated with higher BMI is associated with lower risk for diabetes, require further investigation. The present findings suggest that the A allele at rs373863828 does not modify the effect of BMI on diabetes risk but that it separately lowers diabetes risk while increasing BMI. The association of the A allele with higher HOMA-B suggests that it may result in greater capacity for compensatory insulin secretion, but more detailed physiological studies are required. Regardless of these mechanisms, the study findings demonstrate that genetic variants in CREBRF are important determinants of obesity and diabetes risk in Oceanic populations.