Systems genetics analysis of obesity in a porcine model using WISH network method

High-throughput genotype data have been extensively used to investigate the biological and genetic background of complex traits, such as obesity. However, results explaining complete genetic variation are limited. Systems genetics approaches are increasingly used to elucidate more of the currently limited knowledge of complex traits. Weighted Interaction SNP Hub (WISH) network method is a novel tool which is able to elucidate underlying biology and capture potential causal variants. Here we studied obesity, using systems genetics approaches in an F2 pig population. We created an Obesity Index (OI) based on multi-trait selection indexes containing nine obesity-related phenotypes, to perform a GWAS and WISH network analysis. WISH detected several modules associated with OI. Functional annotation revealed several genes (e.g. NPC2) and pathways (e.g. the insulin signaling pathway) which were (in-)directly related to obesity. This study shows the potential for systems genetics analysis of high-throughput genotype data for obesity. Introduction High-throughput genotype data (HTG) are extensively used for complex traits and diseases in genomewide association studies (GWAS) to find genes associated with the trait of interest. Unfortunately, to date the detected loci explain a limited amount of the genetic variation in many complex traits (Manolio et al. (2009)), as well as do not provide holistic understanding of the systems genetics. Obesity is a complex trait that is associated with various serious diseases (e.g. type 2 diabetes and various types of cancer). Because of the exponential rise in the prevalence of obesity, and its financial, social and welfare consequences, the urge for a decent biological and genetic knowledge is rising. Here, we use a porcine model for human obesity, since the pig has similar metabolic, digestive and cardiovascular features, and it resembles humans more than rodents (Spurlock and Gabler (2008)). We have created an F2 pig population, and demonstrated the potential of this population to study human obesity (Kogelman et al. (2013)). In complex diseases, it has been indicated that gene x gene interactions may have a key role, which is not taken into account in GWAS (Cordell (2009)). Systems genetics approaches, which include the effect of interacting alleles and their functional/biological role, may reveal more of the biological and genetic background of complex traits and diseases (Kadarmideen et al. (2006)); Civelek and Lusis (2013); Do et al. (2014)). Several network approaches have been developed to understand the interaction of genes and their influence on complex diseases (Segal et al. (2003); Diez et al. (2010)). Recently, we published the Weighted Interaction SNP Hub (WISH) network method based on HTG, which makes it possible to identify clusters of highly interconnected SNPs (modules) that are related to disease risks or a phenotype (Kogelman and Kadarmideen (2014)). The WISH network method detects those modules using the interconnectedness between SNPs, based on genotype interactions. Functional annotation and pathway analysis of detected modules may lead to the identification of biologically relevant pathways underlying those phenotypes. In the present study, we applied systems genetics approaches in the F2 pig resource population, to unravel pathways and detect potential causal obesity genes. This included a GWAS and WISH on OI; and pathway profiling and functional enrichment analyses on the resulting variants and SNP modules. These systems genetic approaches revealed novel biological and genetic determinants not otherwise identified only by GWAS. Materials and Methods Experimental design. We have created an F2 pig resource population (Kogelman et al. (2013)), by generating an intercross of two production breeds (Yorkshire and Duroc) and the Gӧttingen Minipig. While the production breeds are prone to leanness and muscle growth due to genetic selection, the minipigs are prone to obesity and share metabolic impairments seen in obese humans (Johansen et al. (2001)). The F2 population (454 pigs) was intensively phenotyped (e.g. weight, body confirmation, DXA scanning and slaughter characteristics) and genotyped using the Illumina Porcine 60K SNP Chip. Obesity Index. Multi-trait selection indexes, based on estimated breeding values (EBVs) for several traits of interest, are intensively used in animal breeding to determine extreme phenotypes (Cameron (1997)). We created the Obesity Index (OI) to determine the degree of obesity of the F2 animals, by calculating selection index weights and combining the EBVs for nine different obesityrelated traits (reported in Kogelman et al. (2013)) into one aggregate total merit index. Traits selected for the OI were: weight and abdominal circumference at slaughter, average daily gain, estimated fat mass and percentage of fat at DXA scanning, back fat thickness at position 1 and position 2, weight of leaf fat and weight of omental fat at slaughtering. The OI followed, as expected, a normal distribution among the F2 animals. Genome-wide association analysis. A GWAS was performed using the R package GenABEL (Aulchenko et al. (2007)) to test associations between the SNPs and the OI. We ran the basic linear model without correcting for any effects, as the OI is constructed based on EBVs which are already corrected for the population structure and other fixed effects: y= μ+g+e, where y = OI, μ = the phenotypic mean, g = the SNP genotypes (codes as 0,1, and 2) and e = the model errors. We calculated the Bonferroni corrected pvalues resulting in a suggestive association at Padj = 1.32E (0.05/SNP number) and a highly significant association at at Padj = 2.64E (0.001/SNP number). Combined linkage disequilibrium linkage analysis. Findings of highly significant regions in the GWAS were confirmed using a combined linkage disequilibrium linkage analysis (LDLA) approach, using the information of the intercross pedigree (paper presented in this conference by Pant et al.). Statistical significance was calculated via a likelihood ratio test of the full model (OI regressed over the phenotypic mean and Identity by descent probabilities of chromosomal segments flanked by successive marker pairs) vs. the null model (containing only the phenotypic mean). The Identity by descent probabilities were estimated using a linkage disequilibrium multi-locus iterative peeling algorithm described in Meuwissen and Goddard (2010). WISH network construction. We previously published and described the WISH network method (Kogelman and Kadarmideen (2014)). This method assumes that SNPs which are highly correlated with each other (clusters of SNPs, here termed modules), will work co-operatively in pathways. Because of computational limits the top 2500 SNPs were selected based on their genome-wide significance (GWAS) and connectivity (sum of correlations of a particular SNP with all other SNPs). The adjacency matrix was calculated as the Pearson’s correlations between SNP genotypes, raised to the power of 5 to create a scale free network. A SNP dendrogram was created based on the topological overlap measure (TOM), representing the relatedness between SNPs, to identify clusters of highly interconnected SNPs (modules). The correlation between the eigenvalue of each module (1 principal component) and the phenotypes of interest (OI and other obesity related phenotypes) was used to create the Genome-wide Module Association Test (GMAT) for selection of potential biologically interesting modules. Systems biology analyses. All significant GWAS top SNPs and identified SNP modules by the WISH network were further analyzed by detecting overrepresented gene ontology (GO) terms and pathways. Genes located at the detected SNPs were obtained using the NCBI2R Rpackage (available at http://cran.r-roject.org/web/packages/ NCBI2R/), using a flanking distance of 20 kB to cover the promotional regions of genes. For downstream analyses we used GOEAST (Gene Ontology Enrichment Analysis Software Toolkit, Zheng and Wang (2008)) and GeneNetwork (http://genenetwork.nl). In GOEAST, the Gene Batch tool was used to import the gene symbols, and following the hypergeometric statistical test, with the Yekutieli multipletesting adjustment method, was performed to identify significantly overrepresented GO terms and corresponding pathways. GeneNetwork was used to identify overrepresented GO terms, KEGG pathways, phenotypes and tissues. The database is constructed using human, mouse and rate expression data. Gene functions were predicted against known pathways and gene sets in various biological databases. Overrepresentation of GO-terms and pathways was tested using the Mann-Whitney U test, and Pvalues were afterwards corrected for multiple testing using the Bonferroni correction. Results and Discussion Genome-wide association analysis. The GWAS results were presented in a Manhattan plot (Fig. 1), showing 1024 SNPs above the suggestive significance threshold (Padj=1.32E) and 404 SNPs above the highly significance threshold (Padj=2.64E). Figure 1. A Manhattan plot of GWA study single-locus P-values. The blue dash line indicates a suggestive significance threshold with adjusted Bonferroni correction at Padj=1.3E and the red line indicates a highly significant threshold with adjusted Bonferroni correction at Padj=2.6E. Three of the most significant SNPs are rs81396056 (P=9.09E), rs81238148 (P=1.13E) and rs81416774 (P=2.56E), which are situated within the NPC2 gene, OR4D10 gene and CACNA1E gene, resp. The NPC2 gene is encoding a protein that plays a role in the regulation of cholesterol transport through the late endosomal/lysosomal system, affecting cholesterol homeostasis (Storch and Xu (2009)). The OR4D10 gene is an olfactory receptor gene, which are responsible for the perception of smell, through neuronal responses, and thereby affecting eating behavior (Nasser (2008)). The CACNA1E gene encodes a protein in a voltage depen