An investigation of obesity susceptibility genes in Northern Han Chinese by targeted resequencing.

Our earlier genome-wide linkage study of body mass index (BMI) showed strong signals from 7q36.3 and 8q21.13. This case-control study set to investigate 2 genomic regions which may harbor variants contributed to development of obesity.We employed targeted resequencing technology to detect single nucleotide polymorphisms (SNPs) in 7q36.3 and 8q21.13 from 16 individuals with obesity. These were compared with 504 East Asians in the 1000 Genomes Project as a reference panel. Linkage disequilibrium (LD) block analysis was performed for the significant SNPs located near the same gene. Genes involved in statistically significant loci were then subject to gene set enrichment analysis (GSEA).The 16 individuals aged between 30 and 60 years with BMI = 33.25 ± 2.22 kg/m. A total of 12,131 genetic variants across all of samples were found. After correcting for multiple testing, 65 SNPs from 25 nearest genes (INSIG1, FABP5, PTPRN2, VIPR2, WDR60, SHH, UBE3C, LMBR1, PAG1, IMPA1, CHMP4, SNX16, BLACE, EN2, CNPY1, LOC100506302, RBM33, LOC389602, LOC285889, LINC01006, NOM1, DNAJB6, LOC101927914, ESYT2, LINC00689) were associated with obesity at significant level q-value ≤ 0.05. LD block analysis showed there were 10 pairs of loci with D' ≥ 0.8 and r ≥ 0.8. GSEA further identified 2 major related gene sets, involving lipid raft and lipid metabolic process, with FDR values <0.12 and <0.4, respectively.Our data are the first documentation of genetic variants in 7q36.3 and 8q21.13 associated with obesity using target capture sequencing and Northern Han Chinese samples. Additional replication and functional studies are merited to validate our findings.

Introduction

There has been a worldwide epidemic in obesity linking increased morbidity and mortality as one of the major public health problems across countries. As a complex disorder, obesity is determined by both genetic and environmental factors, and the genetic influence accounts for 40% to 70% of the individual differences.[ Recent genome-wide association studies (GWAS) of BMI have identified 97 genetic variants.[ Nonetheless, these loci only explained 2.7% of the variance in BMI.[ While epistatic and gene–environment interactions may contribute to the unexplained heritability of obesity, there is possibility that a significant fraction of the missing heritability is due to loci not yet identified or fully characterized.[

The advent of next-generation sequencing (NGS) with high-throughput screening provided both a broad spectrum and a precise vision for the genetic architecture of many diseases.[ Although resequencing projects of whole human genome is still hampered by their high cost, the combination of target genomic region capture with NGS, as a low-cost technology with high efficiency and fidelity, has been used to investigate on several complex disorders and diseases.[

Our previous genome-wide linkage analysis on BMI in 126 dizygotic twins identified a genome-wide significant linkage peak on chromosome 7 with a log10 the odds ratio (LOD) score of 4.06 and 3 suggestive linkage regions on additional regions with LOD score ≥2.2.[ It is notable that our highest linkage region for BMI at 7q36 concurred with the result of a large multicenter linkage study of 4401 twin families from western countries.[

Greatly encouraged by these important findings, this study aims to determine the genetic variants on 7q36.3 and 8q21.13 (the top 2 ranked linkage regions in our previous study) involved in obesity through targeted resequencing technology.

Methods

Study samples

A total of 16 unrelated individuals with obesity were recruited from the Physical Examination Department of Qingdao Diabetes Hospital in October 2015. Information was collected through questionnaire, extraction of blood, together with anthropometric and laboratory measurements by well-trained clinicians face-to-face. BMI was derived by taking body weight (in kilogram) divided by height (in meter) squared. Subjects were included if the following criteria were met: aged 18 to 60 years; Han Chinese; ancestral home is in Shandong Province; BMI ≥ 30 kg/m2; free of hypertension, diabetes, or cardiovascular disease. Those who were pregnant, breastfeeding or taking weight-reducing medication within 1 month were excluded. Written informed consent form was obtained from all participants and the study was approved by the Qingdao University Ethics Committee.

DNA extraction, target genomic region capture and sequencing

Genomic DNA was extracted from whole peripheral blood of the subjects using QIAamp DNA Blood Mini Kit (Qiagen, GmbH, Hilden, Germany). DNA quantification and integrity were determined by the Nanodrop spectrophotometer (Thermo Fisher Scientific, Inc., Wilmington, DE) and the 1% agarose electrophoresis, respectively. A custom capture library (Agilent Technologies, Inc., Santa Clara, CA) on genomic locations of interest: 155100001–159138663 on chromosome 7 and 81874071–82674071 on chromosome 8 (Table 1) was designed. The information about the 2 target capture genomic regions from “UCSC genome browser” was shown in Supplementary Files 1 and 2.

Table 1

Target genomic regions.

Sixteen genomic DNA samples were captured on Agilent SureSelect custom library following the manufacturer's protocol (http://www.chem.agilent.com/library/usermanuals/Public/G7530–90000.pdf). Briefly, approximately 800 ng genomic DNA was sheared to 150 to 200 bp small fragments using sonicator (Covaris, Inc., Woburn, MA). The sheared deoxyribonucleic acid (DNA) was purified and treated with reagents supplied with the kit according to the protocol. Adapters from Agilent were ligated onto the polished ends and the libraries were amplified by polymerase chain reaction (PCR). The amplified libraries were hybridized with the custom probes. The DNA fragments bound with the probes were washed and eluted with the buffer provided in the kit. Then these libraries were sequenced on the Illumina sequencing platform (HiSeq X-10, Illumina, Inc., San Diego, CA) and 150 bp paired-end reads were generated.

Data analysis

We applied raw data filtering using Next Generation Sequencing Quality Control (NGSQC)—Toolkit v2.3.3 software. Raw reads which contained less than 70% high quality bases (Q20) or any N-base were removed. We then removed reads shorter than 70 bp afterwards to obtain clean reads. After quality control, raw data with 11.18 G bases were decreased to 10.41 G.

We used the Churchill software that integrated the following processes to call variants: Burrows–Wheeler aligner (BWA)—0.7.5a for mapping of paired end reads; Picard-tools—1.104 for marking duplicates that originate from PCR amplification (and that map at the same location); Genome Analysis Toolkit (GATK)—3.2 for realigning reads around indels, base recalibrate, and call variants with HaplotypeCaller method. Based on the recommended parameters on official website of GATK (https://software.broadinstitute.org/gatk/), the parameters stand_emit_conf was set to 30 to filter low quality variants.

We used ANNOVAR software to utilize up date-to-date information to functionally annotate genetic variants detected from Genome Reference Consortium GRCh37. And we identified variants documented in specific databases: the 1000 Genomes Project (http://annovar.openbioinformatics.org/en/latest/user-guide/filter/#1000-genomes-project-2015-aug-annotations) for allele frequency (AF) in populations; the LJB∗ databases for calculating SIFT scores, PolyPhen2 HDIV scores, PolyPhen2 HVAR scores, LRT scores, Mutation Taster scores, Mutation Assessor score, FATHMM scores, GERP++ scores, PhyloP scores, and SiPhy scores.

Based on the AF of each mutation in 16 samples and the AF in the 504 East Asians from the 1000 Genomes Project, we used R-3.0.0 software (http://web.mit.edu/people/jhaas/MacData/afs/sipb/project/r-project/arch/sun4x_59/lib/R/library/stats/html/fisher.test.html) for Fisher exact test. The resulting P-value and q-value were used as a basis for screening. The loci with q-value ≤ 0.05 were considered as statistically significant.

Linkage disequilibrium (LD) block analysis was performed for the multiple significant SNPs located near the same gene by using Haploview 4.2. A list of genes involved in statistically significant loci was then submitted to (http://software.broadinstitute.org/gsea/index.jsp) for gene set enrichment analysis (GSEA). False discovery rate (FDR) was calculated to obtain the significant gene sets.

Results

Basic information for 16 individuals with obesity is shown in Table 2. Through target genomic region capture sequencing, we obtained a total of 29,193,094 to 43,910,318 high-quality reads from these patients. The target region capture average ratio was 86.75%. All of the coverage ratios were over 92%, and the average depth of target regions was greater than 1500-fold. Therefore, sequencing coverage was fully adequate to detect gene variants within the majority of the targeted regions.

Table 2

Sex-specific characteristics.

A total of 12,131 genetic variants across the 16 samples were found. After comparison with the 504 East Asians in the 1000 Genomes Project, 65 SNPs involved 25 genes were associated with obesity significantly at q-values ≤ 0.05, as shown in Table 3 and Figs. 1 and 2.

Table 3

Sixty-five independent loci associated with obesity at q ≤ 0.05.

Figure 1

Manhattan plot showing the results of the association with obesity in 7q36.3.

Figure 2

Manhattan plot showing the results of the association with obesity in 8q21.13.

As considering several SNPs were found located near the same gene, LD block analysis was conducted and identified 10 pairs of loci with D′ ≥ 0.8 and r2 ≥ 0.8 (Table 4). The LD pattern of mutations among 65 significant SNPs is shown in Supplementary File 3.

Table 4

Linkage disequilibrium block analysis with D′ ≥ 0.8 and r2 ≥ 0.8.

To identify potential enriched gene sets, a genetic ontology (GO) pathway analysis was performed. The overrepresented gene sets (GO gene sets) were as follows: lipid raft, lipid metabolic process, phosphoric monoester hydrolase activity, cellular protein metabolic process, cellular macromolecule metabolic process, plasma membrane part, protein metabolic process, phosphoric ester hydrolase activity, plasma membrane, signal transduction, cellular lipid metabolic process, membrane part, hydrolase activity acting on ester bonds, as shown in Table 5. As for FDR, although all gene sets were above 0.05, the gene sets of lipid raft and lipid metabolic process had FDR values <0.4.

Table 5

Gene set enrichment analysis.

Discussion

This is our first attempt to explore the obesity-related SNPs in target genomic regions as informed from our previous genome-wide linkage study. We were able to identify 65 SNPs in 25 genes in association with obesity. On pathway level, 2 major gene sets were suggested: LIPID_RAFT and LIPID_METABOLIC_PROCESS.

Among the 25 genes, INSIG1, FABP5, PTPRN2, VIPR2 have been reported to be associated with obesity. INSIG1 encodes an endoplasmic reticulum membrane protein which regulates cholesterol concentrations in cell. Studies[ on INSIG1 gene and obesity suggested the gene plays a critical role in feedback regulation of lipid metabolism and may be involved in obesity development. FABP5 encodes the fatty acid binding protein found in epidermal cells and relevant pathways include glucose/energy metabolism. Canas et al[ have found overweight prepubertal boys showed elevated FABP5. A possible mechanism is that FABP5 regulates diet-induced obesity via GIP (gastric inhibitory polypeptide) secretion from enteroendocrine K cells in response to fat ingestion.[PTPRN2 plays a role in insulin secretion in response to glucose stimuli. Interestingly, most studies in Chinese have shown that the obesity-predisposing alleles were associated with insulin secretion, which is distinctive from the observations reported in Caucasians, among whom the obesity-related loci were primarily associated with insulin resistance.[ GO annotations related to VIPR2 include G-protein coupled receptor activity and vasoactive intestinal polypeptide (VIP) receptor activity. Note that none of the individual genes in the VIP pathway reached the genome-wide significance level in single-marker GWAS on obesity, however one study via pathway-based analysis of GWA-data[ suggested the VIP pathway was important for obesity.

Although there was no strong indication that WDR60, SHH, UBE3C, or LMBR1 polymorphism was the main causal variant of obesity in the population, studies showed that variation in these genes may be part of the multifactorial etiology of this complex condition. WDR60 encodes a member of the WD repeat protein family which plays a role in formation of cilia. The influences of cilia-related genes on adipogenesis via retrograde transport of SHH receptors and SHH signaling have been reported in a most recent review.[LMBR1 and UBE3C were known to be related to Coenzyme Q10 Deficiency, Primary, 2 and Kabuki Syndrome 1, respectively. The clinical symptoms for both diseases include obesity.

The rest of the obesity-related genes we have identified, including PAG1, IMPA1, CHMP4, SNX16, BLACE, EN2, CNPY1, LOC100506302, RBM33, LOC389602, LOC285889, LINC01006, NOM1, DNAJB6, LOC101927914, ESYT2, LINC00689 have not been reported. The biological function of these particular variants remains to be characterized.

Our data are the first documentation via target capture sequencing to identify obesity related rare variants in Northern Han Chinese. The highly interesting genomic regions were derived from our previous genome-wide linkage study in dizygotic twins from the same population, that could serve as important prior information. A major limitation of our study was relatively small sample size and it is desirable to replicate our findings in other studies.

In sum, based on our previous genome-wide linkage study, we identified genes and gene sets associated with adult obesity in 7q36.3 and 8q21.13 chromosome regions through targeted resequencing technology. We believe findings of this study contribute to further replication and functional studies.