Genetic Determinants of Glycated Hemoglobin in Type 1 Diabetes.

Glycated hemoglobin (HbA1c) is an important measure of glycemia in diabetes. HbA1c is influenced by environmental and genetic factors both in people with and in people without diabetes. We performed a genome-wide association study (GWAS) for HbA1c in a Finnish type 1 diabetes (T1D) cohort, FinnDiane. Top results were examined for replication in T1D cohorts DCCT/EDIC, WESDR, CACTI, EDC, and RASS, and a meta-analysis was performed. Three SNPs in high linkage disequilibrium on chromosome 13 near relaxin family peptide receptor 2 (RXFP2) were associated with HbA1c in FinnDiane at genome-wide significance (P < 5 × 10-8). The minor alleles of rs2085277 and rs1360072 were associated with higher HbA1c also in the meta-analysis with RASS (P < 5 × 10-8), where these variants had minor allele frequencies ≥1%. Furthermore, these SNPs were associated with HbA1c in an East Asian population without diabetes (P ≤ 0.013). A weighted genetic risk score created from 55 HbA1c-associated variants from the literature was associated with HbA1c in FinnDiane but explained only a small amount of variation. Understanding the genetic basis of glycemic control and HbA1c may lead to better prevention of diabetes complications.

Introduction

Glycated hemoglobin (HbA1c) reflects the level of blood glycemia and is commonly used to monitor glycemic control in diabetes and to diagnose type 2 diabetes. Genome-wide association studies (GWAS) have revealed >60 genetic variants to be associated with HbA1c in individuals without diabetes (1–6). A recent GWAS of multiple ethnic groups confirmed 18 variants previously associated with HbA1c and added 42 new loci (6). Some of the variants affect biological pathways maintaining glucose homeostasis. For example, single nucleotide polymorphisms (SNPs) in or near glucokinase (GCK) and melatonin receptor 1B (MTNR1B) are known to be associated with HbA1c, fasting plasma glucose, and β-cell function (3). HbA1c is further influenced by genes and pathways that affect erythrocytes. In particular, SNPs in hemochromatosis (HFE) and hexokinase 1 (HK1) (3) have an effect on HbA1c through nonglycemic pathways that affect erythrocyte physiology and survival.

HbA1c reflects the glycemic load over ∼3 months (7), and traits such as hemoglobinopathies, alterations in intracellular glucose metabolism, or defects of glucose transport into erythrocytes affect HbA1c. In addition, diet, drugs, and insulin administration affect HbA1c in individuals with type 1 diabetes (T1D) who have little to no endogenous insulin production and, thus, rely on an external source of insulin to ensure glucose uptake in cells. Despite a strong environmental stimulus, HbA1c levels are correlated between both dizygotic twins and monozygotic twins concordant or discordant for T1D (8). Heritability estimates for HbA1c in individuals without diabetes range from 27% to 62% (8,9), thus providing evidence for a genetic component for the trait.

To date, one GWAS on HbA1c has been conducted in individuals with T1D in the Diabetes Control and Complications Trial (DCCT) (10). The authors found that rs1358030 near sortilin-related VPS10 domain–containing receptor 1 (SORCS1) was associated with HbA1c (P < 2.2 × 10−9). The current study aims to define more genetic variants associated with HbA1c using GWAS of 4,622 individuals with T1D from the Finnish Diabetic Nephropathy (FinnDiane) Study and attempts to replicate GWAS results in five additional cohorts including, altogether, 3,097 individuals with T1D.

Research Design and Methods

Phenotypes

HbA1c

FinnDiane is a nationwide, comprehensive multicenter study launched in 1997. Local ethics committees approved the study protocol, and participants gave informed consent before participation. The study was performed in accordance with the Declaration of Helsinki. HbA1c was measured with standardized laboratory methods locally from a blood sample taken during the visits to the FinnDiane study centers. Additional HbA1c values were collected from laboratory records both before and after the FinnDiane visit. A shift from high-performance liquid chromatography method to immunoassays was seen in many laboratories during the HbA1c collection time (1 January 1991–1 June 2016). In Finland, the HbA1c measurements show a high correlation with the DCCT reference method (11). HbA1c measurements from <1 year of diabetes diagnosis and the measurements after end-stage renal disease (ESRD) diagnosis were excluded. A median of 19 HbA1c measurements per person (range 1–129) were collected (Supplementary Fig. 1). Phenotype in the GWAS was an intrapersonal mean HbA1c, and it roughly followed normal distribution (Supplementary Fig. 2). In the replication cohorts, a mean of natural log-transformed HbA1c was used. Similar phenotypic values were subsequently calculated in FinnDiane to enable comparison of effect sizes (β) in the meta-analysis.

Clinical Parameters

At the FinnDiane visit, each patient underwent a thorough clinical investigation as previously described (12). The GWAS covariates age and diabetes duration were calculated at the intrapersonal mean date of HbA1c measurements. Micro- or macroalbuminuria was defined based on the urinary albumin excretion rate (AER) in two out of three consecutive timed overnight or 24-h urine collections as previously described (12). ESRD was defined as ongoing dialysis or the patient having received a kidney transplant. ESRD diagnosis year was verified from hospital records. Hematological measurements were taken from whole blood samples of 429 individuals using flow cytometry in the Helsinki University Central Hospital Laboratory.

Statistical Methods and Analyses

Patient Selection

Individuals with age at diabetes onset of <40 years, insulin treatment initiated within 1 year of the diabetes diagnosis, and one or more HbA1c measurements were selected from FinnDiane. Relatedness of individuals was calculated with KING 1.40 (13). Among first-degree family members, the individual with the least number of HbA1c measurements (nHbA1c) was excluded, and altogether 4,622 individuals were included in the GWAS.

After also excluding more distantly related individuals according to the estimated genetic relationship matrix (≥0.05), we calculated heritability in 4,319 individuals. Narrow-sense heritability of HbA1c (h2) was the proportion of the phenotypic variance explained by the additive effects of the genotyped SNPs with or without covariates sex, age, duration, genotyping batch, and nHbA1c (GCTA, version 1.24.4) (14).

Power Calculations

Power calculations were performed with Quanto, version 1.2.4 (http://biostats.usc.edu/Quanto.html), using α-level of 5 × 10−8 (genome-wide significance). In the GWAS, we had 94% power to detect variants with a minor allele frequency (MAF) ≥5% and an effect size (β) 0.40 (nontransformed HbA1c [%] [Supplementary Table 1]). In the replication cohort (n = 3,097), we had >57% and 99% power (two-sided α = 0.05, β = 0.30) to replicate variants with MAFs 1–5% or MAFs ≥5%, respectively.

Genotyping and GWAS

SNP genotyping was performed with HumanCoreExome-12 v. 1.0, -12 v. 1.1, or -24 v. 1.0 BeadChip (Illumina, San Diego, CA). SNPs with low genotyping call rate, deviation from Hardy-Weinberg equilibrium, MAF <1%, Mendelian inconsistency, allele frequency difference >20% (MAF ≥5%), or allele frequency difference >5% (MAF <5%) with 1000 Genomes European (EUR) population were filtered.

After quality control, phased haplotypes (SHAPEIT v2.r837) (15) of 316,899 SNPs were used for imputation with 1000 Genomes EUR phase 3, version 5 (16), using Minimac3/Minimac3-omp v1.0.14 (17). A total of 6,019 individuals passed quality control when samples with genotyping rate <0.95 or extreme heterozygosity, sample mix-ups, and genetic outliers were removed.

GWAS was performed with dosages of 8,354,674 variants with imputation quality r2 ≥ 0.7 with PLINK v1.9b3.26 (18) using linear regression with covariates sex, age, diabetes duration, nHbA1c, genotyping batch, and 10 first principal components that were calculated with the EIGENSTRAT, version 3.0 (19). A genomic inflation factor λ was 1.003 (Supplementary Fig. 3).

Altogether, 1,489 and 1,863 individuals were genotyped with the TaqMan technology and ABI Prism 7000/7900HT Sequence Detection Systems and SDS 2.3 software (Life Technologies, Foster City, CA) to verify the imputed genotypes of rs2085277 and rs17076364, respectively. Of the TaqMan-genotyped individuals, 1,259 of 1,489 (84.6%) for rs2085277 and 1,545 of 1,863 (82.9%) for rs17076364 were included in the HbA1c GWAS.

Longitudinal HbA1c

Longitudinal HbA1c data were studied with linear mixed models (nlme package in R, version 3.3.3) in 4,153 individuals with three or more HbA1c measurements. The top SNPs were analyzed with continuous time including diabetes duration (at first HbA1c), sex, age, genotyping batch, and 10 first principal components as covariates assuming random intercept and slope. The best model (lowest Akaike information criterion) had an autocorrelation structure of order 1 with continuous time covariate and a combined covariance structure with exponential variances for fitted model values and separate variances for HbA1c measurements from different decades of diabetes duration.

Phenotypic Comparisons

Groups were compared with Student t test, Welch t test, or Mann-Whitey U test in R. Histograms were used for evaluating normality. Linear regression was used to study phenotypic associations with covariates.

Replication

T1D Cohorts.

Replication of the top SNPs was conducted in DCCT/Epidemiology of Diabetes Interventions and Complications (EDIC), the Wisconsin Epidemiologic Study of Diabetic Retinopathy (WESDR), the Coronary Artery Calcification in Type 1 Diabetes (CACTI) study, the Pittsburgh Epidemiology of Diabetes Complications (EDC) study, and the Renin Angiotensin System Study (RASS) comprising altogether 3,097 individuals with T1D (Fig. 1 and Supplementary Table 2). Individuals with conventional (n = 667) or intensive (n = 637) treatment of diabetes during DCCT were analyzed separately. CACTI was also divided (clinic, n = 265, and nonclinic, n = 253). SNPs with imputation quality r2 ≥ 0.6 and MAF ≥1% in the individual replication cohorts were included in the meta-analysis. These SNPs were replicated in one to seven cohorts comprising 253–3,097 individuals. A mean of natural log-transformed HbA1c was the phenotype used in the replication cohorts, and to get comparable effect sizes (β) the association of the top SNPs and similarly transformed HbA1c was calculated in the discovery cohort. Thereafter, a fixed-effects meta-analysis based on βs and SEs was performed with GWAMA (20).

Figure 1

Outline of the GWAS for HbA1c in the FinnDiane cohort. 1000G, 1000 Genomes.

Populations Without Diabetes and Literature Search.

The top nine independent variants were compared with publicly available data from GWAS meta-analyses of up to 123,517 individuals without diabetes of European ancestry (6). For variants not included in the GWAS by Wheeler et al. (6), we searched for proxy variants (r2 ≥ 0.3) in the 1000 Genomes EUR population using LDlink (21) (https://analysistools.nci.nih.gov/LDlink). Seven out of nine SNPs or proxy variants had HbA1c association data available. In addition, replication of three genome-wide significant SNPs on chromosome 13 was performed in the GWAS data from individuals with East Asian ancestry (n = 23,461), since these SNPs were not available in the Europeans, likely due to their low MAFs. Thereafter, we compared 86 previously reported variants for HbA1c with our GWAS results: 11 variants from the DCCT GWAS (10) and 75 variants from GWAS of individuals without diabetes.

Genetic Risk Scores

Genome-wide significant autosomal SNPs in the largest HbA1c GWAS (6), which had MAF ≥1% and P < 0.01 also in the Europeans, were selected for genetic risk score (GRS). GRS-ALL comprised all 55 SNPs, and like Wheeler et al. (6), we calculated separate GRS for erythrocytic SNPs (n = 19) (GRS-E) and glycemic SNPs (n = 18) (GRS-G). GRS were the sum of expected number of risk alleles weighted by the effect size (β) of the SNP in the Europeans, multiplied by the number of SNPs and divided by the sum of βs. Association of the GRS with HbA1c was studied using linear regression with covariates age, sex, diabetes duration, and nHbA1c.

After the FinnDiane GWAS, a weighted GRS, GRSFD, was created from nine noncorrelated SNPs showing significant or suggestive association with HbA1c (P < 5 × 10−6). GRSFD was tested for associations with diabetes duration, diabetic nephropathy (DN), and laser-treated diabetic retinopathy.

In Silico Analyses

The closest genes were searched (Ensembl transcripts downloaded from https://genome.ucsc.edu/cgi-bin/hg Tables [July 2017]) for variants with P < 0.001 in the GWAS. A total of 9,742 SNPs were mapped to 4,470 unique transcripts and were then connected with the HUGO Gene IDs. A total of 1,175 genes were located <150 kb from the GWAS SNPs. PANTHER database (22) (http://www.pantherdb.org) was used to identify the underlying biological pathways, and STRING, version 10.5 (https://string-db.org), was used to study protein interactions.

For the nearest genes of the top variants, gene expression (Expression Atlas [http://www.ebi.ac.uk/gxa/home]) and protein expression (ProteomicsDB [https://www.proteomicsdb.org]) were studied. The ENCODE data wereaccessed through UCSC Genome Browser (http://www.genome.ucsc.edu) and RegulomeDB v1.1 (23) (http://regulome.stanford.edu/ [July 2017]). Expression quantitative trait loci (eQTL) of the top variants were studied using the Genotype-Tissue Expression project portal (GTEx V7) (https://www.gtexportal.org/home/). PhenoScanner V2 database (24) (http://www.phenoscanner.medschl.cam.ac.uk/phenoscanner) was used to find associations with metabolites.

Results

Mean ± SD HbA1c was 8.5 ± 1.2% (69.5 ± 13.1 mmol/mol) in FinnDiane (Table 1). Narrow-sense heritability of HbA1c was 17.0% (P = 0.002) (unadjusted) or 15.3% (P = 0.003) when adjusted for covariates. Minor alleles of three SNPs on chromosome 13q12.3 showed genome-wide significant association with higher HbA1c values in the GWAS (lead SNP rs2085277: β = 0.42 [95% CI 0.27, 0.56], P = 1.52 × 10−8) (Fig. 2). The three SNPs rs2085277, rs1360072, and rs17076364 are in high linkage disequilibrium (LD) in the 1000 Genomes EUR population and are located between β 3-glucosyltransferase (B3GLCT) (SNPs 232–261 kb downstream) and 146–175 kb upstream of the closest gene, relaxin family peptide receptor 2 (RXFP2) (Supplementary Fig. 4).

Table 1

Baseline demographics of the FinnDiane cohort

n	4,622
Sex (% male)	51.0
Age (years)	38.8 ± 12.4
HbA1c, % (mmol/mol)	8.5 ± 1.2 (69.5 ± 13.1)
Age at diabetes diagnosis (years)	13.9 (8.9–21.6)
Duration of diabetes (years)	21.9 (12.4–31.0)
Number of HbA1c measurements	19 (9–32)
HbA1c collection time (years)*	25.4 (22.4–25.4)
HbA1c measurements (n/year)	0.9 (0.4–1.4)
HbA1c variability, CV (%)	8.9 (6.8–11.9)
BMI (kg/m2)	25.1 ± 3.7
Insulin (IU)	50.7 ± 21.3
Insulin dose (IU/kg)	0.70 ± 0.26
SBP (mmHg)	134 ± 19
DBP (mmHg)	79 ± 10
Triglycerides (mmol/L)	1.02 (0.76–1.45)
DN status, n (%)
Normal AER	2,794 (60.5)
Microalbuminuria	610 (13.2)
Macroalbuminuria	746 (16.1)
ESRD	227 (4.9)
Undefined	245 (5.3)
Retinopathy, n (%)†	1,530 (34.6)
Hard event, n (%)	447 (9.8)

Data are mean ± SD or median (interquartile range) unless otherwise indicated. Coefficient of variation (CV) is defined as CV = SD/mean ×100 and calculated from T1D subjects with three or more HbA1c measurements. Cardiovascular hard events were myocardial infarction, cardiac bypass surgery, stroke, peripheral vascular disease, or amputation at baseline. HbA1c measurements between 1991 or 1 year after diabetes diagnosis year and 1 June 2016 or before the ESRD year were collected. DBP, diastolic blood pressure; SBP, systolic blood pressure.

*HbA1c collection time is calculated as years between the beginning of the HbA1c collection (1 January 1991 or 1 year after diabetes diagnosis year) and the end of the collection (1 June 2016 or ESRD year).

†Laser-treated diabetic retinopathy.

Figure 2

Manhattan plot of the FinnDiane HbA1c GWAS results shows genome-wide significant (P < 5.0 ×10−8) variants on chromosome 13. The lead SNP is rs2085277 (P = 1.52 ×10−8). y-axis shows the –log10 P value, and x-axis tells the chromosomal position of the variants. The two horizontal lines denote the thresholds for genome-wide significant P value of 5 × 10−8 (upper dashed line) and for a suggestive association with P value of 5 ×10−6 (lower dotted line).

Additional genotyping confirmed the imputed rs2085277 and rs17076364 genotypes (Supplementary Table 3). The lead SNP rs2085277 was not significantly associated with insulin dose (IU/kg), coefficient of variation of HbA1c (Supplementary Table 4), or hematocrit or other hematological parameters that were measured in a small subcohort of 429 individuals with T1D (Supplementary Table 5). Altogether, 68 variants had P < 5 × 10−6 in the GWAS. Of these, five tightly correlated SNPs from chr13 locus and eight other independent SNPs with P < 5 × 10−6 were selected for replication and subsequent analyses. Notably, these variants were associated with HbA1c also when longitudinal HbA1c measurements were analyzed in a subcohort of 4,153 individuals with three or more HbA1c measurements (Supplementary Table 6).

Replication

T1D Cohorts

The five T1D cohorts all used the mean of the natural log-transformed HbA1c values as the end phenotype. Two SNPs showed significant replication with ln(HbA1c): rs373883321 in CACTIclinic (chr6:165,403,975, P = 0.006) and rs10459356 near the lead SNP rs2085277 on chromosome 13 in RASS (P = 0.042) (Supplementary Table 7). In a meta-analysis with the discovery cohort, variants rs2085277 and rs1360072 in high LD on chromosome 13 were significantly associated with ln(HbA1c) with β = 0.046 (95% CI 0.030, 0.063) and P = 3.7 × 10−8 and P = 4.2 × 10−8, respectively (Table 2 and Supplementary Table 8).

Table 2

Top loci in the FinnDiane HbA1c GWAS and results from replication in other T1D cohorts

Variant	Position	Closest gene	FinnDiane (n = 4,622)				T1D replication cohorts		Meta-analysis P
			EA/OA	EAF (%)	β (95% CI)*	P	n cohorts; subjects	EAF (%)
rs7514675	1:99,782,478	LPPR4	G/A	10.1	0.23 (0.13, 0.33)	2.7 × 10−6	7; 3,097	12.4	1.5 × 10−2
rs144918527	3:12,740,810	RAF1	T/TTG	53.2	0.13 (0.08, 0.19)	4.4 × 10−6	5; 2,248	51.2	2.5 × 10−4
rs373883321	6:165,403,975	MEAT6	T/–	34.4	0.13 (0.07, 0.18)	3.9 × 10−6	5; 2,248	36.8	2.1 × 10−2
rs478414	7:154,526,318	DPP6	C/T	56.0	0.13 (0.07, 0.18)	3.1 × 10−6	7; 3,097	50.5	8.3 × 10−3
rs113241624	8:109,658,463	TMEM74	A/G	1.2	0.60 (0.35, 0.85)	3.2 × 10−6	1; 596	1.0	8.5 × 10−5
rs4744017	9:93,615,717	SYK	A/G	47.9	0.12 (0.07, 0.17)	8.8 × 10−6	7; 3,097	46.8	8.3 × 10−4
rs2085277	13:32,167,717	RXFP2	T/A	4.1	0.42 (0.27, 0.56)	1.5 × 10−8	1; 253	1.0	3.7 × 10−8
rs552894079	15:89,840,012	FANCI	A/G	97.6	0.40 (0.23, 0.56)	8.1 × 10−6	NA	NA	NA
rs2599441	19:43,989,088	PHLDB3	G/A	72.8	0.14 (0.08, 0.19)	2.0 × 10−6	7; 3,097	76.2	3.2 × 10−3

EA, effect allele; EAF, effect allele frequency; OA, other allele.

*β is the effect size for the mean HbA1c phenotype in the discovery cohort. Replication cohorts are two DCCT/EDIC cohorts (conventional [n = 667] and intensive [n = 637] treatment of HbA1c in the DCCT separately), WESDR (n = 596), RASS (n = 253), CACTI divided into two cohorts (nonclinic [n = 253] and clinic [n = 265]), and EDC (n = 253). Since the mean of natural log-transformed HbA1c was the phenotype used in the T1D replication cohorts, top SNP association with ln(HbA1c) was further calculated in the discovery cohort to get comparable βs and SEs for the meta-analysis. Meta-analysis P value is from fixed-effects meta-analysis of FinnDiane and all replication cohorts restricted to variants with MAF ≥1% in the individual replication cohorts. Results for the variants rs1360072 and rs17076364 on chromosome 13 locus were similar to the top variant, rs2085277: FinnDiane, P = 1.8 × 10−8 and P = 2.6 × 10−8, and meta-analysis, P = 4.2 × 10−8 and P = 6.5 × 10−8, respectively. rs113241624 on chromosome 8 had MAF ≥1% only in WESDR and rs2085277 on chromosome 13 only in RASS. SNP rs552894079 was not available in the T1D replication cohorts.

Of note, the MAFs for rs2085277 and rs1360072 on chromosome 13 were 4% in FinnDiane and 1% in RASS but only 0.2%–0.5% in the DCCT/EDIC and WESDR cohorts that were thus not included in the meta-analysis for these SNPs (Supplementary Fig. 5). Because of low imputation quality of the chr13 top variants (r2 < 0.60), CACTI and EDC cohorts were also excluded from the meta-analysis of these SNPs. Low MAFs reduced the statistical power to replicate the findings, but in RASS the effect was in the same direction as in FinnDiane: lead SNP rs2085277, β = 0.059 [95% CI –0.079, 0.198] for ln(HbA1c), P = 0.402. The MAF for rs2085277 is similar to our study in the 1000 Genomes populations: 5.6% in FIN (Finnish in Finland), 0.5% in GBR (British in England and Scotland), and 0% (monomorphic) in the CEU (Northern Europeans from Utah). The variant showed a markedly higher frequency in East Asians and South Asians with MAFs of 38.5% and 9.2%, respectively.

Cohorts Without Diabetes

The association of the thirteen SNPs with HbA1c was then queried in GWAS meta-analysis data on multiple cohorts without diabetes. As expected due to low MAF in Europeans, no data were available for the chr13 top SNPs or for any proxy SNPs with r2 ≥ 0.3 in the European cohorts. Therefore, replication of the chr13 variants was sought in the meta-analysis results of individuals of Asian ancestry. The SNPs were missing from the South Asian cohorts, but the chr13 top SNPs rs2085277 and rs1360072 were associated with HbA1c (%) in East Asians with β = 0.02 (95% CI 0.006, 0.03), P = 0.005, and β = 0.02 (0.003, 0.03), P = 0.013, respectively. (Table 3). None of the other top SNPs or the proxy variants were associated with HbA1c in the Europeans (Supplementary Table 9).

Table 3

Three genome-wide significant variants on chromosome 13: LD structures and HbA1c association in a population without diabetes

Variant	EA < OA	MAF (%)*				r2†				HbA1c association in a Finnish T1D cohort		HbA1c association in an East Asian population without diabetes
		FD	EAS	SAS	EUR	FD	EAS	SAS	EUR	β (95% CI)	P	β (95% CI)	P
rs2085277	T < A	4.2	38.5	9.2	1.2					0.42 (0.27, 0.56)	1.5 × 10−8	0.02 (0.006, 0.03)	0.005
rs1360072	C < T	4.1	34.6	9.0	1.2	1.0	0.75	0.95	1.0	0.41 (0.27, 0.56)	1.8 × 10−8	0.02 (0.003, 0.03)	0.013
rs17076364	T < C	3.1	2.5	0.2	1.2	0.97	0.04	0.02	1.0	0.45 (0.29, 0.61)	2.6 × 10−8	NA	NA
n		4,622	504	503	489	4,622	504	503	489	4,622		23,461

EA is the effect allele in the FinnDiane GWAS and OA is the other allele.

*MAF frequency in the FinnDiane cohort (FD) and in East Asian (EAS), South Asian (SAS), and EUR populations in the 1000 Genomes phase 3 data.

†LD (r2) with rs2085277 in the FinnDiane and in the 1000 Genomes phase 3 populations. Position of the SNPs in the human genome build GRCh37 (hg19) was rs2085277, chr13:32,167,717; rs1360072, chr13:32,138,331; and rs17076364: chr13:32,163,399. HbA1c association statistics for these variants were not available for the European population without diabetes (n = 123,517) or South Asian population without diabetes (n = 9,908) (Wheeler et al. [6]).

In Silico Analyses

Interestingly, minor T allele for rs2085277 was associated with higher levels of ketone bodies acetoacetate and 3-hydroxybutyrate (25) (PhenoScanner) (P = 0.003 and P = 0.020, respectively). No direct ENCODE data were available for rs2085277, and only minimal transcription factor–binding evidence (RegulomeDB scores 4–6) was seen for other top SNPs in the chr13 locus. Although none of the chr13 top variants were included in the GTEx database, many SNPs near rs2085277 affect the expression of RXFP2 in adrenal gland, which suggests that rs2085277 might also affect RXFP2 expression.

RXFP2 belongs to the positive regulation of cAMP biosynthetic process pathway (Gene Ontology term in the STRING protein network analysis [Fig. 3]). Interestingly, when the search was restricted to databases as opposed to all interaction sources (databases, text mining, experiments, coexpression, neighborhood, gene fusion, and co-occurrence), gastric inhibitory polypeptide receptor (GIPR) and glucagon (GCG)—proteins relevant for glucose homeostasis and diabetes—occurred in the same protein interaction network with RXFP2.

Figure 3

Protein networks around RXFP2 (red). A: Two most enriched Gene Ontology (GO) biological pathways were “cellular response to glucagon stimulus” (GNB1, GNB4, GNB2, and GNB3) and “positive regulation of cAMP biosynthetic process” (RXFP2, RLN2 [relaxin], AVPR2, and INSL3), with false discovery rates 4.6 ×10−5 and 2.8 ×10−4, respectively. Protein network was generated with STRING, version 10.5, selecting all active interaction sources, minimum required interaction score 0.90 (highest confidence), and no more than 10 interactions. B: Only database interaction sources were selected with minimum required interaction score 0.90 (highest confidence) and no more than 10 interactions. Six out of 11 proteins in this network belong to “positive regulation of cAMP biosynthetic process” (GO biological process), which was the most enriched pathway in this network, connecting proteins RXFP2, GIPR, GCG, ADCYAP1, MC1R, and FSHR (false discovery rate 1.3 ×10−8). Of note is that gastric inhibitory polypeptide receptor (GIPR) and glucagon (GCG) present in this protein network have a high relevance in glucose homeostasis and diabetes.

Genes involved in cell-cell adhesion processes and developmental processes were overrepresented in GWAS results (PANTHER, P < 0.05, corrected for multiple testing [Supplementary Table 10]). PANTHER pathway analysis pointed out cadhenin signaling, endogenous cannabinoid signaling, and Wnt signaling pathways (P < 0.05, corrected for 160 independent classes).

Duration of Diabetes and Chronic Diabetes Complications

Our GWAS top variants had a consistent effect on HbA1c irrespective of diabetes duration. In all 10-year diabetes duration bins, the GRSFD (Supplementary Table 11 and Supplementary Fig. 6) was significantly associated with HbA1c (Supplementary Table 12). We confirmed this constant association between GRSFD and HbA1c in a small group of individuals (n = 271) having HbA1c measurements from three consecutive diabetes duration bins: 1–9.9, 10–19.9, and 20–29.9 years.

Subjects with T1D with DN had higher mean HbA1c compared with individuals with normal AER (8.9 ± 1.3% vs. 8.2 ± 1.0%, P < 10−16). Similarly, individuals with diabetic retinopathy had higher HbA1c compared with non–laser-treated subjects with at least 15 years of diabetes duration (8.8 ± 1.2% vs. 8.4 ± 1.2%, P < 10−16). As expected, individuals with DN had higher GRSFD (6.5 vs. 6.4, P = 0.038), while the subjects with diabetic retinopathy did not have higher GRSFD (P = 0.205).

Previous HbA1c Studies

rs1358030 was previously associated with HbA1c in the GWAS of 1,307 individuals with T1D (10). However, this association was not replicated in FinnDiane, although the MAFs were similar (MAFDCCT = 36% vs. MAFFinnDiane = 40%), imputation quality was r2 = 0.99, and the effect was in the same direction (β = 0.034 [95% CI –0.02, 0.09], P = 0.192) (Supplementary Table 13).

Of the SNPs identified in populations without diabetes, four variants were associated with HbA1c in FinnDiane: the minor A allele of the missense variant rs1800562 in HFE was associated with lower HbA1c (P = 0.007) with higher effects size than in the population without diabetes (Fig. 4). The HbA1c values for the genotypes AA (n = 3), AG (n = 250), and GG (n = 4,369) were 7.6 ± 1.1%, 8.3 ± 1.2%, and 8.5 ± 1.2%, respectively. Variants rs906216 and rs7072268 in HK1 and rs12621844 near FOXN2 were also directionally consistently associated with HbA1c in FinnDiane (P < 0.05) (Supplementary Tables 13 and 14). Of note, 16 of 19 SNPs classified as erythrocytic variants by Wheeler et al. (6) showed similar direction of effect in FinnDiane compared with the original study (binomial test for enrichment, P = 0.002 [Fig. 4]).

Figure 4

Effect sizes (βs) of 55 variants associated with HbA1c (%) in the original study by Wheeler et al. (6) (European population without diabetes on x-axis) and in FinnDiane (T1D cohort on y-axis). The SNPs were classified into three groups (erythrocytic SNPs, glycemic SNPs, or unclassified SNPs) in the study by Wheeler et al. based on the literature search of known associations of the SNPs in earlier published GWAS for glycemic or erythrocytic traits. Horizontal dashed line denotes the effect size zero (β = 0) in FinnDiane indicating that variants with β > 0 have the effect on HbA1c in the same direction as in the population without diabetes. SNPs above the diagonal dotted line have a stronger effect on HbA1c in the T1D cohort. Two SNPs with HbA1c association P < 0.05 in FinnDiane are indicated by arrows. Gene names are shown for these variants and for HK1 because this particular SNP (rs4745982) in HK1 was not associated with HbA1c in FinnDiane. Instead, two other SNPs in HK1 (rs906216 and rs7072268 from Paré et al. [1]) showed significant association with HbA1c in FinnDiane. FOXN2, forkhead box protein N2; HFE, hemochromatosis; HK1, hexokinase.

GRS

Next, we tested the combined genetic effect by creating GRS from significant variants found in the largest HbA1c GWAS (6) (Supplementary Fig. 7). A GRS-ALL comprising 55 SNPs was associated with HbA1c in FinnDiane (β = 0.010 [95% CI 0.002, 0.018], P = 0.014) and, along with covariates, explained 2.5% of the variation in HbA1c. However, age, diabetes duration, sex, and nHbA1c alone explained 2.4% of the variation, whereas GRS-ALL explained only 0.1% of the variation in HbA1c. A similar amount of HbA1c variance was also explained by GRS-E including only 19 erythrocytic SNPs (β = 0.021 [95% CI 0.005, 0.037], P = 0.01). In contrast, GRS-G created from the 18 glycemic variants was not associated with HbA1c (β = 0.004 [95% CI –0.01, 0.018], P = 0.573).

Discussion

We identified a locus on chromosome 13 that show genome-wide significant association with HbA1c in 4,622 Finnish individuals with T1D. The three associated SNPs in the locus are located <30 kb from each other and are in high LD, therefore indicating the same association signal. Notably, the association with HbA1c at rs2085277 and rs1360072 remained genome-wide significant in a meta-analysis including individuals with T1D from RASS. In DCCT/EDIC and WESDR, the MAFs for the chr13 top SNPs were ≤0.5%. The low MAFs (≤1%) might explain the lack of significant replication in the T1D cohorts, since the power to detect significant association with MAF ≤1% is limited. On the contrary, both SNP associations were replicated (P ≤ 0.013) in the East Asian population without diabetes, where the SNPs were present with high allelic frequency of 38.5% and 34.5% for rs2085277 and rs1360072, respectively.

rs2085277 is an intergenic SNP surrounded by B3GLCT (261 kb 5′) and RXFP2 (146 kb 3′). More distant B3GLCT is a glucosyltransferase that adds glucose molecules to multiple proteins. RXFP2 seems a more plausible candidate gene, since many SNPs in close vicinity to rs2085277 (although not in LD) have cis-eQTLs for RXFP2 in adrenal gland. RXFP2 is a receptor for insulin-like peptide 3 (INSL3) and relaxin, and RXFP2 is expressed at low levels in many tissues. ProteomicsDB shows RXFP2 expression in adrenal glands and in pancreatic islets. Although the pharmacology and biology of relaxin receptors are relatively well-known (26), we are not aware of any previous studies on RXFP2 in relation to glucose homeostasis. Of note, protein interaction analysis restricted to databases shows that RXFP2 is linked to the glucagon response pathway and to relevant proteins for glucose homeostasis, like GIPR and GCG. Interestingly, adenylate cyclase 5 (ADCY5), which has previously been linked to HbA1c and T2D (6,27), is found in the same protein interaction network. However, there is yet no experimental evidence that would link these proteins. Thus, the biological mechanisms by which RXFP2 possibly affects HbA1c remain to be elucidated.

RXFP2 stimulates adenylate cyclase to increase intracellular second messenger cAMP concentrations. The main ligand for RXFP2 is relaxin hormone, which is naturally increased in pregnancy, and it has both vasodilative and antifibrotic actions (26). Interestingly, relaxin might act along with insulin to control blood glucose homeostasis (28). This was also suggested by a study in mice, where relaxin treatment ameliorated the insulin sensitivity of mice on a high-fat diet (29). The main receptor for human relaxin is relaxin family peptide receptor 1 (RXFP1) but, to a smaller extent, relaxin is also a ligand for RXFP2.

Our main phenotype was mean HbA1c, which was based on a median of 19 measurements per individual. One limitation of the study is that some individuals had only one HbA1c measurement. A single HbA1c measurement is not optimal for T1D studies, since HbA1c values tend to vary within a person over time. When a smaller subcohort of FinnDiane patients with three or more HbA1c measurements was studied with another method, linear mixed models, all top variants were consistently associated with HbA1c, although the effect sizes tended to be smaller. This was maybe due to a smaller cohort (n = 4,153) and deviation from the linear trend over time (individual slopes) in the HbA1c values.

In diabetes, HbA1c values are used as an indicator of glycemic control achieved with insulin treatment or with drugs improving insulin sensitivity, the latter mainly used in patients with type 2 diabetes. The heritability estimate for HbA1c in our study was modest but significant (h2 = 15%) and confirmed that HbA1c in T1D can hold a heritable component (30). In the HbA1c GWAS in DCCT, a variant near SORCS1 was found. We were, however, unable to replicate the finding, although the effect was in the same direction. The lack of replication may be due to the winner’s curse, which causes the effect in the discovery cohort to be overestimated, and therefore the power for replicating the true smaller genetic effect is lower. Here, we had 100% power to replicate rs1358030 with a β = 0.045 [for ln(HbA1c)] and only 62% power with a β = 0.003 corrected for winner’s curse with a bootstrap method at two-sided α = 0.05 (31). On the other hand, four variants in or near HFE, HK1, and FOXN2 originally found in cohorts without diabetes were also associated with HbA1c (P < 0.05) in FinnDiane.

Variants associated with HbA1c have mostly been discovered in large GWAS meta-analyses and have relatively small effect sizes. Therefore, the power to replicate these variants in our study was only ≤33%. However, a GRS comprising 55 SNPs associated with HbA1c in our T1D cohort, although it explained only a small part of the variation of HbA1c. Interestingly, when divided into erythrocytic and glycemic risk scores, GRS-E including only 19 variants classified as erythrocytic SNPs was also associated with HbA1c, while no association was seen for the glycemic risk score. This is expected, as the erythrocyte biology can be expected to be similar in individuals with and without diabetes, but the mechanisms affecting the level of glycemia are presumably different in individuals with T1D compared with the individuals without diabetes. As discussed by Wheeler et al. (6), the genetic factors primarily affecting erythrocyte biology may result in falsely high or low HbA1c levels no longer fully reflecting the level of glycemia and, thus, in false interpretation of HbA1c.

Allele frequencies of risk variants may differ between populations, and therefore it is important to conduct genetic analyses in different populations. Finns represent a population that has been genetically isolated from the rest of Europe, thus presenting an interesting opportunity to explore genetic variants affecting HbA1c in T1D. Indeed, the lead SNP was observed with a markedly higher MAF in the Finnish than in non-Finnish Europeans.

To conclude, we identified a locus on chromosome 13 near RXFP2 that is associated with HbA1c in individuals with T1D. Further studies are needed to find out whether relaxin hormone mediates the effect of RXFP2 on HbA1c. Since effect sizes are expected to be smaller in the population without diabetes, additional replication of the top locus as well as the search for additional genetic factors affecting HbA1c in other T1D cohorts is warranted. As HbA1c has a major influence on the development of diabetes complications, it is important to continue the search for genetic variants affecting HbA1c.