Haplotype-based genome-wide association studies for carcass and growth traits in chicken.

There have been several genome-wide association study (GWAS) reported for carcass, growth, and meat traits in chickens. Most of these studies have been based on single SNPs GWAS. In contrast, haplotype-based GWAS reports have been limited. In the present study, 2 Northeast Agricultural University broiler lines divergently selected for abdominal fat content (NEAUHLF) and genotyped with the chicken 60K SNP chip were used to perform a haplotype-based GWAS. The lean and fat chicken lines were selected for abdominal fat content for 11 yr. Abdominal fat weight was significantly different between the 2 lines; however, there was no difference for body weight between the lean and fat lines. A total of 132 haplotype windows were significantly associated with abdominal fat weight. These significantly associated haplotype windows were primarily located on chromosomes 2, 4, 8, 10, and 26. Seven candidate genes, including SHH, LMBR1, FGF7, IL16, PLIN1, IGF1R, and SLC16A1, were located within these associated regions. These genes may play important roles in the control of abdominal fat content. Two regions on chromosomes 3 and 10 were significantly associated with testis weight. These 2 regions were previously detected by the single SNP GWAS using this same resource population. TCF21 on chromosome 3 was identified as a potentially important candidate gene for testis growth and development based on gene expression analysis and the reported function of this gene. TCF12, which was previously detected in our SNP by SNP interaction analysis, was located in a region on chromosome 10 that was significantly associated with testis weight. Six candidate genes, including TNFRSF1B, PLOD1, NPPC, MTHFR, EPHB2, and SLC35A3, on chromosome 21 may play important roles in bone development based on the known function of these genes. In addition, several regions were significantly associated with other carcass and growth traits, but no candidate genes were identified. The results of the present study may be helpful in understanding the genetic mechanisms of carcass and growth traits in chickens.

Introduction

Single nucleotide polymorphisms (SNP) are the most common type of variant within a genome. They have been extensively used to carry out genome-wide association studies (GWAS). SNP chips have made it possible and affordable to conduct GWAS for complex traits, especially for important economic traits in livestock (Goddard et al., 2016). Therefore, many studies about the successful applications of GWAS in animal breeding and genetics have been reported, and many genes or markers for economically important traits have been identified (Goddard et al., 2016). These results not only supply a number of molecular markers that can be used in prediction/genomic selection but they can also provide important information to help explain the genetic mechanisms that underlie these traits. However, most of these GWAS were based on single SNPs. Single SNP-based GWAS is unlikely to fully capture the variations in regions surrounding the genotyped markers. Instead, haplotype-based GWAS may help to improve this defect and could detect new discoveries of important traits (Howard et al., 2017). In addition, utilization of the haplotype-based approach delivered greater power with no inflation in type I error rate for association studies. The most important process to carry out the haplotype-based GWAS is to construct phasing of the genome, which means that the haplotypes are needed to be constructed. He et al. (2011) developed an efficient approach to accelerate the phasing process and reduce the potential bias generated by unrealistic assumptions in the phasing process. Recently, haplotyped-based GWASs have been conducted and have obtained some useful results (Wu et al., 2014, Sato et al., 2016, Chen et al., 2018). In chickens, GWAS identified genetic variation that has been associated with disease (Raeesi et al., 2017), carcass (Huang et al., 2018), growth (Guo et al., 2017, Pértille et al., 2017), and meat quantitative traits (Moreira et al., 2018). However, nearly all of these GWAS reports were based on single SNP, and no haplotype associations were reported.

The aim of the present study is to identify potentially important genes for carcass and growth traits using a haplotype-based GWAS approach in 2 Northeast Agricultural University broiler lines divergently selected for abdominal fat content (NEAUHLF) for 11 yr. The results of this study may supply useful information for prediction/genomic selection in chicken breeding programs and may also provide important information to explain the genetic mechanisms that underlie carcass and growth traits in chicken.

Materials and methods

Ethics Statement

All animal work was conducted as per 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

Two NEAUHLF were used to carry out the haplotype-based association study (Guo et al., 2011). The population used in the present study included 475 males (203 and 272 birds from the lean and fat lines, respectively) from the 11th generation of NEAUHLF (Li et al., 2013). The birds were weighed at 0, 1, 3, 5, and 7 wk of age (BW0, BW1, BW3, BW5, and BW7, respectively). At 7 wk of age, the metatarsus length (MeL), metatarsus circumference (MeC), keel length (KeL), and chest width (ChWi) were measured before slaughter as previously described (Zhang et al., 2010). Abdominal fat weight (AFW), testis weight (TeW), carcass weight (CW), heart weight (HW), liver weight (LW), spleen weight (SW), and muscular and glandular stomach weight (MGSW) were obtained after the birds were slaughtered.

SNP Genotyping

Genotyping was carried out using the chicken 60 K SNP chip (Illumina Inc., San Diego, CA), which contained 57,636 SNP. After quality control, 48,034 SNP in 475 individuals located on 28 autosomal and Z chromosomes were used in the haplotype-based GWAS. The quality control of the SNP genotypes was described previously by Zhang et al. (2012).

Haplotype-Based GWAS

Haplotypes were constructed by LinkPHASE3 using pedigree information (Druet and Georges, 2015). Missing haplotypes were inferred by DAGPHASE and Beagle, which use an efficient approach based on hidden Markov models (Druet and Georges, 2010). Haplotypes were extracted using every 2 neighboring SNP. Thus, 4 kinds of haplotype (11, 12, 21, and 22) were detected. For the haplotype-based GWAS, we compared each haplotype vs. all others, which means that when haplotype 11 was specified, the individuals with 2 copies of the specified haplotype 11 had the diplotype of AA, the individuals with only one copy of the specified haplotype 11 had the diplotype of AB, and the individuals with no copy of the specified haplotype 11 had the diplotype of BB. In turn, when haplotype 12 was specified, the individuals with 2 copies of the specified haplotype 12 had the diplotype of AA, the individuals with only one copy of the specified haplotype 12 had the diplotype of AB, and the individuals with no copy of the specified haplotype 12 had the diplotype of BB, and so on. The genotype file of all individuals was generated with only 3 diplotypes, AA, AB, and BB. The haplotype-based GWAS was then conducted by Plink v1.07 using a linear regression method (Purcell et al., 2007). Family and Line were used as a 2 fix effects for all the traits to adjust the population structure's effect. BW0 was used as a covariate for BW1, BW3, BW5, and BW7. BW7 was used as a covariate for KeL, MeL, MeC, ChWi, AFW, CW, TeW, HW, LW, SW, and MGSW. A genome-wide 5% type I error after Bonferroni correction was used as the genome-wide significance level. The threshold P-value for declaring genome-wide significance was 0.05/48,005 = 1.04 × 10−6. Manhattan plots of the P-values for all haplotypes associated with carcass and growth were plotted using SNPEVG1, version 2.1 (Wang et al., 2012). Gene locations and information were mined from Ensembl chicken genome galGal3 (https://www.genome.ucsc.edu).

Haplotypes were also extracted using the sliding windows of 3 SNP, 4 SNP, and 5 SNP. The haplotype frequencies were calculated, and the major haplotype was specified, which meant that the individuals with 2 copies of the major haplotype had the diplotype of AA, the individuals with only one copy of the major haplotype had the diplotype of AB, and the individuals with no copy of the major haplotype had the diplotype of BB. Therefore, we got the genotype file of all individuals with only 3 diplotypes, AA, AB, and BB. The haplotype-based GWAS was then conducted by the method described previously.

Results and discussion

Haplotype-Based GWAS for Carcass Trait

For more than 60 yr, broiler chicken breeders have focused on the selection of important economic traits and have made dramatic genetic improvements (Hill and Dansky, 1954, Bedford and Classen, 1992, Demeure et al., 2013). However, long-term intense selection for fast juvenile growth in broiler chickens has increased their abdominal fat deposition and resulted in metabolic changes (Pym, 1987, Emmerson, 1997, Scheele, 1997, Julian, 2005). Excessive deposition of abdominal fat has negative impacts on feed efficiency and carcass quality (Demeure et al., 2013, Ramiah et al., 2014). Therefore, the detection of important genes or markers for abdominal fat content will help to select lean chicken lines. In the present study, haplotype-based GWAS for AFW are carried out to identify genes for abdominal fat content (Figure 1). There were 156 haplotype windows that were significantly associated with AFW (Table 1 and Supplementary Table 1). A total of 132 haplotype windows that were significantly associated with AFW were obtained after combining overlapping windows. The SNP in these significant haplotype windows were concentrated on chromosomes 2, 4, 8, 10, and 26. The 12 regions on these chromosomes were obtained after combining windows that overlapped (Table 2). There were 70 RefGenes located in these 12 regions. Possible candidate genes for abdominal fat deposition include SHH, LMBR1, FGF7, IL16, PLIN1, IGF1R, and SLC16A1. These genes contained a haplotype window or located near a haplotype window with significant effects on AFW (Table 3). Individuals with the major haplotype (Hap1) had significantly lower or higher AFW than the individuals with the other haplotypes (Hap2, Figure 2). These results indicated that SHH, LMBR1, FGF7, IL16, PLIN1, IGF1R, and SLC16A1 are good candidate genes for abdominal fat deposition. SHH (sonic hedgehog) is an obesity susceptibility gene in humans (Wu et al., 2017). This gene can reduce lipid accumulation in adipocytes and decrease the expression of the adipocyte-specific gene (Fontaine et al., 2008). LMBR1 is the limb development membrane protein 1, and the SNP in this gene was significantly associated with obesity in humans (Wu et al., 2017). FGF7 is the fibroblast growth factor (FGF) 7, and the protein encoded by this gene is a member of the FGF family. Most FGF family members could promote the proliferation and differentiation of human preadipocytes by activating a family of receptor tyrosine kinases (Patel et al., 2005). FGF7 was identified as a target of miR-143 in murine adipogenesis and it was plausible that the overexpression of miR-143 could promote adipogenesis by inhibiting its target FGF7 (He et al., 2013). A functional SNP in IL6 gene was strongly associated with waist circumference in a large Dutch study population, which indicated that IL6 may contribute to obesity in humans (van den Berg et al., 2009). Perilipin (PLIN1) is a lipid droplet coat protein that belongs to the lipid droplet–related protein family. Genetic variation in PLIN1 has been significantly associated with adiposity in human (Ruiz et al., 2011), pig (Gandolfi et al., 2011), cattle (Fan et al., 2010), sheep (Gao et al., 2012), duck (Zhang et al., 2013), and chicken (Zhou et al., 2014, Zhang et al., 2015). In mice, knockout of insulin and/or IGF1 receptors (IR/IGF1R) was accompanied by a rapid loss of white and brown fat because of the increased lipolysis and adipocyte apoptosis (Sakaguchi et al., 2017). SLC16A1 is the solute carrier family 16 member 1, which is also known as monocarboxylate transporter 1 (MCT1). MCT1 is abundant in several tissues, including adipose, gut, brain, heart, muscle, liver, and kidney (Hajduch et al., 2000, Pierre and Pellerin, 2005, Iwanaga et al., 2006). It is also a carrier of short-chain fatty acids, ketone bodies, and lactate in several tissues, and MCT1+/− mice displayed resistance to development of diet-induced obesity when fed with high fat diet (HFD) (Lengacher et al., 2013).

Figure 1

Results of haplotype-based genome-wide association studies using PLINK for abdominal fat weight (AFW). The results are presented as Manhattan plots based on haplotype 11-specified, 12-specified, 21-specified, and 22-specified, respectively. The solid line indicates the Bonferroni threshold for multiple test correction with a type I error of 5% (P-value <1.04 × 10−6).

Table 1

Number of haplotype windows with significant effects on each carcass and growth traits in chicken.

Traits	No. of significant windows
	11-Specified	12-Specified	21-Specified	22-Specified	Total
AFW	41	40	50	41	132
BW1	3	2	3	1	9
BW3	4	8	6	3	18
BW5	4	4	7	6	21
BW7	1	2	2	1	6
ChWi	5	3	4	3	14
CW	0	1	1	0	2
HW	1	1	1	0	3
KeL	11	6	9	6	32
LW	0	1	2	0	3
MeC	24	33	35	30	110
MeL	1	2	1	2	6
MGSW	2	0	0	2	4
SW	0	0	0	0	0
TeW	34	33	35	31	123

Abbreviations: AFW, abdominal fat weight; ChWi, chest width; CW, carcass weight; HW, heart weight; KeL, keel length; LW, liver weight; MeC, metatarsus circumference; MeL, metatarsus length; MGSW, muscular and glandular stomach weight; SW, spleen weight; TeW, testis weight.

Table 2

Important chromosome regions for carcass and growth traits.

Chromosome	Start_SNP	Rs#	Start_position	End_SNP	Rs#	End_position	Length	Traits	Genes in the region
2	GGaluGA132691	rs313439121	7896784	Gga_rs15060839	rs15060839	8567871	671,087	AFW	SHH, LMBR1, MNX1, UBE3C
2	Gga_rs14219117	rs14219117	93185343	Gga_rs14219515	rs14219515	93732330	546,987	AFW	/
2	GGaluGA158673	rs312677797	96017288	GGaluGA159074	rs317155927	98822750	2,805,462	AFW	RTTN, MIR1681, TMX3, CDH19, CDH7, MC2R
2	GGaluGA159507	rs315053861	100327421	Gga_rs14224613	rs14224613	100387035	59,614	AFW	/
2	Gga_rs13803296	rs13803296	102079036	GGaluGA160440	rs314547993	103996355	1,917,319	AFW	LAMA1, ZBTB14, AKAIN1, TGIF1, MYL12A
2	Gga_rs16142136	rs16142136	139745278	Gga_rs16141958	rs16141958	140089369	344,091	AFW	/
4	Gga_rs14436487	rs14436487	21729261	Gga_rs14436961	rs14436961	22258381	529,120	AFW	CTSO
8	Gga_rs15906323	rs15906323	8094782	GGaluGA325809	rs431896935	9028904	934,122	AFW	FAM129A
8	Gga_rs14642420	rs14642420	14253680	Gga_rs14642444	rs14642444	14296548	42,868	AFW	ABCD3, ARHGAP29
10	GGaluGA069041	rs317193761	12108078	Gga_rs14008746	rs14008746	14892303	2,784,225	AFW	FGF7, MIR147-1, SLC30A4, BLOC1S6, ITGB1BP3, MIR6596, TRPM7, SPPL2A, GABPB1, HDC, GATM, SCARNA15, FAM103A1, BTBD1, TM6SF1, SH3GL3, EFL1, TMC3, IL16, MESD, ABHD17C, FAH, ZFAND6, BCL2A1, MTHFS, KIAA1024, PLIN1, TICRR, RHCG, FANCI, RLBP1, MFGE8, ACAN, MRPS11, MRPL46, MIR1720, MIR7-2, MIR3529
10	Gga_rs15587351	rs15587351	17309049	Gga_rs14011820	rs14011820	18758907	1,449,858	AFW	NR2F2, MIR1680, MIR1813-2, IGF1R
26	GGaluGA196948	rs314806696	3156806	Gga_rs16203115	rs16203115	3520068	363,262	AFW	KCND3, WNT2B, ST7L, CAPZA1, RHOC, MOV10, SLC16A1, MIR1669, MAGI3
1	Gga_rs13895421	rs13895421	88063956	GGaluGA029830	rs312695192	88670466	606,510	MeC	TBC1D23, TMEM45A, IMPG2, TXNL4B, PCNP
1	GGaluGA031230	rs312759219	92236963	Gga_rs14857266	rs14857266	93018416	781,453	MeC	EPHA3
2	Gga_rs13669384	rs13669384	37647772	Gga_rs14165766	rs14165766	37653833	6,061	MeC	RARB
2	GGaluGA160608	rs318119261	104397754	Gga_rs13794375	rs13794375	104627805	230,051	MeC	/
2	GGaluGA162581	rs431838007	110512137	Gga_rs14232072	rs14232072	111387888	875,751	MeC	MAPRE2, PRKDC, UBE2V2
2	Gga_rs16149569 rs16149569		148367666	GGaluGA173055	rs315266923	150658423	2,290,757	MeC	COL22A1
4	Gga_rs16404447	rs16404447	49102957	Gga_rs14727013	rs14727013	49961124	858167	MeC	MIR1730
6	Gga_rs14593228	rs14593228	32615659	GGaluGA305949	rs312809174	33367322	751,663	MeC	IKZF5, ACADSB, HMX3, BUB3
7	GGaluGA314140	rs315499140	18350464	GGaluGA314144	rs314108745	18364373	13,909	MeC	SP3
7	Gga_rs15862567	rs15862567	24373718	GGaluGA316074	rs312654261	24647515	273,797	MeC	/
8	Gga_rs13663151	rs13663151	7859868	GGaluGA325359	rs316684405	7906123	46,255	MeC	/
8	Gga_rs15910167	rs15910167	10137424	Gga_rs14641638	rs14641638	12972931	2,835,507	MeC	C8H1orf27, AMY1AP, AMY1A, MIR6561, MIR1610, SLC30A7, CDC14A, MFSD14A, SLC35A3, DBT, SASS6, PALMD
21	GGaluGA184599	rs316833978	4781321	Gga_rs15185019	rs15185019	6176965	1,395,644	MeC	MINOS1, NBL1, HTR6, PLA2G2E, PLA2G5, UBXN10, PLA2G2A, DDX19B, DNAJC16, CDA, AGMAT, CTRC, C1orf158, DHRS3, TNFRSF1B, TNFRSF8, PLOD1, CELA2A, NPPC, MTHFR, RNP, NPPA, CLCN6, DRAXIN, MAD2L2, DISP3, GUCA2A, EPHB2
Z	Gga_rs14689552	rs14689552	6673737	Gga_rs14783328	rs14783328	7291085	617,348	MeC	UBE2R2, LOC407092, IFNW1, IFNA3, DCAF12
Z	Gga_rs16129856	rs16129856	9391651	Gga_rs14785793	rs14785793	10575355	1,183,704	MeC	TARS, SLC45A2, AMACR, BRIX1, MIR6613, PRLR, IL7R, LMBRD2, SKP2
1	GGaluGA043278	rs315095993	131588419	Gga_rs13936329	rs13936329	131711376	122,957	TeW	/
2	Gga_rs13534898	rs13534898	4517713	GGaluGA131254	rs314054036	4866215	348,502	TeW	ACAA1, MAPKKK3L, MYD88, MIR6610
2	Gga_rs14240062	rs14240062	121959284	Gga_rs14241677	rs14241677	123386087	1,426,803	TeW	TRPA1, MIR1796, TERF1, RPL7, RDH10, STAU2, UBE2W, ELOC, TMEM70, PI15, CRISPLD1
2	Gga_rs14245700	rs14245700	127443603	Gga_rs13730959	rs13730959	127499632	56,029	TeW	CA3A
3	GGaluGA222074	rs317102159	53276944	Gga_rs10729720	rs10729720	67517313	14,240,369	TeW	GTF2H5, EZR, ADGRG6, CITED2, TXLNB, ABRACL, REPS1, MIR7462, PERP2, PERP1, IFNGR1, MIR6568, PEX7, MAP7, MYB, SGK1, TBPL1, TCF21, EYA4, RPS12, MIR1454, SLC18B1, VNN1, STX7, MOXD1, CTGF, MIR6582, MIR6697, MIR1660, ECHDC1, RSPO3, CENPW, TRMT11, NCOA7, TPD52L1, HDDC2, NKAIN2, FABP7, PKIB, SERINC1, HSF2, GJA1, MCM9, ASF1A, PLN, MIR199B, ROS1, VGLL2, SOT3A1L, RWDD1, FAM26E, HDAC2, MARCKS
10	Gga_rs14002765	rs14002765	5962967	Gga_rs14003104	rs14003104	6635581	672,614	TeW	RORA, ANXA2, GTF2A2
10	Gga_rs14695763	rs14695763	8460335	Gga_rs14722408	rs14722408	13180860	4,720,525	TeW	TCF12, PRTG, PYGO1, DYX1C1, CCPG1, PIGBOS1 RAB27A, RSL24D1, FAM214A, ARPP19, MYO5A, GNB5, BCL2L10, MAPK6, LYSMD2, LEO1, TMOD3, LYSMD2, SCG3, CYP19A1, MIR1744, SLC24A5, MYEF2, DUT, COPS2, GALK2, FGF7, MIR147-1, BLOC1S6, ITGB1BP3, MIR6596, GABPB1, TRPM7,GABPB1, HDC, GATM, SCARNA15, FAM103A1, FAM103A1, TM6SF1, BTBD1, SH3GL3
11	GGaluGA074107	rs312924990	1035483	Gga_rs14958653	rs14958653	1864531	829,048	TeW	CTCF, LOC415664, LOC415664, LOC769668, LOC107080643, LOC415662, AARS, MIR1616, FHOD1, ATP6V0D1, AGRP, SETD6, CNOT1, GOT2, CALB2, HYDIN, VAC14, COG4, ST3GAL2, GLG1

Table 3

Candidate genes for AFW, TeW, and MeC identified from the haplotype-based GWAS results.

Genes	Haplotype window	Near or contained the haplotype window	Chromosome	Major haplotype	Trait
SHH	WIN7856	Near	2	21	AFW
LMBR1	WIN7879	Near	2	12	AFW
FGF7	WIN30447	Near	10	22	AFW
IL16	WIN30558	Near	10	11	AFW
PLIN1	WIN30605	Near	10	12	AFW
IGF1R	WIN30893	Contained	10	22	AFW
SLC16A1	WIN44687	Near	26	12	AFW
TCF21	WIN15421	Near	3	11	TeW
TCF12	WIN30233 and WIN30234	Contained	10	212	TeW
SLC35A3	WIN27613	Contained	8	12	MeC
TNFRSF1B	WIN42161	Near	21	22	MeC
PLOD1
NPPC	WIN42177	Near	21	12	MeC
MTHFR
EPHB2	WIN42231 and WIN42234	Contained	21	1,212	MeC

Abbreviations: AFW, abdominal fat weight; GWAS; genome-wide association study; MeC, metatarsus circumference; TeW, testis weight

Figure 2

The difference of abdominal fat weight (AFW) between the individuals with the major haplotype (Hap1) and the individuals with other haplotypes (Hap2) (t-test). Different alphabets means extremely significantly different (P < 0.01) and the error bar is the standard deviation (SD).

Manhattan plots of haplotype-based GWAS for TeW are shown in Figure 3. A total of 133 haplotype windows significantly associated with TeW were identified (Table 1 and Supplementary Table 1). These significant windows for TeW were mainly distributed on chromosomes 3 and 10. The haplotype windows with a significant effect on TeW on chromosome 3 were concentrated on a 14 Mb region from 53.28 Mb to 67.52 Mb. The significant haplotype windows for TeW on chromosome 10 were concentrated on the 4.72 Mb region from 8.46 Mb to 13.18 Mb. These 2 regions on chromosome 3 and 10 are same as previously detected by the single SNP GWAS (Zhang et al., 2017a). In these 2 regions, 2 transcription factors, including TCF21 and TCF12, were detected as important genes for testis growth and development based on our previous studies (Zhang et al., 2017a, Zhang et al., 2017b). TCF21 gene was located near a haplotype window (WIN15421) that was significantly associated with TeW (Table 3). Individuals with the major haplotype 11 (Hap1) had significantly lower TeW than individuals with the others haplotypes (Hap2) (Figure 4). Previously reported gene expression analysis indicated that TCF21 was differently expressed between lean and fat birds and that its expression level was significantly associated with TeW and TeP (Zhang et al., 2017a). In humans and mice, TCF21 plays important roles in hypertension, gastric cancer, and coronary heart disease (Miller et al., 2014, Fujimaki et al., 2015, Yang et al., 2015). In mice, TCF21 is the first direct downstream target gene of the male sex–determining factor (SRY) (Bhandari et al., 2011, Bhandari et al., 2012). The knockout of TCF21 in mice resulted in male-to-female sex reversal (Cui et al., 2004). SRY could bind to the TCF21 promoter and activate gene expression (Bhandari et al., 2012). In rats, TCF21 and SRY have similar effects on Sertoli cell differentiation and embryonic testis development (Bhandari et al., 2012). Taken together, these results indicated that TCF21 may play an important role in sex differentiation and testis development. TCF12 was located within 2 consecutive haplotype windows (WIN30233 and WIN30234) that were significantly associated with TeW (Table 3). The 3 SNP that constituted these 2 haplotypes were used to construct 3 SNP haplotypes. Individuals with the major haplotype 212 (Hap1) had significantly higher TeW than the individuals with others haplotypes (Hap2) (Figure 4). TCF12 was in the same family as TCF21, which was also identified in the region for TeW on chromosome 10. In our previous study, TCF12 was detected as the important gene for testis growth and development from the SNP by SNP interaction analysis (Zhang et al., 2017b).

Figure 3

Results of haplotype-based genome-wide association studies using PLINK for testis weight (TeW). The results are presented as Manhattan plots based on haplotype 11-specified, 12-specified, 21-specified, and 22-specified, respectively. The solid line indicates the Bonferroni threshold for multiple test correction with a type I error of 5% (P-value <1.04 × 10−6).

Figure 4

The difference of testis weight (TeW) between the individuals with the major haplotype. (Hap1) and the individuals with other haplotypes (Hap2) of TCF21 and TCF12 genes (t-test). ∗means significantly different (P < 0.05) and the error bar is the standard deviation (SD).

For TeW, a single SNP-based GWAS was carried out, previously (Zhang et al., 2017a). The haplotype-based GWAS results were compared with the single SNP-based GWAS, and we found that haplotype-based GWAS identified all significant regions detected by single SNP-based GWAS for TeW. Furthermore, haplotype-based GWAS detected more significant regions for TeW than single SNP-based GWAS. Such significant regions on chromosomes 1 and 11 for TeW in the present study (Table 2) were not detected by single SNP-based GWAS as previously reported (Zhang et al., 2017a). Therefore, from these results we could conclude that the haplotype-based GWAS is a good supplement for single SNP-based GWAS.

For CW, HW, LW, SW, and MGSW, only a couple of haplotypes were significantly associated. Unfortunately, no interesting candidate genes were detected for these carcass traits (Table 1, Supplementary Table 1, and Supplementary Figure 1).

Haplotype-Based GWAS for Growth Trait

Manhattan plots of haplotype-based GWAS for MeC are shown in Figure 5. There were 122 haplotype windows that were significantly associated with MeC (Table 1 and Supplementary Table 1). A total of 110 haploytpe windows were obtained after deleting the overlapped windows. Most of these significant haploptypes were distributed on chromosomes 1, 2, 8, 21, and Z. There were 66 RefGenes located in these regions, and possible candidate genes for bone traits include TNFRSF1B, PLOD1, NPPC, MTHFR, EPHB2, and SLC35A3. These genes were contained within or near a haplotype window that was significantly associated with MeC (Table 3). For each gene, individuals with the major haplotype (Hap1) had significantly lower (or higher) MeC than the individuals with the other haplotypes (Hap2, Figure 6). EPHB2 spanned 2 haplotype windows (WIN42231 and WIN42234). These 4 SNP that constituted these 2 windows were used to construct 4 SNP haplotypes. Individuals with the major haplotype 1212 (Hap1) had significantly lower MeC than the individuals with the other haplotypes (Hap2) (Figure 6C). These results indicated that TNFRSF1B, PLOD1, NPPC, MTHFR, EPHB2, and SLC35A3 are important for bone development. TNFRSF1B is a TNF receptor superfamily member, which could regulate the effects of TNF on osteoclastogenesis (Abu-Amer et al., 2000). The SNP in TNFRSF1B could contribute to the genetic regulation of bone mass (Albagha et al., 2002). PLOD1 is procollagen-lysine, 2-oxoglutarate 5-dioxygenase 1. Variants within this gene have been associated with bone mineral density (BMD) in humans (Spotila et al., 2003, Huang et al., 2009). NPPC is C-type natriuretic peptide 3, which is also known as CNP. Mice that overexpress CNP have longer bones (Chusho et al., 2001). CNP could stimulate chondrocyte proliferation and increase the size of individual hypertrophic chondrocytes (Yasoda et al., 1998, Mericq et al., 2000). CNP has been implicated in the regulation of skeletal growth in transgenic and knockout mice (Bartels et al., 2004). MTHFR is methylenetetrahydrofolate reductase, which catalyzes the conversion of 5,10-methylenetetrahydrofolate to 5-methyltetrahydrofolate, a cosubstrate for homocysteine remethylation to methionine. Variants within this gene have been associated with BMD (Li et al., 2016). EPHB2 is EPH receptor B2. The GWAS meta-analysis of lumbar spine volumetric BMD measured by quantitative computed tomography was carried out and several loci were identified, including rs12742784 within EPHB2, which was associated with higher volumetric BMD and decreased risk of clinical vertebral fracture (Nielson et al., 2016). This noncoding SNP has been associated with increased EPHB2 mRNA expression levels in human bone biopsies (Nielson et al., 2016). The basic function of SLC35A3 is as a UDP-GlcNAc transporter. It has been shown to be expressed in all human tissues examined, including mesodermal derived tissues, skeletal muscle, and bone marrow (Ishida et al., 1999). A missense mutation in SLC35A3 gene has been associated with complex vertebral malformations in bovine and revealed a new mechanism for malformation of the vertebral column caused by abnormal nucleotide sugar transport into the Golgi apparatus (Thomsen et al., 2006). Some other studies have also identified SLC35A3 as having an important role in vertebral malformations (Ghebranious et al., 2006, Ruść and Kamiński, 2007, Chu et al., 2008, Ghanem et al., 2008, Ghanem et al., 2009, Wang et al., 2011).

Figure 5

Results of haplotype-based genome-wide association studies using PLINK for metatarsus circumference (MeC). The results are presented as Manhattan plots based on haplotype 11- specified, 12-specified, 21-specified, and 22-specified, respectively. The solid line indicates the Bonferroni threshold for multiple test correction with a type I error of 5% (P-value <1.04 × 10−6).

Figure 6

The difference of metatarsus circumference (MeC) between the individuals with the major. haplotype (Hap1) and the individuals with other haplotypes (Hap2) (t-test). ∗∗means extremely significantly different (P < 0.01) and the error bar is the standard deviation (SD).

For BW1, BW3, BW5, BW7, ChWi, KeL, and MeL, only a couple of haplotype windows were significantly associated with these traits. No potential candidate genes were detected for these growth traits (Table 1, Supplementary Table 1 and Supplementary Figure 1).

Haplotype-Based GWAS Using Sliding Window of 3 SNP, 4 SNP, and 5 SNP

The GWAS results for carcass and growth traits aforementioned were all based on haplotypes extracted from sliding windows of 2 neighbor SNP. We also constructed haplotypes using 3 SNP in a sliding window, 4 SNP in a sliding window, and 5 SNP in a sliding window. Accordingly, haplotypes-based GWAS were carried out using 3-SNP, 4-SNP and 5-SNP sliding windows, respectively. Manhattan plots of 3-SNP, 4-SNP, and 5-SNP windows for carcass and growth traits are shown in Supplementary Figure 2. These results are similar as the results of 2-SNP window described previously.

In summary, the present study successfully used the haplotype-based GWAS method to detect important chromosome regions that harbor genes associated with carcass and growth traits in chicken. SHH, LMBR1, FGF7, IL16, PLIN1, IGF1R, and SLC16A1 were identified as potential candidate genes for abdominal fat deposition. TCF21 and TCF12, which were also previously detected by single SNP GWAS and epistatic effect analysis, were detected as important candidate genes for testis growth and development. TNFRSF1B, PLOD1, NPPC, MTHFR, EPHB2, and SLC35A3 were potentially important genes for bone development. Only a couple of regions were detected as significantly associated with other carcass and growth traits. The results of this study may be helpful for exploring the metabolic mechanisms of fat deposition and testis growth in chicken.