Exome sequencing followed by genotyping suggests SYPL2 as a susceptibility gene for morbid obesity

Recently developed high-throughput sequencing technology shows power to detect low-frequency disease-causing variants by deep sequencing of all known exons. We used exome sequencing to identify variants associated with morbid obesity. DNA from 100 morbidly obese adult subjects and 100 controls were pooled (n=10/pool), subjected to exome capture, and subsequent sequencing. At least 100 million sequencing reads were obtained from each pool. After several filtering steps and comparisons of observed frequencies of variants between obese and non-obese control pools, we systematically selected 144 obesity-enriched non-synonymous, splicing site or 5′ upstream single-nucleotide variants for validation. We first genotyped 494 adult subjects with morbid obesity and 496 controls. Five obesity-associated variants (nominal P-value<0.05) were subsequently genotyped in 1425 morbidly obese and 782 controls. Out of the five variants, only rs62623713:A>G (NM_001040709:c.A296G:p.E99G) was confirmed. rs62623713 showed strong association with body mass index (beta=2.13 (1.09, 3.18), P=6.28 × 10−5) in a joint analysis of all 3197 genotyped subjects and had an odds ratio of 1.32 for obesity association. rs62623713 is a low-frequency (2.9% minor allele frequency) non-synonymous variant (E99G) in exon 4 of the synaptophysin-like 2 (SYPL2) gene. rs62623713 was not covered by Illumina or Affymetrix genotyping arrays used in previous genome-wide association studies. Mice lacking Sypl2 has been reported to display reduced body weight. In conclusion, using exome sequencing we identified a low-frequency coding variant in the SYPL2 gene that was associated with morbid obesity. This gene may be involved in the development of excess body fat.


INTRODUCTION
Obesity is a major contributor to the global burden of chronic disease, such as type 2 diabetes. 1 A strong genetic impact on obesity has been established. However, the identification of specific genes involved in common forms of human obesity has proven elusive, despite great efforts made through different approaches, including candidate gene studies, genome-wide linkage and genome-wide association (GWA) studies. 2,3 So far, GWA studies have identified around 50 obesityassociated genetic loci. However, these loci together explain only about 1-2% of body mass index (BMI) variations in the studied populations. 4 Therefore a major question is how the remaining so called 'missing' heritability can be explained. 5 There are two principal concerns about the interpretation of GWA results. First, the approach is based on the common disease, common variant hypothesis. The commercial single-nucleotide polymorphism (SNP) arrays used in GWA studies were designed to cover the common genetic variants (minor allele frequency (MAF) 45% in populations). Therefore rare variants could not be directly evaluated using such arrays. Second, imperfect statistical models, which do not account for gene-gene and gene-environment interactions, were used in GWA data analyses. 6 Recently developed high-throughput sequencing technology complements SNP arrays and shows power to detect rare disease-causing variants by deep sequencing of all known exons. [7][8][9] This revolution of technology holds promise to elucidate complex diseases by allowing a genome-wide search for low frequency and rare and putatively functional variants. Exome sequencing might be most efficient when applied to patients with more extreme forms of common diseases, such as morbid obesity; first, it increases the chances to obtain significant association for rare variants with high penetrance; second, genetic variants in the protein coding sequence might be more likely to have a strong impact on the phenotype than variants in intergenic regions.
In this study, we employed exome sequencing to detect low frequency and rare gene variants enriched in adult subjects with morbid obesity and validated them against never-obese elderly adults. We reasoned that in this way we would enrich for obesity genes among cases and filtered away such genes in the controls. Here we report a low-frequency obesity-associated variant in the coding region of the synaptophysin-like 2 (SYPL2) gene. or in association with planned visits to our medical or surgical units for morbid obesity or scoliosis, respectively. Most of them have been described. 10,11 From the obese subset, we selected, among the 10% with the highest BMI, 100 subjects with extreme obesity for exome sequencing (69 women, 31 men, age 41.5 ± 11.5 years, BMI 52.2 ± 3.8 kg/m 2 ). All selected subjects had morbid obesity defined as BMI440 kg/m 2 . As controls for the exome sequencing, we used subjects already collected for an ongoing genetic study on idiopathic scoliosis (76 women, 24 men, age 24.5 ± 12.8 years, BMI 21.9 ± 4.3 kg/m 2 ). For the first genotyping validation, we selected the 494 morbidly obese subjects with the highest BMI in the large sample collection described above, including the 100 subjects who had been used for exome sequencing. As controls, we used 496 never-obese subjects with age 440 years and BMI o30 kg/m 2 . In the second confirmation, we genotyped remaining 1425 morbidly obese subjects and 782 never-obese controls with age 440 years and BMIo30 kg/m 2 .
The non-obese controls for exome sequencing had scoliosis. Otherwise all other control subjects were healthy according to self-report. There was an overrepresentation of women in the studied cohorts. Obese women are more likely to search medical advice for their obesity. Scoliosis is more common among women. Subjects were of European ancestry and living in Sweden. The study was approved by the local Ethics Committees, and all subjects gave their informed consent to participation.

DNA preparation and pooling
Genomic DNA was prepared from PBMCs using the QiAmp DNA blood Maxi kit (cat no. 51194, Qiagen, Hilden, Germany). DNA purity and quality was confirmed by A260/280 ratio 41.8 on Nanodrop (Thermo Fisher Scientific Inc., Waltham, MA, USA) and agarose gel electrophoresis. DNA concentration was measured by Qubit (Life Technologies, Stockholm, Sweden). Subsequently, we took 0.8 μg of each DNA sample and randomly divided them into 10 pools, each containing 10 samples of either obese cases or controls. The concentrations of pooled DNA samples were measured with Qubit, and the samples were run on agarose gel.

Exome sequencing
Exome sequencing was performed at the Science for Life Laboratory (SciLife-Lab), Stockholm, Sweden. Each DNA library was prepared from 3 μg of the pooled genomic DNA. DNA was sheared to 300 bp using a Covaris S2 instrument and enriched by using the SureSelectXT Human All Exon 50 Mb kit and an Agilent NGS workstation according to the manufacturer's instructions (SureSelectXT Automated Target Enrichment for Illumina Paired-End Multiplexed Sequencing, version A, Agilent Technologies, Santa Clara, CA, USA).
The clustering was performed on a cBot cluster generation system using a HiSeq paired-end read cluster generation kit according to the manufacturer's instructions (Illumina, San Diego, CA, USA). The samples were sequenced on an Illumina HiSeq 2000 as paired-end reads to 100 bp/read (Illumina). All lanes were spiked with 1% phiX control library, except for lane 8, which had 2% phiX. The sequencing runs were performed according to the manufacturer's instructions. Base conversion was done using Illumina's OLB v1.9 (Illumina).

Variant filtering and enriching for obesity
We used the following filtering criteria to select low-frequency and rare singlenucleotide variants (SNVs) from exome sequencing data. Those SNVs enriched in morbidly obese subjects were taken forward for validation ( Figure 1). First, we filtered for variants with a depth of 5 × and mapping quality (MQ) ≥ 20. Then we looked for putatively functional variants, that is, variants from exonic regions, splicing sites or 5′ upstream regions. In the next step, we compared the occurrences of a variant between obese and control pools. Only variants called in ≥ 2 obese pools, but with no call or called once in control pools, were used in the further filtering step; this step was applied to avoid potential false-positive variants. The last filtering was based on allele frequency. We looked for lowfrequency and rare SNVs, that is, variants that were not found in public databases or known SNPs with a MAF ≤ 5% in general populations (1000 Genomes project 2011 May release).

Variant validation by genotyping
Variant validation was focused on SNVs enriched by the above filtering steps that could fit into a manageable genotyping panel and was carried out in two steps. First, all selected variants were genotyped in 494 obese subjects and 496 controls (Table 1 and Figure 1) using the Illumina GoldenGate assay performed at the SNP/SEQ Technology Platform at Uppsala University, Sweden (http://www.molmed.medsci.uu.se/SNP+SEQ+Technology+Platform/Genotyping/, Fan et al 14 ). Secondly, genotyped variants showing nominal significant association with obesity were genotyped in an additional 1425 obese and 782 controls. Genotyping was performed using matrix-assisted laser desorption/ ionization time-of-flight mass spectrometry (SEQUENOM, Agena Bioscience, San Diego, CA, USA) at the Mutation Analysis core Facility (MAF) at Karolinska Institutet, Sweden. 15 Multiplexed assays were designed using the MassARRAY Assay Design v4.0 Software (Agena Bioscience). Protocol for allele-specific base extension was performed according to Agena Bioscience's recommendation. All genetic association results in validation cohorts 1 and 2 have been submitted to GWAS Central.

Statistical analysis of association
Subjects were characterized by BMI, and values are reported as mean ± SD for obese and controls separately. Genetic association analysis and odds ratio calculation were performed using PLINK (http://pngu.mgh.harvard.edu/p urcell/plink/, Purcell et al 16 ). Hardy-Weinberg equilibrium (HWE) of the genotype frequencies among cases and controls were checked separately for each run of validations prior to association analysis. A HWE nominal P-value of Values for BMI and age are mean ± SD. Controls used in the first validation and second confirmation were non-obese.
SYPL2 as an obesity-associated gene H Jiao et al 0.01 in controls was used as a cutoff for exclusion of variants from further analysis. BMI was analyzed using standard linear regression implemented in PLINK. Sex and age were used as covariates to evaluate their influence on the association of a variant with BMI.

Discovery of rare variants by exome sequencing
To identify potential functional low-frequency and rare variants associated with extreme obesity, we performed exome sequencing on pooled DNA from 100 subjects (10 pools) with severe morbid obesity and 100 control subjects (10 pools). Thus there was DNA from 10 subjects in each pool. The characteristics of subjects for the exome sequencing are shown in Table 1. The obese subjects were older than the controls, whereas the gender distribution was similar between cohorts, with an overrepresentation of women. In total, we obtained 4100 million reads from each pool (Supplementary Table S1). After trimming low-quality reads and removing PCR duplicates, 95% of reads were mapped to the current human genome reference (assembly hg19, NCBI build 37) except for one control pool that had a little lower mapping rate. The vast majority of reads were uniquely mapped on the same chromosome for both paired reads. Discordance between paired reads was rare. Reads were paired with higher rates (97%) in obese pools than in control pools (80-95%) (Supplementary Table S1).
More than 92% of target regions were covered by 5 × reads, and 480% of target regions were covered by 30 × reads (Supplementary  Table S2).
Potentially interesting variants were identified through several filtering steps (Figure 1). All variants were called using a depth of 5 × and MQ ≥ 20 as the first filtering criteria. In total, we got 114 034-190 130 variants from each of the 20 pools. About 95% of variants were known SNPs, that is, SNPs with rs numbers. The majority of variants were located in exons, introns or intergenic regions, whereas others were in nearby slicing sites or in regulatory regions (Table 2). We focused on putatively functional variants, that is, variants from exonic regions (excluding synonymous variants), splicing sites or 5′ upstream. We obtained 31 520 such variants in obese pools and 26 565 in non-obese pools. Those variants were found in at least one obese pool or control pool, respectively. In all, 20 652 (65.5%) potentially functional variants were present in both obese and control pools, and 10 868 (34.5%) were singletons for obesity. Next we looked for variants appearing more often in obese pools. At this step, there were 5788 variants left, and they were detected in ≥ 2 obese pools but did not appear or appeared just once in control pools. Finally we filtered for low-frequency and rare SNVs, that is, SNVs that were not found in public databases or SNPs with a MAF ≤ 5% in

SNV validation and association analysis
From the 1032 SNV, we excluded insertion or deletion variants to avoid potential technical difficulties in genotyping. Then we filtered for SNVs shared by at least two obese pools with the highest yield. After manually checking alignments of reads carrying the variants by using the IGV software (http://www.broadinstitute.org/igv/) to exclude potential artifacts, we selected 144 SNVs for genotyping (142 dbSNP and 2 unknown SNVs). Sequencing depths and numbers of shared obesity pools of the 144 SNVs are shown in Supplementary Table S3. Most SNVs were found in two (83) or three (31) pools. Ten SNVs had difficulties in primer design for genotyping and were replaced with singleton SNVs. SNV validations were divided into two steps. First, the 144 SNVs were genotyped in 494 morbidly obese subjects, including the 100 exome-sequenced obese subjects, and 496 non-obese controls ( Table 1). Genotyping of 20 SNVs failed. The average per SNV call rate for the remaining 124 SNVs was 97.95%. The concordance rate of MAFs estimated from exome sequencing data and the MAF based on genotype results was 0.83 (Supplementary Figure S1). Among the successfully genotyped SNVs, 21 were monomorphic. The remaining 103 variants were used for association analysis. HWE was checked prior to association analysis. After removal of 5 additional SNVs due to an HWE P-valueo0.01, 98 SNVs were used in further association analysis.
Five out of the 98 SNVs showed significant allelic association with obesity with a nominal P-value 0.05 (Table 3, Supplementary Table S4). These five SNVs were subsequently subjected to the second run of genotyping for confirmation in a large case-control cohort consisting of 1425 morbidly obese subjects and 782 never-obese controls. Only one SNV, rs62623713 (NM_001040709:c.A296G:p.E99G), consistently displayed a significant association with obesity (Table 3). However, in the final joint analysis of all genotyped subjects two SNVs, rs62623713 and rs35923425 (NM_024727:c.G1134C:p.L378F), showed nominally significant association with obesity in 1917 morbid obese and 1278 never-obese controls, with a P-value of 0.0063 (odd ratio (OR) = 1.32, MAF 0.081 in obese and 0.063 in controls) or a P-value of 0.004 (OR = 1.35, MAF 0.076 in obese, 0.058 in controls), respectively (Table 3). rs62623713 located at 110 019 439 base pair (bp) on chromosome 1 (hg19) is a missense variant in exon 4 of the SYPL2 gene. rs35923425 located at 169 569 432 bp on chromosome 3 is also a missense variant in exon 7 of the leucine-rich repeat containing 31 (LRRC31) gene.
By quantitative trait analysis in all genotyped subjects, rs62623713 showed strong association with BMI assuming an additive genetic model (beta 2.13, P-value of 6.28 × 10 − 5 , Table 4). The variant explained 0.5% BMI variation in genotyped subjects. Homozygous subjects for the rs62623713*G allele had higher BMI (mean = 45.56 kg/m 2 ) when compared with subjects who were heterozygous or homozygous for the rs62623713*A allele (Supplementary Table S5). SYPL2 as an obesity-associated gene H Jiao et al The G allele is overrepresented in our obese subjects (MAF 8%), especially in the homozygous form (12 obese to 2 controls). The same tendency of associations was shown in females (beta 1.99, P-value 0.001), males (beta 2.42, P-value 0.016), obese subjects (beta 3.39, P-value 0.0007) and non-obese subjects (beta 2.16, P-value 0.031) ( Table 4). rs35923425 in LRRC31 was also associated with BMI (beta 1.5, P-value 0.60 × 10 − 3 ). However, significant association was seen in females only. The significance of associations of 5 SNVs with either obesity or BMI in the first and the final analysis were not influenced when we excluded 100 subjects used in exome sequencing for analyses (Supplementary Tables S6 and S7).
To further evaluate any confounding effect of age or gender, we added these variables as covariates in regression models and tested their individual effects separately. rs62623713 maintained strong associations with BMI after adjustment for age or gender (Supplementary Table S8). The association between rs35923425 and BMI disappeared after adjustment for age. Furthermore, we also performed a variety of tests to adjust for the analysis of multiple genetic variants. The association of rs62623713 with BMI was still significant with a Bonferroni adjusted P-value of 3.0 × 10 − 4 and a genomic-control corrected P-value of 0.05. The association of rs35923425 became non-significant with genomic controls, whereas the Bonferroni corrected P-value was 0.029 (Supplementary Table S9).

DISCUSSION
Using exome sequencing followed by large-scale genotyping, we identified a low-frequency coding variant, rs62623713 (E99G) in exon 4 of the SYPL2 gene, consistently showing association with morbid obesity. rs62623713 has a MAF of 2.9% for the G allele in the general populations according to 1000 Genomes. The variant is overrepresented among our obese subjects (MAF 8%). rs62623713 was not covered by Illumina or Affymetrix genotyping arrays used in previous GWA studies, and no other variants on those arrays was in strong linkage disequilibrium with rs62623713.
To the best of our knowledge there are only three reports of obesityassociated variants detected by exome sequencing. 9 The hypothesis underlying our study is that more severe forms of common obesity might be due to low-frequency or rare variants with a larger impact on the phenotype. A minor portion of morbidly obese cases has been shown to be due to genetic variants with high penetrance. Variants in the melanocortin 4 receptor (MC4R) gene explain a few percentage of childhood obesity, 19 and a copy number repeat on chromosome 16 20 has a high penetrance in a few percentages of morbidly obese subjects with cognitive defects. We did not detect any variants in exonic or nearby upstream regions of the MC4R in our exome sequencing data, which may be due to the low number of investigated subjects. Our results, that is, the detection of one SNP in SYPL2 remaining significantly associated with obesity after adjustment for multiple testing, support the notion that low-frequency and rare variants contribute to morbid obesity. However, the present study lacks power to estimate to what extent rare genetic variants cause morbid obesity in the population.
The approach we report using pooled DNA for exome sequencing to get obesity-associated variants, which was followed by genotype validation, is a cost-effective method. Our study has its limitations. We did not have access to a reliable family history of our patients and therefore do not know whether their obesity displays patterns of monogenic inheritance. On the other hand, as obesity is common, pedigree information might be difficult to interpret due to heterogeneity, that is, different causes of obesity in different family members. We only genotyped validated variants detected in more than one sequencing pools to limit the number of false positives. Although by  Figure S1). In our study, failure in validation by genotyping seemed not due to low depth as major failed SNVs had an average depth ≥ 20 × (Supplementary Table S3). Imperfect primer designs in multiplex genotyping could be a major cause of genotyping failure. SYPL2 (also called MG29) is expressed in the adipose tissue, brain, kidney, heart and cerebellum. Functionally, the SYPL2 protein has been proposed to participate in cellular calcium ion homeostasis. 21 The protein is expected to be involved in transporter activity and to localize in various compartments, such as synaptic vesicle. 22 Interestingly, mice lacking Sypl2/Mg29 exhibit reduced body weight, abnormal skeletal muscle membranes and irregular skeletal muscle contractility. 23 How the gene might influence adipose tissue formation and its potential role in obesity development is unclear and needs further investigation, which is beyond the scope of this study.
In previous reports, SYPL2 was associated with major depressive disorder in European population. 24 An association between obesity and depression has repeatedly been established. In a meta-analysis, obesity was found to increase the risk of depression. In addition, depression was found to be predictive of obesity. 25 Another study found that the association of obesity with depression was mainly confined to persons with severe obesity. 26 SYPL2 belongs to the synaptophysin family. Synaptophysin regulates activity-dependent synapse formation in cultured hippocampal neurons, 27 and it is required for kinetically efficient endocytosis of synaptic vesicles in cultured hippocampal neurons. 28 Food intake is subject to a complex regulation by the hypothalamus and other brain centers, including the brain stem and the hippocampus. Thus one could hypothesize that SYPL2 is involved in the central regulation of food intake possibly affecting the central reward systems and the hedonic effects of food.
In conclusion, we identified a low-frequency obesity-associated coding variant in the SYPL2 gene using exome sequencing with pooled DNA from 100 morbidly obese and 100 non-obese subjects, which was followed by genotype validation in 3197 case-control subjects. Our results provide evidence of the existence of a coding variant associated with obesity although further functional studies of this genetic variant remain to be performed.