Skip to main content


Identification of genome-wide SNP-SNP interactions associated with important traits in chicken

Article metrics



In addition to additive genetic effects, epistatic interactions can play key roles in the control of phenotypic variation of traits of interest. In the current study, 475 male birds from lean and fat chicken lines were utilized as a resource population to detect significant epistatic effects associated with growth and carcass traits.


A total of 421 significant epistatic effects were associated with testis weight (TeW), from which 11 sub-networks (Sub-network1 to Sub-network11) were constructed. In Sub-network1, which was the biggest network, there was an interaction between GGA21 and GGAZ. Three genes on GGA21 (SDHB, PARK7 and VAMP3) and nine genes (AGTPBP1, CAMK4, CDC14B, FANCC, FBP1, GNAQ, PTCH1, ROR2 and STARD4) on GGAZ that might be potentially important candidate genes for testis growth and development were detected based on the annotated gene function. In Sub-network2, there was a SNP on GGA19 that interacted with 8 SNPs located on GGA10. The SNP (Gga_rs15834332) on GGA19 was located between C-C motif chemokine ligand 5 (CCL5) and MIR142. There were 32 Refgenes on GGA10, including TCF12 which is predicted to be a target gene of miR-142-5p. We hypothesize that miR-142-5p and TCF12 may interact with one another to regulate testis growth and development. Two genes (CDH12 and WNT8A) in the same cadherin signaling pathway were implicated as potentially important genes in the control of metatarsus circumference (MeC). There were no significant epistatic effects identified for the other carcass and growth traits, e.g. heart weight (HW), liver weight (LW), spleen weight (SW), muscular and glandular stomach weight (MGSW), carcass weight (CW), body weight (BW1, BW3, BW5, BW7), chest width (ChWi), metatarsus length (MeL).


The results of the current study are helpful to better understand the genetic basis of carcass and growth traits, especially for testis growth and development in broilers.


Epistasis arises due to interactions, either between single nucleotide polymorphisms (SNPs), genes or quantitative trait loci (QTLs), which result in non-linear effects that control variation in phenotypes. Epistasis can have a large influence on phenotypic variation of traits such as starvation resistance, startle response, and chill coma recovery in Drosophila [1]. Identification of epistatic effects associated with quantitative traits will help us to better understand the genetic architecture that underlies complex variation of phenotypes for both humans and animals [1, 2]. Therefore, more and more attention has been placed on epistasis, which has resulted in some valuable insights [3,4,5,6,7,8].

With the advent of SNP arrays and genomic re-sequencing, it is relatively easy to genotype a wide array of individuals. As a result, many genome–wide association studies (GWAS) have been carried out in the past several years. Most of these studies have focused on single locus additive genetic tests. However, this is not that only type of genetic association. A genome wide SNP-SNP interaction analysis should provide new insights into the genetic architecture that underlies variation in complex traits.

Recently, genome wide SNP-SNP interaction analysis have been conducted in humans [4, 9] and domestic livestock species [10,11,12]. In chicken (Gallus gallus), it has been suggested that epistatic interactions between genes (or QTLs) are important for variation in quantitative traits [8, 13, 14]. However, the study of epistatic interactions at the whole genome level have been limited [15].

Two Northeast Agricultural University broiler lines that have been divergently selected for abdominal fat content (named as NEAUHLF) for more than 10 years were used in the current study. Previously, we reported that 52 pairs of SNPs had significant epistatic interactions that were associated with abdominal fat weight [15]. In the current study, the significant epistatic interactions for carcass and growth traits were identified. The results of this study may provide some helpful information to better understand the genetic basis of carcass and growth traits in broilers.


Experimental populations

Two Northeast Agricultural University broiler lines that have been divergently selected for abdominal fat content (NEAUHLF) were used to identify epistatic interactions. The NEAUHLF lines have been selected since 1996 using abdominal fat percentage (AFP = abdominal fat weight/body weight at 7 weeks of age) and plasma very low-density lipoprotein (VLDL) concentration as selection criteria. The G0 generation of NEAUHLF came from the same grandsire line, which originated from the Arbor Acres broiler, which was then divided into two lines according to their plasma VLDL concentration at 7 weeks of age. The G0 birds were mated (one sire: four dams) to produce 25 half-sib families for each line, with an average of 70 G1 offspring per family in two hatches. From G1 to G11, the birds of each line were raised in two hatches with five birds per cage. Plasma VLDL concentrations were measured for all male birds, which had free access to feed and water at 7 weeks, and the AFP of the male birds in the first hatch was measured after slaughter at 7 weeks. Sib birds from the families with lower (lean line) or higher (fat line) AFP than the average value for the population were selected as candidates for breeding, considering the plasma VLDL concentration and the body weights of male birds in the second hatch and the egg production of female birds in both hatches. These birds were kept under the same environmental conditions and had free access to feed and water. Commercial corn-soybean-based diets that met all National Research Council (NRC) requirements were provided. From hatch to 3 weeks of age, the birds received a starter feed (3,000kal ME = kg and 210 g = kg CP) and from 4 weeks of age to slaughter the birds were fed a grower diet (3100 kal ME = kg and 190 g = kg CP). The birds used in the current study included 475 male individuals from the 11th generation of NEAUHLF [15]. The birds were weighed at 0, 1, 3, 5 and 7 weeks of age (BW0, BW1, BW3, BW5 and BW7). At 7 weeks of age, the metatarsus length (MeL), metatarsus circumference (MeC) and chest width (ChWi) were measured prior to slaughter as described previously [16]. Carcass weight (CW), testis weight (TeW), heart weight (HW), liver weight (LW), spleen weight (SW), muscular and glandular stomach weight (MGSW) were obtained after the birds were slaughtered.

SNP genotyping

Genotyping was carried out using the Illumina Inc. (San Diego, CA, USA) chicken 60 K SNP chip, which contained 57,636 SNPs. After quality control, 48,824 SNPs in 475 individuals were used in the epistatic interaction analyses. The quality control of the SNP genotypes was described previously by Zhang et al. [17].

Genome-wide Pairwise interaction analysis

The EPISNP3 module in epiSNP_v4.2_Windows software package was used to identify significant epistatic effects [18]. The statistical model used to test for epistatic effects associated with carcass and growth traits was as follows: y = Xg + Zb + e, where y is the column vector of phenotypic values, g is the effects of SNP genotypes, X is the design matrix of g, b is the fixed effects of Line and BW7 (or BW0), Z is the model matrix of b, and e is the random error.

The P-values of the epistatic effects were Bonferroni corrected for multiple testing (5.96 × 109 independent tests, with a significance threshold of P < 0.05), which resulted in P < 8.39 × 10−12 as a significance threshold. Significant interactions, including additive by additive (AA), additive by dominance (AD) or dominance by additive (DA) and dominance by dominance (DD), between two SNPs on the same chromosome were deleted because these interactions may potentially be markers for a haplotype effect that contains a single QTL [12]. The remaining significant SNP interactions were further filtered with the criterion that only those interaction with at least 10 animals in every genotype combination were considered [12], which is roughly equivalent to a 15% minor allele frequency for each variant in an additive by additive epistatic interaction.

SNP-SNP network

The figures that illustrate the SNP-SNP networks with the significant epistatic effects for carcass and growth traits were drawn using the epiNet option within the epiSNP_v4.2_Windows software package [18].

Linkage disequilibrium (LD) analysis

The linkage disequilibrium (LD) between SNPs was calculated using Haploview software (version 4.2). The solid spine method within the package was used to define the LD block.

Annotation of SNP-SNP network

Genes within 1 Mb (upstream and downstream) of the SNPs that had significant interactions with another SNP for carcass and growth traits were retrieved from UCSC ( (Galgal4). Functional annotation of genes was performed using DAVID bioinformatics resources 6.8 ( for Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Statistical significance was set at the nominal P-value < 0.05.


Phenotypic and SNPs information

Phenotypic summary statistics for carcass and growth traits are shown in Table 1 and the phenotypic distributions for the traits in the lean and fat lines respectively are shown in Additional file 1: Figure S1. There were extremely significant differences (P < 0.01) in Chwi, MeL, MeC, LW, SW, TeW and significant difference (P < 0.05) in HW between the lean and fat lines. After quality control, 48,824 SNPs were utilized for epistatic interaction analyses (Table 2). These SNPs were distributed on 28 autosomes, Z chromosome, two linkage groups, and SNPs not assigned to any chromosomes in chickens. These markers covered about 1027.01 Mb of the chicken genome, with an average SNP density of 16.08 kb/SNP.

Table 1 The Mean ± Standard deviation (SD) of the carcass and growth traits in lean and fat lines, respectively, and in the combined population
Table 2 Summary information of the genome-wide SNP markers

Epistatic analysis of carcass trait

The pairwise interaction effects between every two SNPs across the whole chicken genome for TeW were calculated using EPISNP3 [18]. After filtering, 421 pairs of SNPs were significantly associated with testis weight (P < 8.39 × 10−12). Of these 421 significant SNP by SNP interactions, 403 (95.72%) exhibited an additive by additive interaction, 18 (4.28%) exhibited an additive by dominance (or dominance by additive) interaction, and no dominance by dominance interactions were detected (Additional file 2: Table S1). The most significant additive by additive effect detected occurred between GGA3 (GGaluGA22768) and GGA10 (Gga_rs14722408). The phenotypic distributions of the four different genotype classes of this additive by additive effect are showed in Additional file 3: Figure S2.

To investigate the complex mechanism of epistatic effects on TeW, networks were constructed using the 421 significant SNP by SNP interactions. The epistatic interaction sub-networks that contained more than three nodes are shown in Additional file 4: Figure S3. Eleven sub-networks were detected (Additional file 4: Figure S3). Sub-network1 was the biggest and contained 372 pairs of SNP by SNP interaction effects. Based on LD information, a simpler sub-network, which was derived from Sub-network1, was obtained, in which numerous SNPs are represented by a single LD block (Fig. 1, Additional file 5: Figure S4). The blocks in the simpler sub-network represented the SNPs that contained in the LD blocks. Two hundred and fifty-five of the interactions in Sub-network1 were between GGA21 and GGAZ, which indicated an interaction between the two chromosomes. The 255 interactions detected between GGA21 and GGAZ involved 24 SNPs on GGA21 and 19 SNPs on GGAZ, which spanned 572 kb (from 76,023 bp to 647,587 bp) and 9.5 Mb (from 37,246,321 bp to 46,745,968 bp), respectively. There were 13 Refgenes in the 572 kb region on GGA21 and 41 Refgenes in the 9.5 Mb region on GGAZ (Table 3). Six GO terms, including protein-arginine deiminase activity, protein citrullination, positive regulation of collateral sprouting, cytoplasm, cell fate determination and protein autophosphorylation, were significantly (P < 0.05) enriched. No significant KEGG pathways were detected. Three genes on GGA21 (SDHB, PARK7 and VAMP3) and nine genes (AGTPBP1, CAMK4, CDC14B, FANCC, FBP1, GNAQ, PTCH1, ROR2 and STARD4) on GGAZ might be important for testis growth and development based on their basic functions.

Fig. 1

A simplified network of Sub-network1 using the LD information. The blocks are from the LD information in Additional file 3: Figure S2 and the lines between two nodes depicts an interaction between these two nodes. The colors of the nodes represent the P-value of an interaction (P < 1.0 × 10−16 = red; P < 1.0 × 10−15 = blue; P < 1.0 × 10−14 = green; P < 1.0 × 10−13 = white). The color of the edge indicates the type of epistatic effect (AA = red; AD = purple; DA = blue; DD = green)

Table 3 The SNPs on GGAZ interacted with the SNPs on GGA21, and the Refgenes in the two important regions

Sub-networks 2, 3, 4, 5, 8 and 11 each had several SNPs in the same LD block on one chromosome, which interacted with a single SNP on another chromosome. For example, in Sub-network2 eight SNPs on GGA10 interacted with the SNP (Gga_rs15834332) on GGA19 (Fig. 2). The eight SNPs on GGA10 were spread across a 4.3 Mb region that contained 32 chicken Refgenes (Fig. 2). The SNP Gga_rs15834332 on GGA19 was located between two genes, C-C motif chemokine ligand 5 (CCL5) and MIR142, which is related to two miRNAs, miR-142-5p and miR-142-3p. Interestingly, five genes in the 4.3 Mb region on GGA10 are predicted to be target genes of miR-142-5p and miR-142-3p (Table 4). Among these five target genes, transcription factor 12 (TCF12) was the only gene that was predicted to be the target gene of miR-142-5p by three different packages, including TargetScan (, miRDB ( and PicTar ( Unfortunately, it was impossible to predict which genes in the other regions of the genome may be good candidates to control testis growth and development.

Fig. 2

Sub-network2 for testis weight (TeW) and the LD information. The color of the node represents the P-value of an interaction (P < 1.0 × 10−16 = red; P < 1.0 × 10−15 = blue; P < 1.0 × 10−14 = green; P < 1.0 × 10−13 = white). The color of the edge indicates the type of epistatic effect (AA = red; AD = purple; DA = blue; DD = green). The genes located in the 4.3 Mb regions of GGA10 were listed

Table 4 Target genes of miR-142-5p and miR-142-3p in the 4.3 Mb region on GGA10 in Sub-network2 for TeW predicted by three packages online

For CW, HW, LW, SW and MGSW, no significant epistatic interactions were detected.

Epistatic analysis of growth trait

For BW1, BW3, BW5, BW7, ChWi and MeL, no significant epistatic interactions were detected. Fifteen pairs of SNPs with significant interaction effects on MeC were detected (Table 5). These 15 interactions were all additive by additive interactions, which implicated an interaction between GGA2 and GGA13 (Table 5). There was a single network that contained all 15 interactions (Fig. 3). The fifteen interactions occurred between five SNPs in a single LD block on GGA2, and three SNPs in a single LD block on GGA13 (Fig. 3). The genes inside the two LD blocks and within 0.5 Mb 5′ and 3′ of the LD blocks were found. There was only a single Refgene (Cadherin-12, CDH12) in the region on GGA2. There were fifteen Refgenes (CDX1, CSF1R, NPY7R, FLT4, CANX, HNRNPH1, DGUOK, RUFY1, MAPK9, RASGEF1C, HNRNPAB, NME5, WNT8A, FAM13B and NPY6R) located in the region on GGA13.

Table 5 Significant epistatic effects on MeC
Fig. 3

Epistatic network among SNPs with significant epistatic effect on metatarsus circumference (MeC). A node represents a SNP. The chromosome in which the SNP is located is shown in the circle. A pair of SNPs connected by an edge had a significant interaction. The color of a node represent the P-value of the interaction (P < 1.0 × 10−12 = red; P < 1.0 × 10−11 = blue; P < 1.0 × 10−10 = green; P < 1.0 × 10−9 = white). The color of the edge indicates the type of epistatic effect (AA = red; AD = purple; DA = blue; DD = green). Genes in the LD blocks were listed


The birds from lean and fat chicken lines used in the current study had significantly different amounts of abdominal fat content and significantly different testis weight, which was an ideal population to study the genetic architecture of abdominal fat deposition and testis growth and development [19]. The birds used in the current study were 7-week-old, therefore, the results of the current study could reflect early testis development. Testis weight has been reported to be controlled by genetic factors in both chickens and mice [20,21,22,23,24]. Roosters with small testes often have poor fertility [25]. There have been previously reported studies on the genetics of testis development in chicken [19]. However, most of the previous studies focused on additive genetic effects, while, epistatic effects were ignored. It is unknown how important epistatic interactions are for development of the testis. In the current study, we used 475 birds from lean and fat lines to conduct epistasis analysis for testis weight and other carcass and growth traits in chickens.

The two lines used in the current study were selected for 11 years and some regions of the genome may be fixed because of selection pressure. If this occurred, it would not be possible to identify epistatic interactions in these regions of the genome. In contrast, if the two lines were crossed, it would be possible to detect epistatic effects in this intercross population. Despite this problem, some epistatic effects were detected in these two lines using epiSNP. Due to the selection applied to these population, it is possible that genomic stratification has occurred. epiSNP, which was used to conduct these analyses is capable of adjusting for family structure [18].

In this study, the SNP by SNP interaction effects for carcass and growth traits were filtered by three criteria. First, to correct for multiple testing comparisons, a Bonferroni correction was used instead of false discovering rate (FDR) method to minimize false positives. Utilization of the Bonferroni correction method will however also decrease the discovery of true interactions compared to the use of a false discovery rate correction. Therefore, the SNP X SNP interaction detected in this study may be fraction of the interactions that affect these. Second, when considering the linkage disequilibrium (LD) between the SNPs and the QTL for traits of interest, an interaction between two SNPs on the same chromosome may detect a haplotype effect but not an interaction. In other words, it is not possible to separate true interaction effects from haplotype effects if the two SNPs are located on the same chromosome. Therefore, in order to increase the power to detect the true interaction effects for the growth and carcass traits, we deleted interactions that occurred on the same chromosome [12]. Third, the smaller the number of birds in any given genotype class, the less likely it is to get a good estimate of the genotype effect. Thus, in order to increase the power to detect the true interactions, only those interactions that contained at least 10 animals in every genotype combination were considered [12]. We carried out the filter of the interactions according to these criteria in order to reduce the chances of obtaining false-positive results (type I errors).

For testis weight, a total of 421 pairs of SNP-SNP epistatic interactions were detected. These pairs of SNPs comprised 211 single SNPs, and none of these individual SNPs were identified in our previous GWAS analysis for testis weight [19]. A similar phenomenon was also detected by Wu et al. for psoriasis in human [26]. These results indicated that all SNPs on the chip should be tested for identifying potential interaction effects. In contrast, testing for interactions between SNPs that have been previously identified in GWAS is not enough.

There were 11 sub-networks with at least three nodes that were identify by using the 421 pairs of interaction effects. Sub-network1 was the biggest one and most of the interactions in Sub-network1 occurred between GGAZ and GGA21. The two regions on GGA21 and GGAZ that may harbor genes important for testis growth and development spanned 572 kb and 9.5 Mb, respectively. Three genes on GGA21 (SDHB, PARK7 and VAMP3) and nine genes (AGTPBP1, CAMK4, CDC14B, FANCC, FBP1, GNAQ, PTCH1, ROR2 and STARD4) on GGAZ might be important for testis growth and development based on their annotated functions. The previous result indicated that the motility and viability of sperm were positively correlated with mitochondrial SDHB [27]. Therefore, SDHB may serve as a marker of sperm quality and male fertility [28]. PARK7 (DJ1) is highly expressed in human testes and has been shown to be essential for sperm maturation and fertilization [29,30,31,32]. VAMP3 has been shown to play an important role in the process of fertilization of sperms in pig [33]. In mice, AGTPBP1 was important for spermatogenesis, moreover, it was important for the survival of germ cells from the spermatocyte stage onward [34]. In mice and rats, the the CAMK4 gene encodes two proteins, Ca2+/calmodulin-dependent protein kinase IV (CaMKIV) and calspermin (CaS) [35,36,37,38]. CaMKIV is highly expressed in mouse testis and ovary and plays a essential role in male and female fertility [38,39,40]. CDC14B mutant mice were less fertile than the wild-type control [41]. Reduced fertility was reported for Fancc−/− mice [42]. The expression of some proteins, including FBP1, were altered and their functions may be damaged in infertile men with unilateral varicocele [43]. The results of a previous study identified that Gnaqd/d male mice were subfertile [44]. The desert hedgehog (Dhh)-null mutant male mice had less mature sperm cells and lower numbers of Leydig cells (LCs), and Dhh played an important role in spermatogenesis by acting in a paracrine manner through the Ptch1 receptor component [45,46,47]. In mice, female Ror2W749FLAG/W749FLAG were fertile, however, Ror2W749FLAG/W749FLAG male mice showed a decreased in fertility [48]. StarD6 was testis-specific expressed which indicated that it may be important for fertility [49]. StarD6 is homology to StarD4, which indicated that StarD4 may have similar function as StarD6. In Sub-network2, eight SNPs on GGA10 all interacted with the SNP (Gga_rs15834332) on GGA19, which could be seen as the hub site of Sub-network2. Thus, Gga_rs15834332 may be the important node which interacts with a 4.3 Mb region on GGA10. The Gga_rs15834332 on GGA19 was located between CCL5 and MIR142. A total of 32 Refgenes were located in the 4.3 Mb region of GGA10. Among these genes, TCF12 was predicted as a target gene for miR-142-5p [50]. We also detected that TCF12 was the only gene that was predicted to be the target of miR-142-5p using three packages online. Therefore, it is proposed that miR-142-5p and TCF12 might work together to regulate the reproductive function of male broilers. Furthermore, TCF12 was a partner of TCF21, which was detected as an important gene for testis growth and development in our previous GWAS result [19].

For MeC, 15 pairs of SNPs with significant epistatic effects were detected, which indicated an interaction between GGA2 and GGA13. It is proposed that CDH12 on GGA2 and WNT8A on GGA13, which are both located in the same cadherin signaling pathway, may be important for bone growth. It had been shown that the pathway was involved in many biological processes, such as development, neurogenesis, cell adhesion, and inflammation, and also involved in many disease, such as cancer [51].


In the current study, a large number of epistatic interactions were found to be significantly associated with testis weight in chicken. It appears that miR-142-5p along with its target gene TCF12, and some other genes in GGA21 and GGAZ (SDHB, PARK7, VAMP3, AGTPBP1, CAMK4, CDC14B, FANCC, FBP1, GNAQ, PTCH1, ROR2 and STARD4) might be important for testis growth and development. In contrast, very few significant epistatic interactions were identified for other carcass and growth traits. These results indicate that epistatic interaction may play very different roles in the control of phenotypic variation for different traits in chickens.



Gene ontology


Genome–wide association study


Kyoto encyclopedia of genes and genomes


Linkage disequilibrium


Quantitative trait loci


Single nucleotide polymorphism


  1. 1.

    Huang W, Richards S, Carbone MA, Zhu D, Anholt RR, Ayroles JF, et al. Epistasis dominates the genetic architecture of drosophila quantitative traits. Proc Natl Acad Sci U S A. 2012;109(39):15553–9.

  2. 2.

    Rajon E, Plotkin JB. The evolution of genetic architectures underlying quantitative traits. Proc Biol Sci. 2013;280(1769):20131552.

  3. 3.

    Hibar DP, Stein JL, Jahanshad N, Kohannim O, Toga AW, McMahon KL, et al. Exhaustive search of the SNP-sNP interactome identifies epistatic effects on brain volume in two cohorts. Med Image Comput Comput Assist Interv. 2013;16(Pt 3):600–7.

  4. 4.

    Hibar DP, Stein JL, Jahanshad N, Kohannim O, Hua X, Toga AW, et al. Genome-wide interaction analysis reveals replicated epistatic effects on brain structure. Neurobiol Aging. 2015;36(Suppl 1):151–8.

  5. 5.

    Young KL, Graff M, North KE, Richardson AS, Bradfield JP, Grant SF, et al. Influence of SNP*SNP interaction on BMI in European American adolescents: findings from the National Longitudinal Study of adolescent health. Pediatr Obes. 2016;11(2):95–101.

  6. 6.

    Carlborg O, Kerje S, Schütz K, Jacobsson L, Jensen P, Andersson LA. Global search reveals epistatic interaction between QTL for early growth in the chicken. Genome Res. 2003;13(3):413–21.

  7. 7.

    Cheng HH, Zhang Y, Muir WM. Evidence for widespread epistatic interactions influencing Marek's disease virus viremia levels in chicken. Cytogenet Genome Res. 2007;117(1-4):313–8.

  8. 8.

    Carlborg O, Jacobsson L, Ahgren P, Siegel P, Andersson L. Epistasis and the release of genetic variation during long-term selection. Nat Genet. 2006;38(4):418–20.

  9. 9.

    Li J, Zhang Q, Chen F, Yan J, Kim S, Wang L, et al. Genetic interactions explain variance in Cingulate Amyloid burden: an AV-45 PET genome-wide association and interaction study in the ADNI cohort. Biomed Res Int. 2015;2015:647389.

  10. 10.

    Kogelman LJ, Kadarmideen HN. Weighted Interaction SNP. Hub (WISH) network method for building genetic networks for complex diseases and traits using whole genome genotype data. BMC Syst Biol. 2014;8(Suppl 2):5.

  11. 11.

    Ali AA, Khatkar MS, Kadarmideen HN, Thomson PC. Additive and epistatic genome-wide association for growth and ultrasound scan measures of carcass-related traits in Brahman cattle. J Anim Breed Genet. 2015;132(2):187–97.

  12. 12.

    Kramer LM, Ghaffar MA, Koltes JE, Fritz-Waters ER, Mayes MS, Sewell AD, et al. Epistatic interactions associated with fatty acid concentrations of beef from angus sired beef cattle. BMC Genomics. 2016;17(1):891.

  13. 13.

    Hu G, Wang SZ, Wang ZP, Li YM, Li H. Genetic epistasis analysis of 10 peroxisome proliferator-activated receptor γ-correlated genes in broiler lines divergently selected for abdominal fat content. Poult Sci. 2010;89(11):2341–50.

  14. 14.

    Ek W, Marklund S, Ragavendran A, Siegel P, Muir W, Carlborg O. Generation of a multi-locus chicken introgression line to study the effects of genetic interactions on metabolic phenotypes in chickens. Front Genet. 2012;3:29.

  15. 15.

    Li F, Hu G, Zhang H, Wang S, Wang Z, Li H. Epistatic effects on abdominal fat content in chickens: results from a genome-wide SNP-SNP interaction analysis. PLoS One. 2013;8(12):e81520.

  16. 16.

    Zhang H, Zhang YD, Wang SZ, Liu XF, Zhang Q, Tang ZQ, et al. Detection and fine mapping of quantitative trait loci for bone traits on chicken chromosome one. J Anim Breed Genet. 2010;127(6):462–8.

  17. 17.

    Zhang H, Hu X, Wang Z, Zhang Y, Wang S, Wang N, et al. Selection signature analysis implicates the PC1/PCSK1 region for chicken abdominal fat content. PLoS One. 2012;7(7):e40736.

  18. 18.

    Ma L, Runesha HB, Dvorkin D, Garbe JR, Da Y. Parallel and serial computing tools for testing single-locus and epistatic SNP effects of quantitative traits in genome-wide association studies. BMC Bioinformatics. 2008;9:315.

  19. 19.

    Zhang H, Na W, Zhang HL, Wang N, ZQ D, Wang SZ, et al. TCF21 is related to testis growth and development in broiler chickens. Genet Sel Evol. 2017;49(1):25.

  20. 20.

    Sarabia Fragoso J, Pizarro Díaz M, Abad Moreno JC, Casanovas Infesta P, Rodriguez-Bertos A, Barger K. Relationships between fertility and some parameters in male broiler breeders (body and testicular weight, histology and immunohistochemistry of testes, spermatogenesis and hormonal levels). Reprod Domest Anim. 2013;48(2):345–52.

  21. 21.

    Chubb C. Genes regulating testis size. Biol Reprod. 1992;47(1):29–36.

  22. 22.

    Lüpold S, Linz GM, Rivers JW, Westneat DF, Birkhead TR. Sperm competition selects beyond relative testes size in birds. Evolution. 2009;63(2):391–402.

  23. 23.

    Soulsbury CD. Genetic patterns of paternity and testes size in mammals. PLoS One. 2010;5(3):e9581.

  24. 24.

    Vizcarra JA, Kirby JD, Kreider DL. Testis development and gonadotropin secretion in broiler breeder males. Poult Sci. 2010;89(2):328–34.

  25. 25.

    John P. Testes development and fertility. Aviagen brief. 2008;0608-AVN-011.

  26. 26.

    Wu X, Dong H, Luo L, Zhu Y, Peng G, Reveille JD, et al. A novel statistic for genome-wide interaction analysis. PLoS Genet. 2010;6(9):e1001131.

  27. 27.

    Xue XP, Shang XJ, Fu J, Chen YG, Shi YC. Detection and significance of succinate dehydrogenase of sperm mitochondria. Zhonghua Nan Ke Xue. 2003;9(8):601–3.

  28. 28.

    Rahman MS, Kwon WS, Lee JS, Yoon SJ, Ryu BY, Pang MG. Bisphenol-a affects male fertility via fertility-related proteins in spermatozoa. Sci Rep. 2015;5:9169.

  29. 29.

    Nagakubo D, Taira T, Kitaura H, Ikeda M, Tamai K, Iguchi-Ariga SM, et al. DJ-1, a novel oncogene which transforms mouse NIH3T3 cells in cooperation with ras. Biochem Biophys Res Commun. 1997;231(2):509–13.

  30. 30.

    Klinefelter GR, Laskey JW, Ferrell J, Suarez JD, Roberts NL. Discriminant analysis indicates a single sperm protein (SP22) is predictive of fertility following exposure to epididymal toxicants. J Androl. 1997;18(2):139–50.

  31. 31.

    Ooe H, Taira T, Iguchi-Ariga SM, Ariga H. Induction of reactive oxygen species by bisphenol a and abrogation of bisphenol A-induced cell injury by DJ-1. Toxicol Sci. 2005;88(1):114–26.

  32. 32.

    Hao LY, Giasson BI, Bonini NM. DJ-1 is critical for mitochondrial function and rescues PINK1 loss of function. Proc Natl Acad Sci U S A. 2010;107(21):9747–52.

  33. 33.

    Tsai PS, Garcia-Gil N, van Haeften T, Gadella BM. How pig sperm prepares to fertilize: stable acrosome docking to the plasma membrane. PLoS One. 2010;5(6):e11204.

  34. 34.

    Kim N, Xiao R, Choi H, Jo H, Kim JH, et al. Abnormal sperm development in pcd(3J)−/− mice: the importance of Agtpbp1 in spermatogenesis. Mol Cells. 2011;31(1):39–48.

  35. 35.

    Ono T, Means AR. Calspermin is a testis specific calmodulin-binding protein closely related to Ca2+/calmodulin-dependent protein kinases. Adv Exp Med Biol. 1989;255:263–8.

  36. 36.

    Ono T, Slaughter GR, Cook RG, Means AR. Molecular cloning sequence and distribution of rat calspermin, a high affinity calmodulin-binding protein. J Biol Chem. 1989;264(4):2081–7.

  37. 37.

    Ohmstede CA, Bland MM, Merrill BM, Sahyoun N. Relationship of genes encoding Ca2+/calmodulin-dependent protein kinase Gr and calspermin: a gene within a gene. Proc Natl Acad Sci U S A. 1991;88(13):5784–8.

  38. 38.

    JY W, Ribar TJ, Cummings DE, Burton KA, McKnight GS, Means AR. Spermiogenesis and exchange of basic nuclear proteins are impaired in male germ cells lacking Camk4. Nat Genet. 2000;25(4):448–52.

  39. 39.

    JY W, Gonzalez-Robayna IJ, Richards JS, Means AR. Female fertility is reduced in mice lacking Ca2+/calmodulin-dependent protein kinase IV. Endocrinology. 2000;141(12):4777–83.

  40. 40.

    JY W, Ribar TJ, Means AR. Spermatogenesis and the regulation of ca(2+)-calmodulin-dependent protein kinase IV localization are not dependent on calspermin. Mol Cell Biol. 2001;21(17):6066–70.

  41. 41.

    Wei Z, Peddibhotla S, Lin H, Fang X, Li M, Rosen JM, et al. Early-onset aging and defective DNA damage response in Cdc14b-deficient mice. Mol Cell Biol. 2011;31(7):1470–7.

  42. 42.

    Whitney MA, Royle G, Low MJ, Kelly MA, Axthelm MK, Reifsteck C, et al. Germ cell defects and hematopoietic hypersensitivity to gamma-interferon in mice with a targeted disruption of the Fanconi anemia C gene. Blood. 1996;88(1):49–58.

  43. 43.

    Agarwal A, Sharma R, Durairajanayagam D, Ayaz A, Cui Z, Willard B, et al. Major protein alterations in spermatozoa from infertile men with unilateral varicocele. Reprod Biol Endocrinol. 2015;13:8.

  44. 44.

    Babwah AV, Navarro VM, Ahow M, Pampillo M, Nash C, Fayazi M, et al. GnRH neuron-specific ablation of Gαq/11 results in only partial inactivation of the Neuroendocrine-reproductive Axis in both male and female mice: in vivo evidence for Kiss1r-coupled Gαq/11-independent GnRH secretion. J Neurosci. 2015;35(37):12903–16.

  45. 45.

    Bitgood MJ, Shen L, McMahon AP. Sertoli cell signaling by desert hedgehog regulates the male germline. Curr Biol. 1996;6(3):298–304.

  46. 46.

    Yao HH, Whoriskey W, Capel B. Desert hedgehog/patched 1 signaling specifies fetal Leydig cell fate in testis organogenesis. Genes Dev. 2002;16(11):1433–40.

  47. 47.

    Morales CR, Fox A, El-Alfy M, Ni X, Argraves WS. Expression of Patched-1 and smoothened in testicular meiotic and post-meiotic cells. Microsc Res Tech. 2009;72(11):809–15.

  48. 48.

    Raz R, Stricker S, Gazzerro E, Clor JL, Witte F, Nistala H, et al. The mutation ROR2W749X, linked to human BDB, is a recessive mutation in the mouse, causing brachydactyly, mediating patterning of joints and modeling recessive Robinow syndrome. Development. 2008;135(9):1713–23.

  49. 49.

    Soccio RE, Adams RM, Romanowski MJ, Sehayek E, Burley SK, Breslow JL. The cholesterol-regulated StarD4 gene encodes a StAR-related lipid transfer protein with two closely related homologues, StarD5 and StarD6. Proc Natl Acad Sci U S A. 2002;99(10):6943–8.

  50. 50.

    Liao R, Sun J, Zhang L, Lou G, Chen M, Zhou D, et al. MicroRNAs play a role in the development of human hematopoietic stem cells. J Cell Biochem. 2008;104(3):805–17.

  51. 51.

    Nelson WJ, Nusse R. Convergence of Wnt, beta-catenin, and cadherin pathways. Science. 2004;303(5663):1483–7.

Download references


The authors would like to acknowledge the members of the Poultry Breeding Group of the College of Animal Science and Technology at the Northeast Agricultural University for managing the birds and collecting the data.


This research was supported by the National 863 Project of China (No. 2013AA102501), the National Natural Science Foundation (No. 31301960), the China Postdoctoral Science Foundation (No. 2015 M581421), Heilongjiang Postdoctoral Financial Assistance (No. LBH-TZ0612), and the University Nursing Program for Young Scholars with Creative Talents in Heilongjiang Province (No. UNPYSCT-2015007).

Availability of data and materials

The data sets supporting the results of this article are included within the article and its additional files. The chicken 60 k SNP data presented in this paper have been deposited into Gene Expression Omnibus ( with the identifier GSE58551.

Author information

HZ analyzed and interpreted the data, drafted, and wrote the manuscript. LLY, JQY and LK participated in the interpretation of data. XYZ and WN participated in the experiments. JR participated in the interpretation of data and helped write the manuscript. HL led the conception and design of the study and helped write the manuscript. All authors submitted comments on drafts, and read and approved the final manuscript.

Correspondence to James M. Reecy or Hui Li.

Ethics declarations

Ethics approval and consent to participate

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.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Additional files

Additional file 1: Figure S1.

Phenotypic distribution of the carcass and growth traits in the lean and fat lines, respectively. The red and light blue bars represent the traits in the lean and fat lines, respectively. The dark blue bars represent the overlap of the traits between the two lines. (PDF 322 kb)

Additional file 2: Table S1.

The 421 pairs of SNPs with significant interaction effects on testis weight (TeW). (XLSX 67 kb)

Additional file 3: Figure S2.

The phenotypic distributions of the four different genotype classes for the most significant AA effect between GGaluGA22768 on GGA3 and Gga_rs14722408 on GGA10. (PDF 122 kb)

Additional file 4: Figure S3.

Epistatic network among SNPs that affect testis weight (TeW). Each node represents a SNP. The chromosome in which a given SNP is located is shown within the circle. A pair of SNPs connected by an edge had a significant interaction. The colors of the nodes represent the P-value of an interaction (P < 1.0 × 10−16 = red; P < 1.0 × 10−15 = blue; P < 1.0 × 10−14 = green; P < 1.0 × 10−13 = white). The color of the edge indicates the type of epistatic effect (AA = red; AD = purple; DA = blue; DD = green). (PDF 1222 kb)

Additional file 5: Figure S4.

The LD blocks of the SNPs on every chromosome in Sub-network1 for testis weight (TeW). This information was used to simply Sub-network1. (PDF 619 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark


  • Carcass and growth traits
  • Testis
  • Epistasis
  • SNP-SNP interaction
  • Chicken