Epistatic Effects on Abdominal Fat Content in Chickens: Results from a Genome-Wide SNP-SNP Interaction Analysis

We performed a pairwise epistatic interaction test using the chicken 60 K single nucleotide polymorphism (SNP) chip for the 11th generation of the Northeast Agricultural University broiler lines divergently selected for abdominal fat content. A linear mixed model was used to test two dimensions of SNP interactions affecting abdominal fat weight. With a threshold of P<1.2×10−11 by a Bonferroni 5% correction, 52 pairs of SNPs were detected, comprising 45 pairs showing an Additive×Additive and seven pairs showing an Additive×Dominance epistatic effect. The contribution rates of significant epistatic interactive SNPs ranged from 0.62% to 1.54%, with 47 pairs contributing more than 1%. The SNP-SNP network affecting abdominal fat weight constructed using the significant SNP pairs was analyzed, estimated and annotated. On the basis of the network’s features, SNPs Gga_rs14303341 and Gga_rs14988623 at the center of the subnet should be important nodes, and an interaction between GGAZ and GGA8 was suggested. Twenty-two quantitative trait loci, 97 genes (including nine non-coding genes), and 50 pathways were annotated on the epistatic interactive SNP-SNP network. The results of the present study provide insights into the genetic architecture underlying broiler chicken abdominal fat weight.


Introduction
Epistasis -the interactions between polymorphic loci, such as SNPs, genes or quantitative trait loci (QTLs) -is a hot topic in quantitative genetics [1]. Epistasis is a major factor that determines variation in quantitative phenotypes [2]. Mapping the epistatic loci affecting quantitative traits will advance our understanding of the genetic architecture of these complex phenotypes for humans and model organisms [2,3]. Recently, the study of epistasis among genes, SNPs or QTLs has progressed rapidly [4][5][6][7][8], as has the study of epistasis in biological experiments [9].
Much progress has been made using genome-wide association studies (GWASs). Performing genome-wide SNPs interaction analysis represents the next step for detecting the variations of quantitative traits, because single-locus tests cannot identify the interactions among SNPs, genes or other genetic or environmental factors [10]. The markers detected by GWAS only explained a fraction of the heritable variance [11][12][13]. Identifying markers that show interactions would help to explain a higher proportion of the heritable variance. From sequence to phenotype, the SNP-SNP and gene-gene interactions can provide new insights into the genetic basis of complex traits.
In chickens (Gallus gallus), it has been suggested that epistatic interactions between genes (or QTLs) are important for quantitative traits such as growth and fatness traits [14][15][16]. Abdominal fat in meat-type chicken is a quantitative trait, and is an adverse economic factor; therefore, many QTLs and SNPs affecting fatness traits have been identified [17][18][19][20][21]. A previous study on chicken abdominal fat traits identified epistatic interactions among 10 candidate genes, and constructed networks of the interacting genes [15]. However, a large proportion of the effects of the epistatic interactions between polymorphic loci on fatness traits remain undetermined, and the whole genetic network is unclear. Hence, it is essential to study the effects of interactions between SNPs, genes and QTLs on fatness traits.
In this study, we aimed to determine the interactions between polymorphic loci affecting chicken abdominal fat weight (AFW). Pairwise epistatic interaction effects among genome-wide SNPs were detected. A network of significant SNP pairs was analyzed and annotated. The results provide further insight into the genetic network controlling abdominal fat deposition in chickens.

Ethics Statement
All animal work was conducted according to the guidelines for the care and use of experimental animals established by the Ministry of Science and Technology of the People's Republic of China (Approval number: 2006-398) and was approved by the Laboratory Animal Management Committee of Northeast Agricultural University.

Experimental Populations
The experimental chickens were from the Northeast Agricultural University broiler lines divergently selected for abdominal fat content (NEAUHLF) [22]. A total of 475 male individuals, 203 in the lean line and 272 in the fat line, derived from the 11 th generation population of NEAUHLF, were used in this study.

SNP Genotyping
Genotyping was carried out using the Illumina Inc. (San Diego, CA, USA) chicken 60 K SNP chip, which contains 57,636 SNPs. Markers with minor allele frequencies less than 5% or monomorphic loci were filtered out. Individuals were removed who had 5% or more missing SNP genotypes. Finally, 45,611 SNPs in 475 individuals were used for the pairwise interaction analyses. SNP genotypes were provided by Zhang et al. [23].
The EPISNP's statistical model for testing epistatic effects for AFW is where y is the dependent variable (AFW), m is the population mean, SNP 1 and SNP 2 are the two single-locus genotypic effects, SNP 1 6SNP 2 is the two-locus interaction effect, F (family effect) is a random effect, and e is the random error.
The two-locus interaction effect was partitioned into four individual epistatic effects using the extended Kempthorne model, which allows Hardy-Weinberg disequilibrium and linkage disequilibrium (LD): Additive6Additive (AA), Additive6Dominant (AD), Dominant6Additive (DA), and Dominant6Dominant (DD) epistatic effects [24]. An F-test was used to test the significance of the two-locus interaction effect.
The P-value of model (1) was adjusted by a Bonferroni correction (4.16610 9 ; independent tests using the same data set, with a significance threshold of P,0.05), and the significance threshold of the test was determined as P,1.20610 211 .

Contribution Rate Calculation
The contribution rate of every significant epistatic interactive SNP pair to phenotypic variation was calculated by the formula, where SS SNP1 6 SNP2 is the variance of the significant SNP 1 6SNP 2 interactive effect (AA, AD, DA or DD), and Var p is the phenotype variance.

SNP-SNP Network
The figure showing the SNP-SNP network with significant epistatic effects was drawn using the Cytoscape 2.8 software package [25]. The interactive SNPs affecting AFW with a pairwise interaction P-value of P,1.20610 211 were loaded into Cytoscape to visualize the network. Based on entropy theory [26], the importance of the different subnets was evaluated by the formula, where w is the importance of the subnet, n is the edge number (i.e. the number of SNP pairs), p i is the P-value of the interactive effect being tested for the i th SNP pair, and c i is the contribution rate of the i th SNP pair.

Annotation of SNP-SNP Network
The purpose of this step was to annotate significant SNPs and mine QTLs or gene networks affecting abdominal fat traits. First, the significant interactive SNPs were put into the Chicken QTL database (http://www.animalgenome.org/cgi-bin/QTLdb/GG/ i-ndex), and mapped in terms of QTLs. Second, the network was annotated with genes and pathways. A mean value of r 2 .0.8 [27] was observed in pairwise distances of approximately 0.2 Mb, as calculated by HAPLOVIEW v4.1 [28], and genes whose physical distance to the significant SNP was smaller than 0.2 Mb were selected. The positions of genes were acquired from UCSC Genome Bioinformatic site (http://genome.ucsc.edu/). Genes in the same network were put into the KEGG program (http://www. genome.jp/kegg/) to identify pathways that might affect the traits.

SNP Makers
A total of 45,005 SNPs on 28 autosomes, the Z chromosome, linkage groups and 606 SNPs not assigned to any chromosomes in chickens were included in this study (Table 1). These markers covered 1026.23 Mb of the genome, with an average of 16.63 kb between adjacent markers.

Pairwise Interaction Effects Analysis
EPISNP3 [24] was used to analyze the pairwise interaction effects among genome-wide SNPs for AFW across both lines. Fiftytwo pairs of significant SNPs (Table 2), containing 45 pairs of AA and seven pairs of AD epistatic effects, were detected. No DA or DD epistatic effects were detected. The P-value of the epistatic effect between SNP Gga_rs13569377 on GGA18 (Gallus gallus, GGA) and Gga_rs14988623 on GGA13 reached 2.54610 214 , which was the most significant effect detected. There were six pairs of SNPs on the same chromosome (one pair on GGA2, four pairs on GGA3 and one pair on GGA10). The other 46 pairs of SNPs were distributed on different chromosomes. The contribution rate of the significant epistatic interactive SNP pairs ranged from 0.62% to 1.54%, with 47 pairs having a contribution rate of more than 1%. The 52 pairs of significant SNPs comprised 68 single SNPs covering 18 chromosomes. The position of each SNP is shown in Table S1. Many SNPs were detected to interact with the same locus; for example, seven SNPs on GGA3 all had epistatic interactions with Gga_rs14303341 on GGA27, and four SNPs on GGA18 interacted with Gga_rs14988623 on GGA13. The SNP pairs having the same SNPs are shown in Figure 1.

SNP-SNP Network Analysis
To investigate the complex mechanism of epistatic effects on AFW, a network of SNPs having epistatic interaction affecting AFW was constructed. SNP epistatic interaction subnets containing more than three nodes are shown in Figure 1.
The degree of nodes refers to the number of edges connecting it, and an edge stands for a two-way interaction. For example, the degree of SNP Gga_rs14988623 was 4, and the degree of SNP Gga_rs14303341 was 7 ( Figure 1). SNP Gga_rs14303341 on GGA27 with the most two-way interactions, whose degree was the maximum, could be seen as the hub site of subnet A in the network of epistatic SNPs. SNP Gga_rs14988623 on GGA13 connected to four SNPs on GGA18 in subnet E was also a hub site. Subnet C contained eight edges, which was the maximum. The eight edges all happened between GGAZ and GGA8, so subnet C implies that that the interaction occurs between the two chromosomes. Interactions between two chromosomes were also hinted at by subnets A, D, E, F, G, and H. Interactions within the same chromosome were suggested by subnet I. Subnet B contained SNP pairs from four chromosomes, which was the largest number of chromosomes in any subnet.
Interestingly, in some subnets, the SNPs were adjacent to one another on the same chromosome. For example, in subnet A, seven SNPs on GGA3 were physically close to one another. This phenomenon was also observed for subnets C, D, E, G and H. Thus, there seemed to be a concentration of adjacent SNPs in these subnets. This indicated that the interaction takes place in the corresponding regions, especially for subnet C.
The importance of the subnets, as evaluated by formula (3), is shown in Table 3. The subnets with higher scores are more important. According to formula (3), subnets should get higher scores and be more important if they have more significant Pvalues, larger contribution rates and more edges. That a subnet is more important implied that the subnet contained more information on the SNP epistatic effect affecting AFW.

SNP-SNP Network Annotation
First, the significant epistatic interactive SNPs were mapped to QTLs affecting chicken AFW (Table S2) Second, fragments of 0.4 Mb were defined, and the fragments' centers were the 68 single SNPs in the significant epistatic interactive SNP pairs. The distances between some SNPs were small, and the fragments may overlap. A union set was used when the fragments overlapped. Regions were determined according to the fragments and union sets ( Table 4). The 97 genes (including nine non-coding genes) in these regions are listed in Table 4. Information on protein coding genes is listed in Table S3.
Third, genes in every subnet were submitted to the KEGG database. The 97 genes were located in 50 pathways (Table S4). Some genes were located in the same pathway, for example GRB2, PDPK1, PIK3CA and SOCS3 are in the insulin-signaling pathway; and GRB2, PIK3CA, SOCS3 and IL21R are in the Jak-STATsignaling pathway.

Discussion
Abdominal fatness traits in chickens are important quantitative traits, and are under complex genetic control. Understanding the molecular mechanisms underlying them will help efficient growth selection in broiler chickens. In this study, analyses using genomewide SNPs were performed to examine the genetic contribution of epistatic effects to the phenotypic variation of chicken abdominal fatness traits.
GWASs have identified many significant association markers. However, many studies showed that the identified markers explain only a fraction of the heritable variance [11][12][13]. This discrepancy has been ascribed to the insufficient power of GWAS to identify gene-gene interactions [13]. Zuk et al. [29] showed that the ''missing'' heritability of Crohn's disease is 62.8%, and that genetic interactions could account for 80% of this missing heritability. Recently, most genome-wide association analyses in chickens have focused on individual SNP marker association analysis [20,30]. However, quantitative traits often arise from the combined effects  of multiple loci [14][15][16]. In the present study, the contribution rates of the significant epistatic interactive SNP pairs to AFW ranged from 0.62% to 1.54%. Two SNP pairs showed the largest contribution rate of 1.54%, and the contribution rates of 47 pairs were more than 1%. In addition, we calculated the contribution rates of single SNPs in the GWAS (Table S5). The single SNP Gga_rs14276105, which was significant in GWAS (Table S5), showed the largest contribution rate 1.46%. Comparing the results showed that the contribution rates of significant SNP pairs was larger than that of the single SNPs contained in the pairs. Although the contribution rates of the significant epistatic interactive SNP pairs might be overestimated for model (1), there were no high dimension interactions; therefore, we concluded that epistasis has an important effect on the phenotypic variance and  cannot be ignored. Our results also showed that Gga_rs14276105 was the only SNP that was significant in an association study among the significantly interacting SNPs (Table S5). Wu et al. [31] also demonstrated that the majority of significantly interacting SNPs did not show marginal association in their data. This suggested that testing interactions of SNPs with only one significant locus association would lose the majority of interactions. The features of scale-free networks are several nodes with large degrees and many nodes with few connections. Further, the nodes with large degrees should be the hub sites of the network. Our SNP-SNP network displayed the features of a scale-free network. The degree of one node was 7, and of another was 4; and 47 nodes had a degree of 1 in our SNP-SNP network. Thus, the nodes (Gga_rs14988623 and Gga_rs14303341) with larger degrees should be the important ones. Scale-free behavior has been demonstrated in many biological interaction networks, such as protein or gene interaction networks [32].
The formula for evaluating the importance of subnets proposed in the present study was based on the entropy theory [26]. The evaluation indices comprised the edge number, P-value and the contribution rate of every SNP pair in the subnets. The importance evaluation score stands for the information quantity of a subnet. A subnet containing more information would be more important. We suggested that subnets A and C are the most important; as they received the top two importance evaluation scores (Table 3).
Subnet A was radial, with seven SNPs connecting the center. A previous study concluded that many significant epistatic effects involving one locus were less likely to be random than was a single epistatic effect [33]; i.e., the network is less likely to be random. Combined with the results of the scale-free network analysis, we identified Gga_rs14303341 as the important SNP. The annotated information indicated that Gga_rs14303341 maps to QTLs (11809, 11817) [34], and the genes GOSR2, GJC1, CCDC43, WNT3, NSF and EFTUD2 are within the defined 0.4 Mb region. These QTLs and genes should receive attention in a future study.
Subnet C was a loop graph, and contained six SNPs: three on GGAZ and three on GGA8. Between the two chromosomes, nine (363) pairs of SNP interactions were obtained, and eight of these interactions were significant. In addition, the SNPs were adjacent to one another on GGAZ and GGA8. The result indicated that the interaction occurs between the two regions on GGAZ (59,561,642-60,090,315) and GGA8 (28,415,925,822).
The epistatic interactions of SNPs might provide the basis for identifying interactions among genes, QTLs, or even pathways; therefore, we annotated the network. The result showed that 35% (24/68) of the SNP markers are located in regions containing 22 QTLs associated with chicken AFW (Table S2). The radial network based on subnet A suggested interactions between QTLs. A QTLs network with a similar topological structure was found in a study of chicken growth traits by Carlborg et al. [16]. They identified five significant epistatic QTL pairs that also formed a radial network and that accounted for a significant amount of the phenotypic variance.
Gene-by-gene interaction plays a major role in genetic studies of quantitative traits; however, the detection of gene-gene interaction has been traditionally assessed by SNP interactions. Single SNPs cannot capture the total variation of a gene; therefore, SNP interactions cannot represent gene-gene interactions. Thus, Cui et al. [35] and Li et al. [36] extended the idea of single SNP interactions to haplotype interactions, and proposed a novel statistical approach for capturing variations in genes and potential interactions between two genes. In the present study, we found an interesting phenomenon in some subnets, in which SNPs on the circumference were very near in physical distance. Therefore, we believe that the definition of block interaction may provide additional biological insights, compared with the analysis of SNP interactions.
In this study, we detected 26 regions on 12 chromosomes. Ninety-seven genes (including nine non-coding genes), including BMPR1A, GIPR, GRB2, LITAF, SOCS3, WNT3 and PDPK1, were identified in these regions. Based on the literature, some genes are associated with obesity; for example, BMPR1A is associated with human obesity [37]. GIPR on GGA18 is associated with fat droplet formation [38]. LITAF plays an important role in the regulation of tumor necrosis factor alpha (TNF-a), which is important in adipose tissue [39,40]. WNT3 is part of the Wnt-signaling pathway, which is one of the pathways operating in fat cells [41]. PDPK1 affects the fat-pad mass in mice [42]. We identified 50 pathways related to AFW, according to gene annotation analysis (Table S4). Some of the pathways are related to obesity, such as the Jak-STAT-signaling pathway and the insulin-signaling pathway [43,44]. their molecular mechanisms, and will contribute to faster progress in the artificial genetic selection of broiler chickens.

Supporting Information
Table S1 SNPs contained in the significant SNP pairs. (DOC)