Genome-wide landscape establishes novel association signals for metabolic traits in the Arab population

While the Arabian population has a high prevalence of metabolic disorders, it has not been included in global studies that identify genetic risk loci for metabolic traits. Determining the transferability of such largely Euro-centric established risk loci is essential to transfer the research tools/resources, and drug targets generated by global studies to a broad range of ethnic populations. Further, consideration of populations such as Arabs, that are characterized by consanguinity and a high level of inbreeding, can lead to identification of novel risk loci. We imputed published GWAS data from two Kuwaiti Arab cohorts (n = 1434 and 1298) to the 1000 Genomes Project haplotypes and performed meta-analysis for associations with 13 metabolic traits. We compared the observed association signals with those established for metabolic traits. Our study highlighted 70 variants from 9 different genes, some of which have established links to metabolic disorders. By relaxing the genome-wide significance threshold, we identified ‘novel’ risk variants from 11 genes for metabolic traits. Many novel risk variant association signals were observed at or borderline to genome-wide significance. Furthermore, 349 previously established variants from 187 genes were validated in our study. Pleiotropic effect of risk variants on multiple metabolic traits were observed. Fine-mapping illuminated rs7838666/CSMD1 rs1864163/CETP and rs112861901/[INTS10,LPL] as candidate causal variants influencing fasting plasma glucose and high-density lipoprotein levels. Computational functional analysis identified a variety of gene regulatory signals around several variants. This study enlarges the population ancestry diversity of available GWAS and elucidates new variants in an ethnic group burdened with metabolic disorders.


Introduction
The post-oil era of the Arabian Peninsula has witnessed a substantial increase in the prevalence of metabolic traitrelated disorders, such as obesity, dyslipidemia, hypertension, and type 2 diabetes mellitus (T2D), in its population. Despite the high prevalence of metabolic-related disorders in the Arabian Peninsula (Abuyassin and Laher 2015;Al Rasadi et al. 2016;Al Sifri et al. 2014;Channanath et al. 2013;Klautzer et al. 2014;Ng et al. 2014;Tailakh et al. 2014), there is a lack of convincingly identified genetic determinants for metabolic traits in people from this region. The global genome-wide association studies (GWAS) performed for metabolic diseases and traits overrepresent people of European ancestry (Mills and Rahal 2019). Only a few GWAS from the Arabian Peninsula are reported in the literature, including those on unrelated individuals from Saudi Arabia (Ram et al. 2017;Wakil et al. 2016) and Lebanon (Ghassibe-Sabbagh et al. 2014), on an extended family from the United Arab Emirates (Al Safar et al. 2013), and our own studies from Kuwait (Hebbar et al. 2017a(Hebbar et al. , b, 2018(Hebbar et al. , 2020. Although these studies identified few novel risk variants, they achieved virtually no success in replicating risk variants that had been discovered in non-Arab populations. While the inability to reproduce established genetic risk variants in these studies may possibly be attributed to different genetic architectures in terms of gene-environment interactions (Fahed et al. 2012) and a different pattern of metabolic disorders (such as dyslipidemia (Al Rasadi et al. 2016)) in Arab populations, it is more likely to be caused by technical factors, such as insufficient genome coverage provided by single nucleotide polymorphism (SNP) arrays that are used to genotype the discovery cohorts (Hebbar et al. 2019).
Imputation and meta-analysis have become standard practices in GWAS delineating genetic risk loci associated with complex disorders (Marchini et al. 2007;Pei et al. 2010). The imputation of genotypes for genetic variants that are untyped in the array increases the information provided by each microarray by accurately evaluating the evidence for association at genetic markers that are not directly genotyped (Li et al. 2009). Current publicly available imputation reference panels, many of which are based on haplotypes identified by the 1000 Genomes (1000G) Project (1000Genomes Project Consortium et al. 2015, accurately predict genotypes both for common variants (minor allele frequency [MAF] ≥ 5%) and low-frequency variants (0.5% ≤ MAF < 5%) across diverse populations (Mitt et al. 2017;Vergara et al. 2018). For example, Ghassibe-Sabbagh et al. (Ghassibe-Sabbagh et al. 2014) derived a marker set of > 5 million SNPs (either directly genotyped or imputed using the 1000G Project reference panels), showing an association with two established T2D loci (rs7766070/CDKAL1 and rs34872471/TCF7L2) at genomewide significance in a Lebanese population sample. Metaanalysis improves the capability to detect associations (Zeggini and Ioannidis 2009) by combining summary results from independent GWAS (usually with imputed genotypes) and has helped to markedly enlarge the catalog of GWA-identified risk loci from single studies for disorders such as T2D (Mahajan et al. 2018;Saxena et al. 2012).
Since the advent of GWAS, many risk loci for metabolic traits have been globally identified, concentrating mainly on the European population (Bustamante et al. 2011;Need and Goldstein 2009). Although there has been an increase in the number of GWAS from Asia (mostly East Asia) and Africa, the Eurocentricity remains prominent (Popejoy and Fullerton 2016). Studies have shown that (i) while a large number of Eurocentric risk loci associations replicate in one or the other non-European population, novel risk loci are often discovered in ethnic populations; and (ii) Eurocentric associations can have different effect sizes in ethnic populations (Carlson et al. 2013). There is a lack of studies that evaluate the transferability of established risk loci in Arab populations. Furthermore, Arab populations are characterized by a high level of inbreeding due to the practice of consanguineous marriages, often between first cousins. The consideration of such ethnic populations can lead to the identification of novel risk loci.
In the present study, we imputed two genome-wide genotype datasets from two independent cohorts of Arab ethnicity (from Kuwait) using the 1000G Phase 3 haplotype reference panel. We also performed statistical tests for associations with well imputed common variants (≥ 5% frequency) of 13 quantitative metabolic traits. We combined summary statistics from the imputed datasets using metaanalysis and compared the resulting associations with those reported in global GWAS (as listed in the GWAS Catalog, which compiles the results of published GWAS as a curated resource of SNP-trait associations (Welter et al. 2014)) or in other genetic studies published in the literature. Our study identified 70 risk variants from nine genes associated with metabolic traits at genome-wide significant p values. Many other risk variants were observed at p values borderline to genome-wide significance or with suggestive evidence of association. Of the identified variants, 349 were established variants reported in global GWAS or in the literature. The established markers seen replicating in our study cohort were characterized in comparison with those that are not seen in our study cohort for effect size and projection of the required sample size. Identified variants were characterized by computational functional analysis. We performed finemapping analysis to identify candidate causal variants. Table 1 presents the demographic and clinical characteristics of the study participants. The quality control (QC) procedures resulted in retaining 1434 samples and 118,793 SNPs using the MetaboChip and 1298 samples and 674,131 SNPs using the OmniExpress BeadChip. The two cohorts were genotyped in different times but at the same site. A total of 34,800 SNPs was genotyped in both the BeadChips. Only in the case of 15 SNPs, the proportional test p values for differences in allele frequencies were significant; however, the significance vanished when Bonferroni correction was applied. We performed principal component analysis by merging the SNPs with data from 1000G Project populations and performing linkage disequilibrium (LD) pruning to estimate the heterogeneity of each study cohort. The LDpruned dataset included 43,537 SNPs for the cohort genotyped using the MetaboChip and 66,645 SNPs for the cohort genotyped using the OmniExpress BeadChip. Scatterplots presenting the first three principal components derived from a merged dataset of the samples genotyped on each of the two BeadChips and representative populations from the 1000G Project are presented in Supplementary Figure S1 (A-F). These scatterplots depict similar positioning for the individuals from both of the study cohorts. Figure 1 presents the flow of data from the two genotyping platforms through the imputation, meta-analysis, and QC stages to the final choice and characterization of association signals.

Imputation and quality score
Our imputation strategy used 1000G Project Phase 3 data (1000Genomes Project Consortium et al. 2015. Utilizing the Michigan Imputation Server (MIS), we imputed 48,424,667 markers for the 1298 samples genotyped with the OmniExpress BeadChip and 46,597,973 markers for the 1434 samples genotyped with the MetaboChip (Supplementary Table S1). Supplementary Figure S2A presents a plot of the imputation mean quality score (Rsq) for the imputed markers at various MAF values. As expected, the Rsq increased as the MAF increased, and at every MAF value, the Rsq was higher with the OmniExpress BeadChip platform, which has a higher array SNP content than the MetaboChip platform.

Final dataset of imputed markers
Supplementary Figure S2B shows the proportion of imputed variants falling into different MAF ranges for each of the two platforms before and after filtering the variants at Rsq > 0.50. We observed that a large proportion of variants with an MAF < 5% had poor Rsq values, whereas those with an MAF ≥ 5%, constituting a smaller proportion of the total variants, were confidently imputed with an average Rsq > 0.75 on both platforms. Using a threshold of ≥ 0.5 for the imputation mean quality score, 12.2 million and 2.16 million SNPs were imputed for the OmniExpress BeadChip and MetaboChip, respectively. To create our final dataset, we filtered these for common markers (MAF ≥ 5%), resulting in an imputed OmniExpress dataset of 6.08 million markers and an imputed MetaboChip dataset of 1.52 million markers (mean Rsq = 0.822 and 0.759, respectively) (Supplementary Table S1).

Concordance of imputed genotypes
In an attempt to examine the quality of imputation using 1000G in our study, we validated the genotypes imputed in our study by way of comparing the genotype frequencies at the imputed markers with those genotyped in populations from the region as reported in either our in-house Kuwaiti Arab Exome Variant database  or in the Greater Middle East (GME) Variome Project ) (which presents GME genetic variations). Our inhouse exome variant database, derived from 291 samples, contained 5864 imputed markers (with a frequency ≥ 5% and Rsq ≥ 0.5), of which only five (0.09%) showed significantly different genotype frequencies (Fisher's exact test, p adjusted ≤ 0.05). All five variants were either multiallelic or indels. Similarly, the GME database contained 8550 imputed markers (with a frequency ≥ 5% and imputation quality of Rsq ≥ 0.5). Of them, only 62 imputed variants (0.73%) showed significantly different genotype frequencies (Fisher's exact test, p adjusted ≤ 0.05). Interestingly, these included 33 variants located in the major histocompatibility complex region, which is known to exhibit poor imputation quality (Matzaraki et al. 2017). The remaining 29 variants were either multiallelic or indels. It is also possible that the higher fraction of inconsistency observed with GME can be attributed to population differences within the Greater Middle Eastern region. These observations indicate that the quality of imputation using 1000G is good and is in agreement with genotyped data from Kuwait and the region. We removed all variants that showed significant differences in genotype distribution from further analyses.

Summary of associations observed between metabolic traits and imputed variants
Statistical analyses of the associations between the variants (genotyped and imputed) and the 13 metabolic traits identified 978 unique associations involving 821 unique variants way of comparing to association signals published in GWAS Catalog for the 313 search terms relating to metabolic traits from 251 gene loci (Supplementary Dataset S1). Variants from 237 of the 798 association signals reported in Supplementary Dataset S1 exhibited higher effect allele frequencies in Arabs as compared in Europeans and the difference was statistically significant as per proportional test for comparison; of these 237 signals, 172 (73%) were not seen replicated in global GWA studies. When the association signals were not at a level of genome-wide significance or at least borderline to genome-wide significance, they only involved the established variants (annotated in the GWAS Catalog for traits relating to metabolic processes). Of these 821 variants, (i) 70 variants (from nine genes) formed 72 associations at genome-wide significance (p value < 5.0 × 10 −08 ), (ii) 440 variants (from 76 genes) formed 455 associations of borderline to genomewide significance (p value < 1.0 × 10 −06 and > 5.0 × 10 −08 ), and (iii) 319 variants (from 181 genes) formed 451 established associations (reported in the GWAS Catalog for the 313 traits relating to the 13 study-specific metabolic traits) at p values of suggestive evidence of association (> 1 × 10 −06 and ≤ 0.05). LD pruning at r 2 = 0.10 suggested nine LD-independent variants (two variants of genome-wide, two variants of borderline, and five variants of suggestive associations) among 821 variants. Table 2 presents the overall summary statistics of the extent of associations observed between variants and the tested metabolic traits, along with comparisons with the established 7668 SNP associations from the GWAS Catalog (see "Methods" for the protocols used to search the GWAS Catalog with search terms relating to metabolic traits and disorders, as listed in Supplementary Table S2).

Associations observed at genome-wide significance in meta-analysis
As mentioned above, 70 unique variants (15 established risk variants and 55 novel risk variants) from nine gene loci (seven genomic regions) were found to be associated at a genome-wide significant p value. Of them, 63 were associated with high-density lipoprotein (HDL), one with triglycerides (TGs), one with low-density lipoprotein (LDL), three with systolic blood pressure (SBP), and two with each of fasting plasma glucose (FPG) and diastolic blood pressure (DBP). The associated gene loci are listed in the footnote to Table 2. SNP quality information of these top variants is presented in Supplementary Table S3. These 70 variants contributed to 72 associations with the traits (Supplementary Table S4). The quantile-quantile plots depicting the expected and observed − log10 (p values) of all variant associations for all the 13 traits for the two cohorts and for the meta-analysis along with the values for genomic inflation (λ) are presented in Supplementary Figure S3. Manhattan plots depicting the − log 10 (p values) from the GWAS for all the traits in our meta-analysis are presented in Supplementary Figure S4. Regional association plots for regions of 500 Kb centered at the risk variants associating with the metabolic traits at genome-wide significance are presented in Supplementary Figure S5.

Fine-mapping analysis
Fine-mapping analysis of the seven identified genomic regions revealed 95% credible causal variants for nine SNPtrait association signals, corresponding to eight lead SNPs. Table 3 presents a detailed list of 95% credible causal variants along with their association statistics. Among these eight lead SNPs giving credible sets, three (rs1864163-HDL, rs7838666-FPG, and rs2835788-SBP) showed a posterior inclusion probability (PIP) ≥ 0.5, another three (rs112861901 for HDL, rs66505542 for TG, and rs2920844 for SBP and DBP), showed PIP > 0.3 to < 0.5; and the remaining two lead SNPs (rs76018028-HDL and rs10635970-LDL) showed a very low PIP.

SNP heritability analysis
Regional SNP heritability analysis explained an estimated variance of 2.5% in HDL for the 1-Mb region comprising rs1864163, 1.7% in TGs for the region comprising rs66505542, 1.4% in HDL for the region comprising rs76018028, and 1.1% in SBP for the region comprising rs2835788. All the remaining top signals showed ≤ 1% trait variance ( Table 3).

Association of novel variants in our meta-analysis with metabolic traits
The meta-analysis identified 472 unique novel SNP variants (from 76 genes) in our study cohort that were associated with one or more of the tested 13 metabolic traits, either at genome-wide or borderline to genome-wide significance. While none of these 472 variants were associated with any of the 313 traits relating to metabolic processes in the GWAS Catalog, polymorphisms from 57/76 genes were associated with metabolic traits (see Supplementary Dataset S1). Literature evidence supported the involvement of eight additional genes in metabolic processes (see Supplementary Dataset S1). Thus, the evidence supported a role in metabolic processes for 65/76 genes, while 11 genes (mostly uncharacterized) showed no apparent link to metabolic traits.

Observed association signals and consistency in traits associated in GWAS Catalog
Trait consistency in association signals with the GWAS Catalog was demonstrated in 83 variants (from 42 genes), comprising the direct category (see "Methods"). Such variants were associated with 11 traits (namely, height, fasting Table 2 Summary statistics on observed associations between variants (genotyped and imputed) and the tested 13 metabolic traits. Data from comparing these associations with the 7668 established SNP associations from GWAS Catalog (as accessed in February 2019) is also included # Considered are only those variants that are seen common between the two imputed data sets and are with MAF > 0.05 and HWE > 1E−05 and imputation mean score ≥ 0.5

Discovered association signals and overlap with broad classes of metabolic traits in GWAS Catalog
Interestingly, many of the variants that were found associated with the 13 study-specific metabolic traits in our cohort exhibited associations with members of the broad metabolic trait classes in the GWAS Catalog (see Fig. 2). The complete dataset of all the variant associations, with grouping and trait classes, is provided in Supplementary Dataset S2. The extent of the overlap in trait associations among the different classes of metabolic traits were as follows: (i) Among the variants that were associated with anthropometric traits in our cohort, 9 variants (7 genes) were also associated with lipid traits, 11 variants (eight genes) with cardiometabolic phenotypes, 11 variants (10 genes) with blood pressure traits, and 7 variants (6 genes with glycemic traits in the GWAS Catalog. (ii) Among the variants associated with lipid traits in our cohort, 48 variants (32 genes) were also associated with anthropometric traits, 22 variants (18 genes) with blood pressure traits, 48 variants (25 genes) with cardiometabolic phenotypes, and 31 variants (26 genes) with glycemic traits in the GWAS Catalog. (iii) Among the variants associated with glycemic traits in our cohort, 11 variants (9 genes) were associated with anthropometric traits, 2 variants (2 genes) with blood pressure traits, four variants (4 genes) with lipid traits, and 5 variants (5 genes) with cardiometabolic phenotypes in the GWAS Catalog. (iv) Among the variants associated with blood pressure traits in our cohort, 17 variants (14 genes) were associated with anthropometric traits, 2 variants (2 genes) with lipid traits, 2 variants (2 genes) with glycemic traits, and 7 variants (5 genes) with cardiometabolic phenotypes. The associated variants and genes shared with the different classes of metabolic traits are illustrated as Venn diagrams in Supplementary Figure S6 (for variants) and Supplementary Figure S7 (for genes). The associations comprising the broad category may illustrate that rather than trait inconsistency, pleiotropy and a common genetic basis exists among different classes of metabolic traits.

Transferability of established association signals for metabolic traits to the Arab population
We observed 349 unique SNPs (from 187 genes) established in global populations as associated with metabolic traits replicated in the Arab population through direct, indirect, or broad relationships. Of these unique variants, 83 (from 42 genes) were replicated through a relationship of direct association, 164 (from 77 genes) through indirect association, and 203 (from 131 genes) through broad association. The counts of transferable association signals at such variants from each global population to the Arab population are presented in Fig. 3, and similar counts at the gene level are presented in Supplementary Figure S8. The presented association signals were generally well illustrated in various studies (see Fig. 2) and in multiple global populations. Table 4 shows populations illustrating the replicability of associations from the direct relationship class. For example, the replicability of associations from the direct class for the traits of BMI, TG, and height has been observed in many populations, with BMI associations in Northern and Western European (CEU), East Asian (EAS), South Asian (SAS), Hispanic, Filipino, and Seychellois populations; TG associations in CEU, EAS, and Han Chinese in Beijing, China (CHB), populations; and height associations in CEU and EAS populations. Figure 4 presents the trend in effect sizes against increasing MAF values of the established variants that were seen replicating in our study cohort (Fig. 4a) and those that were not replicating in our study cohort (Fig. 4b). Study variants, replicating the established variants, showed a high effect size. On the other hand, non-replicating variants showed a low effect size.

Estimation of sample size for replicating and non-replicating variants
The estimated power and sample size (see "Methods" for power calculation) for the variants of different effect sizes and MAFs for different traits are presented in Supplementary Figure S9 (the left panel presents the results of established variants that were replicated in our study, and the right panel presents the results of established variants that were not replicated in our study). The figure illustrates that most of the established variants that were seen replicating in our study would attain genome-wide significant p values (5.0 × 10 −8 ) for most of the traits at 80% and even more power with a sample size of 10,000, while non-replicated markers would be less likely to replicate at genome-wide significant p values, even at a sample size of 20,000.

List of gene loci implicated in metabolic processes by our meta-analysis
This study established the association of 66 genes with anthropometric traits, 42 genes with blood pressure traits, 25 genes with glycemic traits, and 133 genes with lipid traits (Fig. 5) at different p value thresholds for the SNP-trait 1 3 associations (genome-wide significance of < 5.0 × 10 −8 , borderline to genome-wide significance of < 1.0 × 10 −6 and > 5.0 × 10 −8 , and suggestive p values of > 1.0 × 10 −6 ). An association at a suggestive p value was only considered if the variant was also listed as associated with metabolic traits in the GWAS Catalog. There were 24 genes associated with multiple anthropometric traits, five with both SBP and DBP, one with both FPG and HbA1c, and 28 with multiple lipid traits. Many of these genes (9/66 associated with anthropometric traits, 25/42 associated with blood pressure traits, 20/25 associated with glycemic traits, and 101/133 associated with lipid traits) were also found associated with metabolic traits in the GWAS Catalog.

Functional consequences of discovered variants
The distribution of the variants as per their functional consequences on gene structures as gathered (see "Methods") is presented in Supplementary Figure S10A. Up to 52% of the variants were intronic, and up to 30% were intergenic. Up to 6% were located upstream of a gene. Approximately 3% of the variants were seen in regulatory regions, and approximately 1.5% of variants led to amino acid changes in the encoded proteins. The distribution of the variants as per their proximity to the transcription start sites (TSSs) is presented in Supplementary Figure S10B. The majority of the 821 study variants were situated close to the TSS, with 56 variants within a 1-Kb region from the TSS, and 107 variants within a 2-Kb region from the TSS.

Functional prioritization of discovered variants
Analysis for functional prioritization of the variants led to valuable knowledge regarding the potential functional signatures around the study variants, and such signatures plausibly drive the course of downstream gene expression. As many as 752/821 variants could be ranked based on the functional niches. These 752 variants, along with their functional niches, are presented in Supplementary Dataset S3. Our study identified 33 of the top 100 ranked variants as novel risk variants. The top-ranked variants included the following: (i) the rs12740374 variant of the CELSR2 gene (associating with non-HDL, LDL, and TC at p values borderline to genome-wide significance), which carried a high score for DNase hypersensitive sites (HSs), DNase footprint, and transcription factor binding sites (TFBS) and was ranked 1; (ii) the variant rs3749147 of GPN1,ZNF512 (associating with TG at a suggestive p value), which overlapped with the promoter, located very close to the TSS, carried a high score for DNase HSs and TFBS, and was ranked 2; (iii) the variants rs7670 and rs4705745 of the DCP2 gene (associating with WC at p values borderline to genome-wide significance), which carried high scores for DNase HSs and were ranked 3 and 62, respectively; and (iv) the variants rs62355943 and rs72758038 of the MAP3K1 gene (associating with TG at p values borderline to genome-wide significance), which overlapped with the promoter and CpG shore, carried high scores for DNase HSs and TFBS, and were ranked 5 and 6, respectively.

The discovered risk variants and their ability to regulate expression of genes
A functional assessment of the variants (and the genes harboring the variants) in terms of their ability to up-and/or downregulate genes (as deduced by examining eQTL's from GTEx) is illustrated in Supplementary Figure S11 (A Figure S12).

Discussion
Arabian Peninsula is at the nexus of Africa, Europe, and Asia and has been implied in early human migration route out of Africa and in early inter-continental trade routes. The post-oil rich era has seen a drastic shift from nomadic style of living to more sedentary lifestyle and rapid nutrition transition resulting in a dramatic increase in the prevalence of metabolic-related disorders. However, the global genetic studies for metabolic studies were mostly performed on European populations followed by African American, East Asian, and South Asian populations. Populations from Arabian Peninsula have not been included in these global studies and there is limited genetic association data for metabolic traits in Arab population. Close-kin marriage and large families are cultural factors in the region. The practice of consanguineous marriages and living in isolation by community expectedly leads to increased endogamy, homozygosity, and accumulation of deleterious recessive alleles in the gene pool. The Arab population has been vulnerable to a plague of recessive genetic disorders. Consanguinity and inbreeding can also play an important role in the etiology of complex disorders (Rudan et al. 2006). Given such a unique genetic profile, studying such an Arab population is expected to augment  international efforts to identify genetic regulation of metabolic traits. Although published GWAS have reported several thousand genotype associations with metabolic traits, they have mostly focused on people of European ancestry and, to an extent, East Asian ancestry, and it has remained unclear whether these findings could generalize to every other population. The frequencies of risk alleles associated with a range of traits in GWAS can differ substantially between continental populations (Adeyemo and Rotimi 2010). Thus, it is crucial to assess how well these associations can be extended to populations with different continental ancestry. An increasing number of studies aimed at generalizing these associations to various ethnic populations are being pursued globally Langlois et al. 2016;Liu et al. 2012;Lu and Loos 2013). Fine-mapping of the replicating risk gene loci is performed to identify the true causal variants Mahajan et al. 2018;Wu et al. 2013). Determining the transferability of established risk loci and variants to different ethnic populations is essential for generalizing and reconfirming the findings from previous global collaborative studies as well as for enabling the successful transfer of the knowledge, research tools, resources, and drug targets generated by the global studies to a broad range of populations with different ethnicities.
This study performed genome-wide imputation of genotypes for variants that were untyped in the bead chips used in our previous GWAS, namely, the high-density OmniExpress BeadChip array (Hebbar et al. 2017a and lowdensity Cardio-MetaboChip array (Hebbar et al. 2017b). An assessment of the imputation quality obtained in this study (based on the imputation quality score and the proportion of variants above the threshold for the quality score at every bin of MAF, incremented by 5%), suggested that the imputation quality was unprecedentedly good with common variants compared to low-frequency variants (MAF < 5%). Furthermore, it was observed that even using the low-density SNP array (Cardio-Metabo BeadChip), a mean Rsq of 0.75 could be obtained for common variants. Although the 1000G ALL panel was used as the reference panel for imputation, a high concordance was observed between the genotype frequencies of the imputed markers and their frequencies as obtained from the GME Variome Project  or from our in-house Arab Exome Variant database . This suggests that the imputation of both highdensity and low-density SNP data using 1000G haplotypes as the reference panel successfully imputes common variants in the Arab population.
Our study identified 978 association signals (involving 821 variants from 251 genes) for the selected metabolic traits. Up to 95% of the 251 identified gene loci are supported by evidence from either annotation in the GWAS Catalog or from experimental studies reported in the literature (see Supplementary Dataset 1). We found 70 variants from nine gene loci associated with metabolic traits at a genomewide significance; these gene loci include BUD13, CETP,CSMD1,DYRK1A,HERPUD1,INTS10,LPL,and RTN4, which are known to be involved in metabolic processes. Of them, 17/70 variants were established risk variants for metabolic traits and are listed in the GWAS Catalog. The genes harboring the remaining variants were associated with metabolic processes as seen in the GWAS Catalog or from the literature.
Our study inferred three of the observed risk variants (rs112861901/[INTS10, LPL] associated with HDL, rs1864163/CETP associated with HDL, and rs7838666/CSMD1 associated with FPG) as causally implicated. Furthermore, the regional SNP heritability analysis demonstrated that the regions corresponding to the top signals showed low values, typically < 2%, except for the cholesteryl ester transfer protein (CETP) region (2.5%). These low values indicate that the proportion of metabolic trait variance that is attributable to these identified genetic variations is low among the study population. Although studies reporting a higher explained variance due to genetic variation in CETP exist [for example, Blauw et al. (2018) reported a high value of 16% for the explained variance in CETP concentration due to a set of 3 CETP variants], we cannot rule out the possibility of environmental interactions being responsible for larger relative contributions (Visscher et al. 2008) to the variance of traits in our study cohort.
The GWAS era has successfully associated thousands of genetic variants with risk for complex disorders and traits; for example, examination of the GWAS Catalog (accessed Feb 2019) found 7668 associated genetic variants for the 13 metabolic traits considered in this study. However, it still remains a challenge to translate these statistical associations to knowledge on how they functionally manifest to alter the biology underlying the disease risk (Edwards et al. 2013). A large part of this difficulty results from the fact that > 90% of disease-associated variants are located in non-protein coding regions of the genome (such as introns and intergenic regions) and not within the 1.5% coding part of the human genome (Hindorff et al. 2009;Maurano et al. 2012;Schaub et al. 2012). In line with these reports, we observed that up to 88% of the 821 study variants were located in introns or intergenic regions or upstream of genes as opposed to a mere 1.5% of the variants leading to amino acid changes in the encoded proteins. It has been proposed that many such variants located in non-protein coding regions might influence disease risk by way of altering the regulation of the target genes and they might form part of regulatory motifs (Gallagher and Chen-Plotkin 2018;Maurano et al. 2012;Schaub et al. 2012). The analysis of regulatory elements that we present in this study demonstrates that only around 3% of the 821 study variants were seen in regulatory motifs and that at the most 13% of the 821 study variants were seen within a 2-Kb region of TSSs. Variant prioritization based on the observed functional niches could uncover in our study a variety of regulatory signals around novel variants from genes such as DCP2, MAP3K1,HHIP,KANK2,DAGLA,MAFF,CETP,LOC646576,and APOA5. It is now established that GWAS hits for common diseases are enriched for eQTLs and that disease-associated risk variants are likely to regulate expression levels of target genes than would be expected by chance (Albert and Kruglyak 2015). In line with these reports, we found that up to 510 (62%) of the 821 study variants were annotated in GTeX as regulating 464 genes in 49 tissues; further use of the criteria of q value ≤ 0.05 (a measure of significance in terms of false discovery rate) to filter the eQTL variants, reduced the number of such variants to 62 (7.6% of the 821 study variants) which differentially regulated expression of 39 genes in 38 tissues. Deciphering the functional role of the GWAS-associated risk variants is a challenge and requires comprehensive computational and experimental functional analysis studies.
It is well recognized that several genes have pleiotropic effects on multiple disorders, such as among T2D, obesity, and dyslipidemia ); among obesity, cardiovascular disease outcomes, and cardiovascular risk factors (Rankinen et al. 2015); and among lipid metabolism and metabolic syndrome (Park et al. 2011). Our approach of using 7668 established variant associations for metabolic traits to classify the associations identified in our study aided in the identification of genetic loci with shared impacts on different metabolic processes. Some examples of such exemplary variants and genes include the following: (i) rs780093 of GCKR associated with TG shares associations with all five broad classes of metabolic traits (anthropometric, blood pressure, cardiometabolic, glycemic, and lipid traits); (ii) rs1260326 (GCKR), rs780094 (GCKR), and rs2925979 (CMIP) share associations with four broad classes of metabolic traits (anthropometric, cardiometabolic, glycemic, and lipid traits); and (iii) ten variants share associations with three broad classes of metabolic traits, including (a) rs6795735 (ADAMTS9-AS2) and rs6857 (PVRL2) with anthropometric, glycemic, and lipid classes; (b) rs15285 (LPL) with blood pressure, cardiometabolic, and lipid classes; (c) rs11977526 (LOC102723446) and rs198846 (HIST1H1T) with blood pressure, glycemic, and lipid classes; (d) rs2820315 (LMOD1) and rs564398 (CDKN2B-AS1) with anthropometric, cardiometabolic, and glycemic classes; (e) rs645040 (RPL31P23-PCCB) with anthropometric, cardiometabolic, and lipid classes; (f) rs11556924 (ZC3HC1) with anthropometric, blood pressure, and cardiometabolic classes; and (g) rs7248104 (INSR) with anthropometric, blood pressure, and lipid classes.
The replicability of GWAS-identified signals for consistency in allele frequencies and effect sizes has been assessed in different major continental populations (Haiman et al. 2012;Kraft et al. 2009;Li and Keating 2014;Marigorta and Navarro 2013;Marigorta et al. 2018;Ntzani et al. 2012). We, for the first time, evaluated the impact of effect sizes on the replicability of established risk variants in ethnic Arab populations. Our results confirm the common notion that large effect associations usually replicate well. Our estimation of the sample size for an Arab cohort required to replicate established associations at genome-wide significance revealed that even with a sample size of 10,000, the established variants that we found replicating in our study (at different p values) would attain genome-wide significance. Furthermore, even with a sample size of 20,000, the established variants that we found as non-replicating in our study were less likely to replicate at genome-wide significance.
Our previous publications using one or the other of the two study cohorts could identify risk variants at genomewide significance only under genetic models based on recessive mode of inheritance. Examples include: the studies using the OmniExpress cohort identified a recessive marker from RPS6KA1 associated with FPG (Hebbar et al. 2020) and 6 recessive markers from RPS6KA1, LAD1, Or5v1, CTTNBP2-LSM8, PGAP3 and RP11-191L9.4-CERK associated with TG ; and the study using the CardioMetabo cohort identified a recessive marker from ZNF106 associated with HbA1c (Hebbar et al. 2017b). The only exception to identifying recessive markers is the study wherein we identified an additive marker from TCN2 associated with WC using the OmniExpress cohort (Hebbar et al. 2017a). Adopting the approach used in the current study to perform meta-analysis on the summary statistics from both the cohorts led to identification of several additive markers as opposed to recessive markers at genome-wide significance-such as those from LPL (LDL), CETP (HDL), BUD13 (TG), CSMD1 (FPG), RTN4 (DNP, SBP), DYRK1A (SBP). It is further the case that some of the additive markers that we identified at suggestive p values in our previous studies (Hebbar et al. 2017b) appeared in the current study at genome-wide significance as in the case of markers from CETP (HDL) and BUD13 (TG). The ability to detect additive risk variants enabled us to explore the transferability of established markers (which are often identified through additive genetic model) to Arab population.
Ours study has certain limitations. (i) The association analysis considered only common variants (MAF ≥ 5%), and we set this criterion upon considering the small sizes of the datasets and the rapid degradation in imputation accuracy when low-frequency variants are considered. Thus, the transferability of low-frequency variants could not be examined in this study. However, we are considering "ROHregion-based multi-marker association tests in an attempt to improve analysis power" as follow-up to the current work. (ii) Because the two cohorts used in our study included both diabetic patients and healthy participants, the identified associations retained significance when the models were adjusted for the covariate of diabetes status. It is often the case that quantitative trait associations are performed on entirely non-diabetic participants or on entirely diabetic patients. In our earlier study (Hebbar et al. 2020), using one of the two cohorts mentioned here, we showed that 3/4 markers identified from the full cohort performed better in terms of retaining significance in the sub-cohort of entirely diabetic patients, and the fourth marker performed better in terms of retaining significance in the sub-cohort of participants free of diabetes. (iii) Although the Haplotype Reference Consortium (HRC) panel is more recent, we used the 1000G Project panel for imputation. There is no gain in HRC imputation in our study cohort as the HRC Release 1 (https ://www.haplo type-refer ence-conso rtium .org/parti cipat ing-cohor ts) used by the MIS server does not contain Middle Eastern samples. The 1000G panel has been used more often in the literature, and it includes comparatively more diverse parental populations. However, the choice between HRC and 1000G panels does not matter, since our study considers only common variants. Furthermore, it has been demonstrated that in individuals of African ancestry, the 1000G panel also showed higher performance compared with the HRC panel in terms of the number of imputed variants with high accuracy (Vergara et al. 2018), and in Estonian individuals, the accuracy and sensitivity of imputation for common variants seemed to be similar between the 1000G and HRC panels (Mitt et al. 2017). Arab-specific haplotype reference panels are still under development in the region. (iv) The Cardio-Metabo BeadChip used to genotype one of the two cohorts is a fine-mapping array with much fewer markers as opposed to the OmniExpress BeadChip used to genotype the second cohort. As a result, there are much fewer genotyped and imputed SNPs resulting for the cohort using the Cardio-Metabo chip (at 1.52 million as opposed to 6.08 million obtained for the other cohort with OmniExpress BeadChip) leading to spotty coverage of the genome. (v) The study lacks formal pathway and gene-set analyses-it is our intention to perform these analyses as part of our future studies to translate the identified statistical associations to functional manifestations. (vi) While the values for genomic inflation factor (λ) were consistently close to 1.0 in the GWAS summary statistics from each of the two cohorts for each of the 13 traits, summary statistics from the meta-analysis exhibited slight inflation to the order of λ = 1.09 in the cases of HDL, SBP and Height. This can be a concern because of the moderate sample size of 2732. Considering that we have taken care of cryptic relatedness among samples in each of the two cohorts as well as between the cohorts, it is possible to attribute that these traits are probably associated with elevated polygenicity. (vii) The sample size was small. However, our study has succeeded in achieving its aim of finding causal variants and identifying potential novel hits specific to the Arab population.
In conclusion, this is the first presented study using an Arab population to perform genome-wide imputation and meta-analysis. We identified novel variants from seven established genes and two novel gene loci at genome-wide significance in the Kuwaiti population. We further identified 83 established SNP-metabolic trait association signals from the GWAS Catalog in their exact forms (the same variant associated with the same trait as in the GWAS Catalog) in the Kuwaiti population. We observed that of the identified association signals (for the four classes of anthropometric, blood pressure, lipid, and glycemic traits), a set of 330 SNPs (from 182 genes) were annotated in the GWAS Catalog as associated with other traits from the same class or from different classes of metabolic traits. These gene loci are indicative of a common genetic basis across the different metabolic processes. The replication of such established association signals provides further external validation and demonstrates the transethnic translatability of the risk loci associated with metabolic traits. The report conveys the following overall messages: the study (i) examined an understudied population characterized by inbreeding; (ii) demonstrated imputation from low density arrays to 1000 genomes haplotype reference panel; (iii) identified significant risk loci for metabolic traits, albeit previously known, in this novel population of Arabs; and (iv) indicated important metabolic pathways and genes.

Study participants
The study participants were selected from two independent cohorts recruited in our previous studies as part of two approved research projects at Dasman Diabetes Institute, Kuwait: (i) a cohort of 1965 samples resulting after QC steps based on genotype data (Hebbar et al. 2017b) and (ii) a cohort of 1913 samples resulting after QC steps . The participants were randomly recruited under protocols that were approved by the scientific and ethics advisory boards at Dasman Diabetes Institute. The study consists of two participant groups. The first group included a random representative sample of adults (> 18 years of age) of Arab ethnicity across the six governorates of the State of Kuwait. A stratified random sampling technique was used to select native Kuwaiti participants from the computerized register of the Public Authority of Civil Information, a government body that keeps and maintains personal information records of both Kuwaiti citizens and expatriates (including citizens of other Arab countries from the region). The second participant group comprised patients with diabetes or prediabetes seeking tertiary medical care in clinics at the Dasman Diabetes Institute, visitors to our nutrition programs and fitness center, visitors to our open day events (where various diagnostic screening services are offered), and visitors to our campaigns at primary health centers and blood banks in each of the six governorates of the State of Kuwait. As previously mentioned in our studies (Hebbar et al. 2017b, upon confirming that the participants had fasted overnight, they signed consent forms, blood samples were collected, vital signs were measured, and the nationality and ethnicity of each participant were confirmed via a rigorous questionnaire that addressed parental lineages up to three generations. Details of medications taken by the participants for lowering lipid levels, diabetes, and hypertension were collected and used in correction procedures with the association statistics. Data on illnesses (e.g., diabetes and cardiovascular complications) were also recorded. The guidelines of the Institutional Ethical Review Committee were followed for the proper collection of blood samples and the measurement of vital signs. All clinical assays for measurements of trait outcomes for both the cohorts were performed at a CAP (College of American Pathologists) accredited laboratory at Dasman Diabetes Institute.

Genome-wide genotyping
We previously genotyped 1965 individuals of Arab ethnicity using the Illumina Human Cardio-Metabo BeadChip (MetaboChip) (Hebbar et al. 2017b) and 1913 individuals of Arab ethnicity using the Illumina Human OmniExpress BeadChip (OmniExpress) (Hebbar et al. 2017a. The protocols that were followed to perform genome-wide genotyping and QC are also described in these studies. We used PLINK 1.9 tools (Purcell et al. 2007) for data management and QC. Genotypes were converted to the National Center for Biotechnology Information's Genome Reference Consortium Human Build 37 (hg19) to ensure consistent SNP phasing for each genotyping array. The genotyping QC procedures included the following: (i) samples with a call rate > 95% were retained; (ii) samples were checked against relatedness (both within each of the two cohorts and between the cohorts) and ancestry mismatch: one randomly selected representative of each set of related samples was retained, and samples with ancestry mismatch were removed; (iii) samples with heterozygosity > median + 3* interquartile range were excluded; and (iv) SNPs with a call rate > 98%, Hardy-Weinberg equilibrium (HWE) > 10 −6 , and MAF > 1% were retained. In a further QC step, samples that were common to both the study cohorts were removed from one of the sample sets.

Phenotyping study participants for metabolic traits
For each participant, measurements were available for 13 quantitative metabolic traits from four classes: (i) anthropometry: height, weight, BMI, and WC; (ii) blood pressure: SBP and DBP; (iii) glycemia: FPG and HbA1c; and (iv) serum lipids: LDL, HDL, TC, TG, and non-HDL cholesterol (obtained by subtracting the HDL from the TC level).

Imputation
Genotype imputation was performed on the MIS (Das et al. 2016) using the 1000G ALL reference panel from 1000G Project Phase 3 V5 for imputation, Eagle v2.3 for phasing, and Minimac3 as the algorithm for imputation . Prior to data submission to the MIS, the input data was checked against the 1000G reference SNP list using the HRC or 1000G Imputation Preparation and Checking Tool (https ://www.well.ox.ac.uk/~wrayn er/tools /HRC-1000G -check -bim.v4.2.5.zip). Strand designations were corrected to the forward strand, and reference/alternate (REF/ ALT) allele designations were corrected using PLINK2 and the design files for OmniExpress and MetaboChip. Variants from chromosome X were excluded from imputation because the MetaboChip had very few SNPs from this chromosome. The imputation quality was quantified using the Rsq score (which estimates the squared correlation between imputed and true genotypes). Imputation was separately performed for genotypes from OmniExpress (1298 samples) and MetaboChip (1434 samples).

Validating imputed genotypes
To evaluate the quality of imputation we evaluated whether genotype distribution at imputed variants in our study are consistent with distribution at the same variants genotyped in other studies. We examined consistency in genotype distribution at imputed variants versus genotype distribution at the same variants in our in-house Kuwaiti Arab Exome Variant database , and in genotype distribution at imputed variants versus genotype distribution at the same variants in GME Variome dataset . We used Fisher's exact test to check the genotype distribution consistency for the above comparisons, and then we adjusted the derived p values using the Benjamini-Hochberg procedure. Adjusted p values < 0.05 were considered as genotype proportion inconsistency.

Sensitivity analysis and transformations performed on the trait measurements
There was a concern that the glycemic, lipid, and blood pressure trait values that were determined in individuals who were receiving medication would not represent the naturally observed values in the population. We addressed this by applying adjustments to the measurements of the phenotypes by adding an average effect size to the phenotypes in the subjects under drug treatment. For FPG and HbA1c, we examined the in-house Knowledge-Based Health Records, an electronic health record system maintained at our institute, to calculate the average effect sizes among the native Kuwaiti T2D patients taking glucose-lowering medication who had improved their glucose levels (average reduction, − 2.59 mmol/L; p < 0.001). For lipid corrections, we used the recommendations of the Global Lipids Genetics Consortium and the Genetic Investigation of Anthropometric Traits Consortium (Willer et al. 2013). For antihypertensive medication, we used data reported in the literature (Tobin et al. 2005). The corrections implemented were as follows: (i) for participants taking lipid-lowering medication, we preadjusted TC to TC/0.8 and LDL to LDL/0.7, and non-HDL was determined by subtracting HDL from the preadjusted TC; (b) for participants taking antihypertensive drugs, we preadjusted SBP and DBP by adding 15 mmHg and 10 mmHg, respectively; and (c) for participants taking glucose-lowering medication, we preadjusted FPG and HbA1c values by adding 2.59 mmol/dl and 1%, respectively. Before generating the residuals, all the traits were adjusted (regular correction) for age, age 2 , and sex and for population stratification using the first four principal components (derived for each of the two study cohorts). Quantitative trait association tests require that the residual is normally distributed. Hence, we performed an inverse normal transformation on the raw residuals for the traits (except for TG, FPG, and HbA1c, for which log inverse transformation was performed, as it yielded a better Gaussian distribution). Association tests were performed with these inverse trait distributions.

Association tests and meta-analysis
Genotype associations for the 13 quantitative metabolic traits were separately tested for each array platform using all genotyped and imputed SNPs that passed the QC threshold metrics (Rsq > 0.05 and MAF ≥ 5%). RVTESTS software (Zhan et al. 2016) was used to perform association tests using linear regression analysis for an additive genetic model. The imputed markers that were passed to the subsequent meta-analysis stage were required to meet the following QC criteria: (i) consistency in the allele and genotype frequencies of the markers between the cohorts genotyped with the two BeadChips; (ii) HWE < 10 −6 ; (iii) p value for association < 0.05 in each of the two datasets; and (iv) consistency in the direction of association between the two datasets. Meta-analysis was performed using METAL software (Willer et al. 2010).

p value thresholds for associations
The p value threshold was set at 5.0 × 10 −8 for genome-wide significance, < 1.0 × 10 −6 and > 5.0 × 10 −8 for borderline to genome-wide significance, and > 1.0 × 10 −6 and < 0.05 for suggestive evidence of association. Associations at suggestive p values were only considered if they were listed in the GWAS Catalog.

Examining the NHGRI-EBI GWAS Catalog for reported association signals and variants
Reported traits for metabolic processes in the GWAS Catalog are highly diverse. The search terms relating to the 13 metabolic traits tested in our study belonged to five broad classes of curated trait search terms relating to anthropometry and obesity, blood pressure and hypertension, glycemia and diabetes, lipid profiles, and cardiometabolic phenotypes. We extracted 313 search terms (Supplementary Table S2) relating to the 13 study-specific metabolic traits from the NHGRI-EBI GWAS Catalog v1.0 (Buniello et al. 2019) (as accessed in February 2019): 99 search terms for the trait class of anthropometry and obesity, 46 for blood pressure and hypertension, 75 for glycemia and diabetes, 55 for lipid profiles, and 38 for cardiometabolic phenotypes. The GWAS Catalog contained 7668 SNP association signals for these terms. We compared the results of our meta-analysis to the known association signals from the GWAS Catalog and further removed any association that showed an inconsistent direction of effect from our subsequent analysis.

Classifying variants associated with the tested 13 metabolic traits
Every variant that was associated with one or the other of the 13 metabolic traits (from the classes of anthropometry, glycemia, lipid, and blood pressure) in our meta-analysis was classified as either established or novel, depending on whether they appeared in the list of the 7668 known SNP association signals.

Classifying the meta-analysis associations involving established variants
We classified meta-analysis associations involving established variants as direct, indirect, or broad, by comparing the associated metabolic traits in our study with known SNP association signals from the GWAS Catalog (as accessed in February 2019) as follows: (i) direct: the same trait was associated with the variant in the study cohort and the GWAS Catalog (e.g., a variant associated with BMI in our study cohort that was associated with BMI in the GWAS Catalog); (ii) indirect: the associated trait in the GWAS Catalog was a member of the class of traits related to the associated trait in our meta-analysis (see Supplementary Table S2 for classes of search terms relating to the groups of the 13 tested traits) (e.g., a variant associated with BMI in the study cohort that was associated with any trait described by the 99 search terms listed under the anthropometry and obesity class); or (iii) broad: the trait associated with the variant in our meta-analysis was related to any of the 313 metabolic traits (listed under the classes of anthropometry and obesity, blood pressure and hypertension, glycemia and diabetes, lipids, or cardiometabolic phenotypes) in the GWAS Catalog (e.g., a variant associated with BMI in the study cohort that was associated with adiponectin levels in the GWAS Catalog).

Prioritizing variants based on a fine-mapping approach
In this study, we explored the identified top-associating signals (p value ≤ 5 × 10 −08 ) using FINEMAP software to identify the plausible causal variants and regional heritability (Benner et al. 2016). FINEMAP efficiently explores sets of the most plausible causal variant configurations from a given genomic region using a shotgun stochastic search algorithm based on summary statistics from the meta-analysis and pairwise LD correlation statistics between the variants. We included variants within a ± 500 Kb region around the lead SNP with inverse-variance weighted summary statistics, computed the LD correlation matrix for both datasets independently using LDstore software (Benner et al. 2017), and then averaged them according to the following formula: meta_LD_matrix = (n1 × R1 + n2 × R2)/(n1 + n2), where R1 and n1 are the LD matrix and sample size of imputed data from OmniExpress, respectively, and R2 and n2 are the LD matrix and sample size of imputed data from Metabo-Chip, respectively. We created 95% credible causal variant sets by ranking the variants from the region based on their decreasing posterior probability of association and estimated regional heritability to understand the proportion of phenotypic variance of the complex trait attributable to a credible set of SNPs.

Estimation of power for replicating and not replicating variants at genome-wide significance
We performed power calculation for replicating and not-replicating variants to estimate power required to expect them at p value ≤ 5E−08 using R scripts (available at https ://githu b.com/kaust ubhad /gwas-power ) that were developed using formulas as per Visscher et al. (2008). The power_beta_maf function from the script was used to calculate power based on a range of values of beta and MAFs (while values for sample size (n) and pval were fixed).

Functional variant prioritization
We used the Ensembl Variant Effect Predictor tool (McLaren et al. 2016) to annotate the functional consequences of variants. We assessed the proximity of the variants to the TSS using ChIPpeakAnno (Zhu et al. 2010) and TSS.human. GRCh38 data (available from the Bioconductor software suite at https ://bioco nduct or.org/packa ges/). We used the tools of SuRFR (Ryan et al. 2014) (SNP Ranking by Function R package), which interacts with SAILR (SNP Annotation Information List R package) to prioritize variants based on the likelihood of their function using features such as promoter, CpG, CpG shore, DNase HSs, DNase footprints, TFBS, TSS, chromatin states (histone acetylation and methylation), conserved sequences, and enhancers. To aid in the variant prioritization, we used background variants from European populations and used a pre-trained weighting model for complex disease variants, with default parameter values (Ryan et al. 2014).

Expression quantitative trait loci analysis
We examined genotype-tissue expression data using GTeX v8 (https ://www.gtexp ortal .org) to assess the involvement of the discovered variants in the regulation of gene expression. Two levels of analysis were carried out: (i) all eQTL's with p value ≤ 0.05 were considered; and (ii) only those eQTL's with p value ≤ 0.05 AND q value ≤ 0.05 were considered. The q value is similar to the well known p value, except it is a measure of significance in terms of the false discovery rate rather than the false positive rate (Storey and Tibshirani 2003). Beta distribution-adjusted empirical p values from FastQTL were used to calculate q-values.

3
Funding No external funding was obtained for this study.
Availability of data and material All the data on the variants and their associations with metabolic traits are provided as supplementary data sets. GWAS summary statistics is available to researchers upon writing to the corresponding authors.

Compliance with ethical standards
Conflict of interest The authors affirm that there is no conflict of interest to declare.

Ethics approval
The study was approved by the Ethical Advisory and Review Committee at Dasman Diabetes Institute, Kuwait.
Informed consent Written informed consent forms were signed by the participants.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creat iveco mmons .org/licen ses/by/4.0/.