Genome wide screening of candidate genes for improving piglet birth weight using high and low estimated breeding value populations.

Birth weight is an economically important trait in pig production because it directly impacts piglet growth and survival rate. In the present study, we performed a genome wide survey of candidate genes and pathways associated with individual birth weight (IBW) using the Illumina PorcineSNP60 BeadChip on 24 high (HEBV) and 24 low estimated breeding value (LEBV) animals. These animals were selected from a reference population of 522 individuals produced by three sires and six dam lines, which were crossbreds with multiple breeds. After quality-control, 43,257 SNPs (single nucleotide polymorphisms), including 42,243 autosomal SNPs and 1,014 SNPs on chromosome X, were used in the data analysis. A total of 27 differentially selected regions (DSRs), including 1 on Sus scrofa chromosome 1 (SSC1), 1 on SSC4, 2 on SSC5, 4 on SSC6, 2 on SSC7, 5 on SSC8, 3 on SSC9, 1 on SSC14, 3 on SSC18, and 5 on SSCX, were identified to show the genome wide separations between the HEBV and LEBV groups for IBW in piglets. A DSR with the most number of significant SNPs (including 7 top 0.1% and 31 top 5% SNPs) was located on SSC6, while another DSR with the largest genetic differences in F ST was found on SSC18. These regions harbor known functionally important genes involved in growth and development, such as TNFRSF9 (tumor necrosis factor receptor superfamily member 9), CA6 (carbonic anhydrase VI) and MDFIC (MyoD family inhibitor domain containing). A DSR rich in imprinting genes appeared on SSC9, which included PEG10 (paternally expressed 10), SGCE (sarcoglycan, epsilon), PPP1R9A (protein phosphatase 1, regulatory subunit 9A) and ASB4 (ankyrin repeat and SOCS box containing 4). More importantly, our present study provided evidence to support six quantitative trait loci (QTL) regions for pig birth weight, six QTL regions for average birth weight (ABW) and three QTL regions for litter birth weight (LBW) reported previously by other groups. Furthermore, gene ontology analysis with 183 genes harbored in these 27 DSRs suggested that protein, metal, ion and ATP binding, viral process and innate immune response present important pathways for deciphering their roles in fetal growth or development. Overall, our study provides useful information on candidate genes and pathways for regulating birth weight in piglets, thus improving our understanding of the genetic mechanisms involved in porcine embryonic or fetal development.

Introduction

Birth weight is an economically important trait closely associated with piglet survival and growth after birth. For example, low birth weight in piglets often results in high mortality rates, slow growth and poor carcass quality 1-5. Generally speaking, the birth weight of a newborn piglet reflects the progress of its embryonic and fetal development, which mainly depends on the nutritional supply to the embryo/fetus and its ability to use the available substrates 6. On the other hand, individual birth weight (IBW) has a substantial genetic component. For example, the direct heritability for IBW ranges from 0.13 to 0.36 5, 7. Therefore, identification of IBW related candidate genes or their pathways are imperative for developing marker assisted selection to rapidly improve this phenotype in pigs.

In the past, several research groups have focused on candidate genes associated with piglet IBW, such as myogenin (MYOG) 8, myostatin (MSTN) 9, dopamine beta-hydroxylase (DBH) 10, POU class 1 homeobox 1 (POU1F1) 11, solute carrier family 27 (fatty acid transporter), member 4 (SLC27A4) 12, solute carrier family 6 (amino acid transporter), member 14 (SLC6A14) 13, proliferation associated nuclear element 1 (PANE1) 14, and secreted phosphoprotein 1 (SPP1) 15. Recent completion of the pig draft genome assembly and discovery of more than twenty million single nucleotide polymorphisms (SNPs) in the genome (http://www.ncbi.nlm.nih.gov/SNP/) have allowed the community to further advance genome wide association studies for the identification of QTLs related to many economically important traits in the species.

For instance, the Illumina PorcineSNP60 BeadChip is a powerful tool for identifying SNPs or genes associated with economically important quantitative traits, such as growth, body composition and structural soundness 16-17, reproduction and fertility 18-19, disease resistance 20-21, farrowing phenotypes 22, meat quality and fat deposition traits 23-25. Due to commercial services available to genotype these chips, such genome wide association studies can be done in a short period of time, but the cost can be high, particularly when a large number of samples are involved.

To date, only one group has focused on the applications of the Illumina PorcineSNP60 BeadChip in investigating quantitative trait loci (QTL) for litter birth weight (LBW) and average birth weight (ABW) 22. Using a total of 1,152 first parity gilts, the team found 33 QTL for LBW and 65 QTL for ABW, respectively 22. Therefore, the objective of this study was to use both high estimated breeding value (HEBV) and low estimated breeding value (LEBV) piglets for IBW to perform genome-wide analysis using the PorcineSNP60 BeadChip for identification of candidate genes and pathways associated with this economically important trait. Our results provide novel insights in molecular mechanisms involved in genetic differences of fetal development in piglets.

Materials and Methods

Animals and DNA Isolation

A composite population was developed at the Washington State University Swine Center during 2003-2004. When piglets were born, tails were docked and IBW were recorded. Tail docking is an industry routine practice as removal of part of the pigs' tail can reduce the risk of tail biting in older pigs. Tail biting is painful and can lead to serious injury, which might result in symptoms of physical or mental stress. As such, these previously collected tails were used for DNA extraction in the present study. The tissue was lysed and DNA prepared with the GenEluteTM Mammalian Genomic DNA Miniprep kit (Sigma, St. Louis, MO, USA) according to the manufacturer's instructions. The concentration and purity of DNA samples were measured using the NANODROP® spectrophotometer ND-1000 (Thermo Scientific, Waltham, MA, USA). Overall, the DNA concentrations, A260/280 ratios and A260/230 ratios ranged from 40.6 to 76.3 ng/μl, from 1.85 to 1.94 and from 1.98 to 2.31, respectively. Based on the pedigree information, we selected a total of 522 piglets derived from three sires and six dam lines to form our reference population. Both sire and dam lines were crossbreds, representing six breeds: Yorkshire, Landrace, Chester White, Large White, Hampshire, and Duroc.

Estimated Breeding Value (EBV) Analysis and Selection of Piglets for Genotyping on the Illumina PorcineSNP60 BeadChip

The EBV analysis for IBW was performed using the following mixed linear model:

Where y is the observation of birth weight, b is a vector including the systematic effects of year-season (5), sex (2), parity (3), dam line (6) and sire (3). The vectors d, c, and e represent the direct additive genetic effects, common environmental litter effects (56), and the environmental residual effects, respectively. All computations were performed using AI-REML method with DMU v6.0 software 26. Based on the distribution of estimated breeding values for piglet IBW within 6 dam lines sired by 3 crossbred boars, 48 piglets with high (HEBV) and low estimated breeding values (LEBV) were selected to perform genome wide screening of candidate genes and pathways associated with piglet IBW. A student's t-test was used to determine the significant difference between HEBV and LEBV groups. Therefore, a total of 48 DNA samples were genotyped with standard procedures at GeneSeek (Lincoln, NE, USA) on the PorcineSNP60 BeadChip. The genotype quality control was done according to our previous report 27.

Detection of Differentially Selected Regions (DSRs) for IBW between HEBV and LEBV Piglets

As White et al (2012) 28 suggested, Hardy-Weinberg Equilibrium (HWE) of polymorphic SNPs was set to 0.000001 with PLINK v1.06 data analysis tools 29. The DSR algorithm was described previously 30-31, but with slight modifications: 1) raw values were ranked and used to identify regions; 2) Fisher's exact test was executed in R 2.14.0 to compare the allele frequencies between HEBV and LEBV, and SNPs with P values < 0.05 were considered statistically significant; 3) estimation of SNP FST were based on the model proposed by Nicholson et al (2002) 32 and Flori et al (2009) 33, and then the significant SNPs with 0.1% or 5% highest FST values were selected; and 4) by placing the top significant SNP(s) (0.1%) in the center of a DSR, adjacent SNP markers were collected to determine the region boundaries until more than three consecutive SNPs were not in the top significant 5% threshold. When DSRs overlapped, we combined them. In addition, when a region contained more than ten SNPs in the top significant 5%, but none were in the top significant 0.1%, we also considered it as a DSR region in our study.

SNP-specific FST values were smoothed over each chromosome with a local variable bandwidth kernel estimator 34. Genes in these DSRs were identified using the Sscrofa 10.2 assembly (http://www.animalgenome.org/cgi-bin/gbrowse/pig). Gene ontology information was collected from the NCBI database based on the summary of human gene functions.

Results and Discussion

Illumina PorcineSNP60 BeadChip Genotyping Basics

As described above, a total of 48 animals were selected for genotyping using the Illumina PorcineSNP60 BeadChip, including 24 HEBV and 24 LEBV piglets. The average IBW was 1.64 kg for all 522 piglets in our reference population. As shown in Table 1, the HEBV piglets were 0.65 kg heavier while LEBV animals were 0.63 kg lighter than the population average. As such, the 1.28 kg difference in birth weight between HEBV and LEBV was significant (t-test, P < 0.01). Genotyping revealed that among 61,565 SNPs on the Illumina PorcineSNP60 BeadChip, only 54,238 SNPs were exactly assigned to a unique chromosome position based on the Sus scrofa Build 10.2 assembly. Meanwhile, we observed that 3,725 SNPs failed to be scored on at least 95% of all the individuals, 7,235 SNPs had MAF ≤ 0.05 in the whole dataset, and 21 SNPs showed significant deviations from Hardy-Weinberg equilibrium (HWE) in any one of two groups. By excluding the SNPs described above, the remaining 43,257 SNPs, including 42,243 autosomal SNPs and 1,014 SNPs on chromosome X, were used for further analysis (Additional File 1: Figure S1). The HWE test was not performed on chromosome X because the SNPs on chromosome X in males carry only one copy.

Table 1

EBVs and IBW in the reference population.

Population	Number of piglets	EBVs	IBW (kg)
All pigs	522	-0.03 ± 0.40	1.64 ± 0.35
HEBV	24	0.73 ± 0.17a	2.29 ± 0.14a
LEBV	24	-0.80 ± 0.12b	1.01 ± 0.18b
P		< 0.01	< 0.01

HEBV: high estimated breeding value population; LEBV: low estimated breeding value population; IBW: individual birth weight. Values are shown as mean ± standard deviation. A student's t-test was used to compare the significant difference between HEBV and LEBV in EBVs and IBW, respectively. P value < 0.05 or 0.01 is considered statistically significant or extremely significant level, respectively. The different lowercase letters between HEBV and LEBV indicate that the difference reached the significance level of P < 0.05.

DSRs Discovered for IBW between HEBV and LEBV Groups

As shown in Table 2 and Figure 1A, a total of 27 DSRs across the entire pig genome were differentiated between the HEBV and LEBV groups for IBW, which included 43 significant top 0.1% SNPs (33 SNPs located on the autosomal chromosomes and 10 SNPs on chromosome X). Based on the current pig genome assembly (Sscrofa 10.2), these DSRs harbored a total of 183 unique genes: 168 of them are found within the autosomal DSRs and 15 genes were included in the chromosome X DSRs (Additional File 2: Table S1). Among these 27 DSRs, five supported the QTL regions previously reported for the birth weight in pigs: QTL#5232 on SSC1, QTL#369 on SSC4, QTL#3189 on SSC5, QTL#1012 on SSC6, QTL#187 and QTL#5197 on SSC 7 (Table 2 and Figure 2, http://cn.animalgenome.org/cgi-bin/QTLdb/index). More interestingly, six and three of these 27 DSRs for IBW were also located in relative close vicinity of QTLs identified for average birth weight (ABW) and litter body weight (LBW), respectively 22 (Additional File 2: Table S2). These results clearly indicated that some genome regions might be commonly responsible for IBW, ABW and LBW in piglets.

Table 2

Differentially selected regions between high and low estimated breeding value populations.

Region	PorcineSNP60 BeadChip Position (Mb) based on Sscrofa 10.2 assembly	Peak SNP (FST)	Top 0.1%	Top 5%	Candidate genes	Birth weight QTL number
1	chr.1:163614679..163648076	ASGA0004841(0.1031)	2	0	NA	5232
2	chr.4:33947049..34338016	ALGA0024390(0.1129)	2	6	NA	369
		ASGA0019178(0.1129)
3	chr.5:3829983..5355995	M1GA0007314(0.1073)	1	10	NA	NA
4	chr.5:6836955..7331522	M1GA0007436(0.1409)	1	3	TMEM184B	3189
5	chr.6:13486031..13812780	H3GA0017575(0.1031)	1	4	NA	NA
6	chr.6:21038055..21721506	ALGA0034850(0.1082)	1	7	NA	NA
7	chr.6:37040339..37139826	H3GA0054966(0.1293)	1	2	NA	NA
8	chr.6:62369260..64915744	ASGA0095271(0.1300)	7	31	TNFRSF9, CA6	1012
		ALGA0117367(0.1300)
		MARC0067004(0.1300)
9	chr.7:11223472..11802513	ASGA0031182(0.1031)	1	1	NA	NA
10	chr.7:57441888..58806078	ALGA0042164(0.0813)	0	12	CCDC37	187, 5197
11	chr.8:20367972..21000928	MARC0112253(0.1268)	2	3	TBC1D19	NA
12	chr.8:93481336..94303744	MARC0093552(0.1163)	1	1	NA	NA
13	chr.8:130415971..131295589	ASGA0039827(0.1143)	2	3	ADH4	NA
14	chr.8:139592943..139593441	MARC0054584(0.1163)	1	1	NA	NA
15	chr.8:142867557..143732188	ALGA0050003(0.1282)	2	7	WDFY3	NA
16	chr.9:60085650..61556267	ASGA0102114(0.1185)	3	10	NA	NA
17	chr.9:81214756..82621184	ALGA0053850(0.1258)	1	3	LOC100738652	NA
18	chr.9:95825746..97051426	ALGA0054166(0.0857)	0	16	HDAC9	NA
		CASI0007446(0.0857)
		ASGA0043971(0.0857)
19	chr.14:9010230..10254612	ALGA0074932(0.1014)	1	10	NA	NA
20	chr.18:30963795..31174385	ALGA0097763(0.1259)	1	1	NA	NA
21	chr.18:32516193..33286269	INRA0055670(0.1420)	1	2	MDFIC	NA
22	chr.18:55026817..55614254	ALGA0098742(0.1204)	1	2	NA	NA
23	chr.X:8818096..10005053	MARC0040504(0.1077)	1	9	NA	NA
24	chr.X:15592452..16389800	ASGA0080878(0.1282)	2	2	NA	NA
25	chr.X:19357189..21123895	H3GA0051592(0.1675)	1	10	NA	NA
26	chr.X:39671416..40276593	ASGA0081063(0.1007)	5	3	BCOR	NA
		M1GA0023676(0.1007)
		H3GA0051711(0.1007)
		H3GA0051713(0.1007)
		MARC0069431(0.1007)
27	chr.X:136986704..137562054	MARC0046821(0.1300)	1	1	NA	NA

Note: The candidate genes are given within the SNPs of top significant 0.1% or peak SNP for each region.

Figure 1

Genome-wide distribution of SNPs between high and low estimated breeding value populations in piglets. (A) A global view of FST value distributions on pig autosomes 1 - 18 and sex chromosome X. A representative DSR (region 8 of SSC6) with 7 top significant 0.1% and 31 top significant 5% SNPs was indicated. (B) Smoothed FST estimates in region 8 of SSC6 showing the strong selection signals between the HEBV and LEBV groups. (C) Gene allele frequencies of each locus in region 8 of SSC6 between the HEBV and LEBV groups. HEBV: high estimated breeding value; LEBV: low estimated breeding value.

Figure 2

The distribution of selected genomic regions (DSRs) and QTL for birth weight across the pig genome (based on Sus scrofa 10.2 Assembly). All QTL information for birth weight was retrieved from the Animal QTLdb release 21 (http://cn.animalgenome.org/cgi-bin/QTLdb/index). A total of 44 QTL from Animal QTLdb database and 27 DSRs of this study for IBW were placed on each pig chromosome based on their locations. The QTLs were painted in blue along each chromosome and DSRs were presented in red bars, respectively.

DSR with the Richest Top SNP for IBW Located on SSC6

A DSR region with 7 top significant 0.1% and 31 top significant 5% SNPs was found on SSC6 (named region 8, 62.37 to 64.92 Mb) (Table 2 and Figure 1B). The allele frequency differences among SNPs along the region are illustrated in Figure 1C. A major QTL for ABW was also discovered in this region using microsatellite markers 35 (Figure 2) as well as using the genome-wide association study 22 as described above. Two important candidates, namely tumor necrosis factor receptor superfamily member 9 (TNFRSF9) and carbonic anhydrase VI (CA6), are located in this region. TNFRSF9 is a member of the tumor necrosis factor receptor superfamily and it functions as an inducible co-stimulatory molecule for T cells 36. Recently, TNFRSF9 was found to have an important function during the implantation period of pregnancy in mice, which is a necessary process for the mother to adequately provide nutrients for fetal development 37. Moreover, TNFRSF9 plays an important role in reversing the hypo-responsiveness of pathogen-reactive maternal CD8+ T cells in the placenta 38. CA6 serves an important physiological role in the homeostasis of oral tissue by catalyzing the hydration of carbon hydroxide in body fluid and is the only isoform secreted in saliva, tears and milk of mammals. CA6 polymorphisms are known to be associated with the concentrations of secreted CA6 in young people, as well as salivary buffer capacity and dental plaque pH in children 39-40. There is currently limited information about the roles of TNFRSF9 and CA6 in livestock, and thus confirming the roles of these two genes in porcine fetal development warrants further investigation.

DSR with the Highest FST Signal for IBW Located on SSC18

Among autosomal DSRs, the highest FST signal (INRA0055670, FST = 0.1420, region 21; Table 2) was located at 33.29 Mb on SSC18 (Sscrofa 10.2 assembly) where the MyoD family inhibitor domain containing (MDFIC) gene resides. As a transcriptional regulator, this gene is a member of a family of proteins that include a specific cysteine-rich C-terminal domain for transcriptional regulation of viral genome expression 41 and have a modulatory role in immune cells 42, which, taken together, suggest that the function of this gene might be related to immune system capabilities. However, little is known about the role of MDFIC in livestock to date. Our results suggest that MDFIC might affect IBW, but additional studies are needed to confirm this speculation.

DSR with Four Imprinting Genes for IBW Located on SSC9

As shown in Additional File 2: Table S1, a DSR region of 81.21 - 82.62 Mb on SSC9 (Region 17; Table 2) harbors 4 putative candidate genes: PEG10 (paternally expressed 10), SGCE (sarcoglycan, epsilon), PPP1R9A (protein phosphatase 1, regulatory subunit 9A) and ASB4 (ankyrin repeat and SOCS box containing 4). These are imprinting genes that are associated with placental or embryonic growth in mouse 43. Therefore, we present here for the first time, to our knowledge, a potential imprinting region for IBW in pigs on SSC9. PEG10 and SGCE are candidate genes for embryonic lethality 44-45, while PPP1R9A and ASB4 play important roles in early development of multiple tissue types, which affect embryo growth and development in mouse and cattle 46-48. For pig, PEG10 may regulate the expression status in diploid parthenogenetic embryos and alter methylation at an imprinting control region 49. Furthermore, imprinting has been previously demonstrated for PEG10 and SGCE in porcine placenta at day 30 50.

Important candidate genes involved in DSRs on SSC9 and SSCX

Two DSRs, one located in region 95.83 - 97.05 on SSC9 and another in region 39.67 - 40.28 on SSCX, appear to be important in early development. The former DSR contains HDAC9 (histone deacetylase), while the latter contains BCOR (BCL6 corepressor). HDCA9 is a member of the HDAC family of proteins that function in deacetylation of histones, which change chromatin structure and gene transcription, and thus regulate different development and differentiation processes 51-52. This gene is sensitized to hypertrophic signals and exhibits stress-dependent cardiomegaly 53 and can suppress the transcriptional activity of myocyte enhancer factor (MEF2), which is a key regulator of muscle development in embryonic mice that forms a negative-feedback loop in the transcriptional circuitry of muscle differentiation 54. Moreover, the gene is also strongly expressed in skeletal muscle and heart during murine embryogenesis 54. Taken together, these data suggest that HDAC9 might play a role in cardiogenesis and/or myogenesis in pigs as well. The BCOR gene is a key transcriptional regulator of embryonic stem cell differentiation and early embryonic development, in particular in the differentiation of multiple tissue lineages or embryonic stem cells into ectoderm, mesoderm and hematopoietic lineages 55-56. Furthermore, in situ analysis showed that BCOR expression is significantly increased throughout the embryo after embryonic turning in mouse 55, suggesting a significant role in embryonic development. In addition, BCOR is differentially expressed in diverse tissue lineages such as eye, brain, neural tube and branchial arches during organogenesis and fetal period, which correlate with tissues affected in oculofaciocardiodental (OFCD) patients 55. Most recently, a mutation in BCOR was found to be associated with OFCD syndrome in human 57. Interestingly, our results suggested HDCA9 and BCOR genes might be associated with IBW, but how the tissue-specific expression patterns of these two genes specifically affect IBW during piglet embryogenesis requires further investigation.

A View of Gene Ontology from Candidate Genes

As indicated above, our present study revealed a total of 183 genes that reside in 27 DSRs for IBW in piglets (Additional File 2: Table S1). Gene ontology analysis revealed that these genes are mainly involved in functional or biological processes related to protein/metal/calcium/zinc ion binding, ATP binding, viral process and innate immune response (Figure 3 and Additional File 2: Table S3 and S4). This is not surprising because these processes are necessary to meet the requirements of embryonic or fetal development. For example, zinc is an essential element during early murine morphogenesis 58, and its deficiency causes congenital malformations in rat fetuses 59. Zinc deficiency can affect birth weight of offspring of rats or ewes when the degree of deficiency is severe 60-61. Moreover, zinc supplementation had significant effects on birth weight in pregnant adolescents who were underweight multiparas or from poor urban communities 62-64. Adequate maternal calcium intake is significantly correlated with appropriate birth weight in newborns 65 and a suitable concentration of calcium may meet fetal bone requirements in human 66. Our results provide further evidence that calcium/zinc are important factors associated with IBW in pigs. However, the influence of maternal calcium level on human infant growth remains unclear because other studies did not show this positive correlation 67-68. Therefore, additional studies need to be conducted to confirm the exact role of calcium ions in pig fetal development. On the other hand, viral process and innate immune response are also important to developing fetuses. Actually, low birth weight is often related to intrauterine growth restriction (IUGR), which is defined as impaired growth and development of the mammalian embryo or other organs during pregnancy and is found in multiple mammalian species 69. Up to now, several studies have discovered that IUGR decreased the levels of proteins related to immune function and reduced the immune response in pigs 70-71. In human, many immune response related genes are down-regulated in low birth weight newborns compared to normal birth weight newborns 72. Interestingly, some genes in our study, including RPL34, RPS17 and RPL9 (Additional File 2: Table S1), were the same down-regulated genes previously reported in human 72, suggesting that these genes might be important in antiviral innate immune response affected by a maternal viral infection in fetal development. Given their known developmental role in Rab GTPase activator activity, metal ion binding, and zinc ion binding (Additional File 2: Table S3), additional candidate genes (TBC1D19, WDFY3, and ADH4) were also identified in this study and suggest that these genes might be involved in embryonic or fetal development in pigs, respectively. These data provide fodder for future functional research that may aid in the discovery of genetic mechanisms of embryonic or fetal development in pigs.

Figure 3

Gene ontology analysis related to the candidate genes from DSRs between high and low estimated breeding value populations. These 183 functional genes harbored in 27 DSRs for IBW in piglets were identified to represent a total of 264 molecular functions or 706 biological processes. Here the top 10 molecular functions and the top 10 biological processes are shown by brown and blue bars, respectively. Vertical axis represents gene ontology categories, while horizontal axis indicates the number of genes in each ontology category.

In summary, we selected two groups of animals with high and low estimated breeding values (HEBV and LEBV) for IBW from a reference population of 522 piglets and revealed a total of 27 DSRs between them using Illumina PorcineSNP60 BeadChips. Twelve of these DSRs were consistent with QTL regions identified for piglet birth weight, average birth weight and litter birth weight reported previously 22. Moreover, several putative candidate genes and pathways in these regions based on Sscrofa 10.2 assembly were identified. The innovative feature of this study is that it provides an economic and effective way to confidently determine important genes and pathways related to economically important traits in pigs. More importantly, this study provides a novel genome-wide view of the genetic foundation for IBW in pigs.

Figure S1.

Click here for additional data file.

Table S1-S4.

Click here for additional data file.