Journal of Clinical Diabetes
Open Access

Diabetic nephropathy is the most common cause of chronic kidney disease (CKD) in developed countries [1]. The clinical characteristics do not fully predict development of nephropathy in diabetes patients. Epidemiological findings suggested genetic background plays an important role in the development of this renal disease [2,3].

There have been several genome-wide association studies (GWAS). In a recent meta-analysis by Pattaro et al. on associations of estimated glomerular filtration rate (eGFR) based on serum creatinine (Scr), cystatin C, and CKD (defined as eGFR based on Scr <60 ml/min/1.73 m2) with about 2.5 million autosomal single-nucleotide variations (SNVs) in 133,413 individuals of European ancestry in the stage 1 discovery analysis, 29 previously identified loci were confirmed and 48 independent novel loci newly identified [4]. In their trans-ethnic analyses in 42,296 Asian participants, 7 out of 24 newly identified loci with eGFR based on Scr achieved direction-consistent significance. If the GWAS discovery analysis had been started in an Asian population, some different SNVs may have been found. Here we found 77 SNVs upstream of NBEA that affected renal function in a Japanese population with type 2 diabetes mellitus (T2D).

Study population

A cross-sectional GWAS was performed in participants at ages from 35 to 69 years in the Japan Multi-Institutional Collaborative Cohort (JMICC) study. The participants in the J-MICC study were invited from 12 different locations in Japan between 2004 and 2013 (Number of participants: 14,539). The J-MICC study is a cohort study started in 2005 to examine interactions between gene and environment factors in lifestyle-related diseases. Minutes of the J-MICC Study were published elsewhere [5,6]. Briefly, participants answered a questionnaire regarding lifestyle and medical information, and gave a blood specimen at the time of the baseline study. The J-MICC study participants encompasses community citizens, first-visit patients to a cancer institution and health checkup examinees. Written informed consent was obtained from all participants in this study, and the study protocol was approved by the Ethics Committees of Aichi Cancer Center, Nagoya University Graduate School of Medicine, and the other participating institutions in the J-MICC study. Our present study was conducted in keeping with to the principles expressed in the World Medical Association Declaration of Helsinki.

Out of 14,539 participants, 448 were excluded by the GWAS screening described below. Of the remaining 14,091 participants, we identified 1,109 with T2D by the criteria described below. Of the 1,109 participants with T2D, 154 were excluded due to missing data. As a result, we analyzed the data of 955 participants with T2D.

Questionnaire and measurements

The J-MICC study questionnaire consisted of questions on medical history, weight, height, smoking and alcohol drinking habits. Measurements of participants’ body and blood drawing were done as part of the health examination or for research purposes at the J-MICC study participating institutions [5]. Weight and height were measured to the nearest 0.1 kg and 0.1 cm, respectively. Body mass index (BMI) was calculated by dividing body weight in kilograms by the square of height in meters. The HbA1c percentage was obtained with use of a latex aggregation immunoassay (Japan Diabetes Society [JDS] value). The HbA1c percentage was estimated as the National Glycohemoglobin Standardization Program (NGSP) equivalent percentage obtained according to the following formula: HbA1c (NGSP (%))=1.02 × HbA1c (JDS [%])+0.25% [7]. We defined T2D as a fasting blood glucose concentration ≥ 126 mg/dL, or ≥ 200 mg/dL if less than 8 h after meals, or NGSP HbA1c ≥ 47 mmol/mol (6.5%), or participants diagnosed as having DM. Scr was measured enzymatically. The eGFR of each participant was obtained based on SCr, age, and sex with use of the following Japanese eGFR equation proposed by the Japanese Society of Nephrology: eGFR (mL/min/1.73 m²)=194 × Scr (mg/dL)-1.094 × age-0.287 (× 0.739 if female) [8]. Semi-quantified dip-stick urine protein data were available in a part of the patients (N=516) that were classified as 0, 15, 30, 100, or ≥ 250 mg/dL.

Genotyping and quality control filtering

Buffy coat fractions were prepared from blood specimen and stored at −80°C at the central office for J-MICC Study. DNA was prepared from all buffy coat fractions with use of a BioRobot M48 Workstation (Qiagen Group, Tokyo, Japan) at the central study office. For the specimen from two locations (Fukuoka and Kyushu-KOPS), DNA was prepared locally from specimen of whole blood with use of an automatic nucleic acid isolation system (NA-3000, Kurabo, Co., Ltd, Osaka, Japan). The 14,539 study participants from the 12 locations of the J-MICC study were genotyped at RIKEN Center for Integrative Medicine Sciences with use of a HumanOmniExpressExome-8 v1.2 BeadChip array (Illumina Inc., San Diego, CA, USA). Inconsistent sex information in 26 participants between an estimate from genotyping and the questionnaire was detected, and these participants were excluded. The identity-by-descent method supplemented in the PLINK 1.9 software [9,10] identified 388 close relationship pairs (pihat> 0.1875) and one sample from each pair of the 388 was deleted. Principal component analysis (PCA) [11,12] with a 1000 Genomes reference panel (phase 3) [13,14] found 34 participants whose estimated ancestries were outside the Japanese population [15]. These 34 participants were excluded. In the remaining 14,091 participants, SNVs with a Hardy-Weinberg equilibrium exact test P-value <1 × 10-6, and/or a genotype call rate <0.98, a departure from the allele frequency computed from the 1000 Genomes Phase 3 EAS samples, or a low minor allele frequency (MAF) <0.01, were excluded. These filtering resulted in 14,091 individuals and 570,162 SNVs.

Genotype imputation

Genotype imputation was conducted with use of SHAPEIT [16], followed by Minimac3 [17] software (with the 1000 Genomes phase 3 as the reference panel [13]). Strict quality control filters were applied after genotype imputation; i.e., variants with an R2<0.3 were excluded, resulting in 12,617,547 variants. Finally 4, 233,797 variants with MAF <0.01 in T2D patients were removed, resulting in 8,383,750 variants for the analysis. We used Dosage Convertor software [18] to convert dosage files in VCF format from Minimac3 to PLINK formats.

Association analyses between genetic variants and eGFR in T2D

Associations between all imputed variants and eGFR were analyzed with adjustments for age and sex with use of PLINK 1.9 software. Variants achieving genome-wide significance (P<5 × 10^-8) were regarded to be eGFR-associated variants. An R package for creating QQ plot, GWASTools was used [19]. For scatter plots of P values derived from genome-wide scan results for eGFR in T2D, a software, Haploview was used [20]. To visualize regions of interest, we used the LocusZoom program [21].

We performed replication analysis for previously-reported SNVs with use of the J-MICC samples [4,22].

Confounding factor adjustment

In 519 patients with urine protein data, the following additional association analyses were performed. Model 1: associations between imputed variants in chromosome 13 and eGFR were tested with adjustment for age and sex with use of PLINK 1.9 software. Model 2: analysis with model 1 variables plus urine protein data. Model 3: model 1 variables plus PCA components 1 to 3, which were obtained in 14,091 participants with the use of Eigensoft 6.0.1 [10,11]. This analysis was performed since the structure of Japanese population was reported as not homogenous [15]. Model 4: model 1 variable plus BMI, alcohol consumption (g/day), and smoking (current, or past smoking vs. never smoking).

Functional annotations

Genomic locations of variants found in this study were examined based on the UCSC [23] and Ensembl [24] genome browsers. CiseQTL pairs of variants and genes were obtained from the GTEx [25].

Baseline characteristics

Baseline characteristics of the participants with T2D are shown in Table 1. The mean age of the participants was 59.1 ± 7.4 y and the percentage of women and prevalence of current smokers were 34% and 25%, respectively. Mean of BMI, Scr concentration, HbA1c, eGFR and alcohol consumption were 24.8 ± 4.1 kg/m2, 0.76 ± 0.30 mg/dL, 55 ± 11 mmol/mol (7.2 ± 1.5%), 89.3 ± 27.1 ml/min/1.73 m2, and 17.7 ± 27.5 g/day, respectively. Baseline characteristics of the participants with T2D who had urine protein data are also shown in Table 2. The mean age of the participants was 58.5 ± 7.6 y and the percentage of women and prevalence of current smokers were 33% and 24%, respectively. Mean of BMI, Scr concentration, HbA1c, eGFR, alcohol consumption and urinary protein were 25.4 ± 4.0 kg/m2, 0.76 ± 0.22 mg/dL, 53 ± 10 mmol/mol (7.1 ± 1.4%), 88.8 ± 29.5 ml/min/1.73 m2, 18.4 ± 29.5 g/day and 10.3 ± 34.1 mg/dL, respectively.

Table 1: Background characteristics of the study participants.

Table 2: SNVs associated with eGFR in T2.

Genome-wide association study

Among the 8,383,750 variants, genome-wide analyses with adjustment for age and sex, we found 77 SNVs upstream of the NBEA gene that were significantly associated with eGFR in participants with T2D with P values <5 × 10^-8. A representative 8 SNVs are shown in Table 2. Information for all 77 SNVs is shown in the Supplementary Information (Table S1). The quantile-quantile plot of the observed P values is depicted in Figure 1. The inflation factor of the genome-wide scan was 1.031, which indicates that the population structure was welladjusted. Figure 2 depicts scatter plots of P values obtained from genome-wide scan results for eGFR in T2D, which found that 77 close SNVs at 13q35.43-35.46 met significance (P<5 × 10^-8). In the 13q35.43-35.46 locus, stronger significances were observed upstream of the NBEA gene (Figure 3).

Figure 1: The quantile-quantile plot of the observed P values.

Figure 2: Depicts scatter plots of P values obtained from genomewide scan results for eGFR in T2D.

Figure 3: The 13q35.43-35.46 locus, stronger significances were observed upstream of the NBEA gene.

Confounding factor adjustment

The results of confounding factor adjustment analysis for the representative 8 SNVs in 519 T2D patients with urine protein data are shown in Table 3. Compared to the results in all 955 T2D patients, P values were more significant and β was larger in the age sex adjusted analysis in 519 T2D patients with urine protein data (model 1). With the addition of urine protein in model 1, P values were attenuated very slightly (model 2). With the addition of PCA components to model 1, P values were somewhat attenuated (model 3). With the addition of BMI, smoking, and alcohol consumption to model 1, P values were slightly attenuated (model 4) (Table 3).

Table 3: Confounding factor adjustment analyses in SNVs in chromosome 13 associated with eGFR in T2D with urine protein data.

Functional annotations for SNVs associated with eGFR in T2D in the 13q35.43-35.46 locus

We examined the expression quantitative trait loci (eQTL) relationship between the 77 significant SNVs and NBEA gene. From the GTEx database [25], 31 eQTL hits were found, of which expression levels were significantly associated with the 13q35.43-35.46 variants (P<1.0E-04; Table 4 and Supplementary Information (Table S1)).

Table 4: Functional annotations for SNVs associated with eGFR in T2D in the 13q35.43-35.46 locus.

Replication of previously reported SNVs in the Asian population

We performed a replication study on the 7 out of 24 newly identified loci associated with eGFR based on Scr that achieved directionconsistent significance in the trans-ethnic analyses of Asian participants in the meta-analysis conducted by Pattaro et al. [4]. Only one SNP (rs163160) at the intron of the KCNQ1 gene was statistically significant (P=0.032, EA FRQ=0.132), however, the direction was inconsistent. We could not replicate nine SNVs at the ELMO1 gene that were reported previously in a Japanese population [22]. However, 8 SNVs at the ELMO1 gene in the present study showed significant association with eGFR (P<0.025) (Table 5).

Table 5: Replication analysis using the J-MICC samples for previously-reported SNVs.

In the present GWAS of participants with T2D in a Japanese population, we found 77 novel SNVs located upstream of the NBEA gene that were associated with eGFR. The NBEA gene involves a member of a large, diverse group of A-kinase anchor proteins that designate the activity of protein kinase A to specific subcellular locations by binding to its type II regulatory subunits. NBEA is expressed in blood cells, the brain, internal systems such as the kidneys and colon, and secretory systems such as the pancreas and adrenal glands [26]. A wide variety of diseases are associated with NBEA including migraine in bipolar disorder [27], idiopathic autism [28], schizophrenia [29], major depression [30], substance abuse [30], and multiple myeloma [31,32]. Olszewski et al. found a significant association for two upstream SNVs in NBEA , rs17775456 and rs7990537 with BMI as a continuous quality and trends for weight among the overweight adult men. They also found that NBEA+/2 mice result in moderately elevated body weight during early adulthood. Increased insulin concentrations are consistent with this phenotype [33]. Despite these, there are no reports indicating an association of this gene with renal function in diabetes or in general. Our findings warrant replication studies and further functional study in addition to support the observed association between SNVs located upstream of the NBEA gene locus and eGFR in T2D.

Proteinuria is a marker of nephropathy in diabetes. Plasma proteins leak into the urine in increased glomerular permeability. These proteins taken up by proximal tubular cells can initiate an inflammatory process and interstitial scarring, ending up in fibrosis [34]. An important role of advanced glycation end products in the pathogenesis of proteinuria and degenerative changes in kidney was recently recognized [35]. In one of our confounding factor adjustment analyses, the addition of urine protein in the model did not attenuate the associations of newly found SNVs located upstream of the NBEA gene with eGFR. This indicated that the gene has an effect on eGFR independently of urinary protein.

The other confounding factor adjustment analyses, adjusting for PCA scores, or for BMI, smoking and alcohol intake did not result in an alteration of the association between the gene and eGFR. The P values were more significant and β was larger in the age sex adjusted analysis in 519 T2D patients with urine protein data compared to the results in all 955 T2D patients, despite a reduction in the number of participants. The reason for this finding is not clear. It is possible that data quality in T2D with urine protein data were better than in T2D without protein urine data.

The KCNQ1 gene was significant in our replication of prior reported SNVs in the Asian population. This was one of 7 replicated loci that were associated with eGFR in the trans-ethnic analyses of Asian participants in the meta-analysis study by Pat taro et al. The KCNQ1 gene involves a voltage-gated potassium channel that is required for the repolarization phase of the cardiac action potential. Disease conditions associated with KCNQ1 include Long QT Syndrome [36-41] and gestational diabetes mellitus [42]. Although we could not replicate the nine SNVs at the ELMO1 gene that were reported previously to be associated with nephropathy in diabetes in a Japanese population by Shimazaki et al. [22], we found 8 SNVs at the ELMO1 gene in the present study that showed significant association with eGFR (P<0.025). Success in replication studies of ELMO1 in non- Japanese populations was reported [43-47].

We should state the limitations of this study. First, we did not make a replication study in a different population. Second, there is not a large number of participants with T2D, because J-MICC is not a study specific for T2D. Third, we only had semi-quantitative urine protein data in a part of the patients, and we lacked important variables such as duration of T2D.

In conclusion, we have discovered that the 13q35.43-35.46 locus is associated with eGFR in T2D participants among a Japanese population. Future studies are needed to examine the biological mechanism that relates the locus and renal function in T2D.

We would like to thank all the staff of the Laboratory for Genotyping Development, Center for the Integrative Medical Sciences, RIKEN, and the staff of the BioBank Japan project.

The following grants supported fully or in part this study: Grantsin- Aid for Scientific Research for Priority Areas of Cancer (No. 17015018) and Innovative Areas (No. 221S0001), and JSPS KAKENHI Grants (No. 16H06277 and 15H02524) from the Japanese Ministry of Education, Culture, Sports, Science and Technology, a Grant-in-Aid from the Japan Society for the Promotion of Science (JSPS) KAKENHI Grant B Numbers 24390165, 20390184, 17390186; (Grant-in-Aid for Young Scientists) Numbers 15K19236 and 17K15840, funding for the BioBank Japan Project from the Japan Agency for Medical Research and development since April 2015, and the Ministry of Education, Culture, Sports, Science and Technology from April 2003 to March 2015.