Genome resequencing and bioinformatic analysis of SNP containing candidate genes in the autoimmune vitiligo Smyth line chicken model

Background The Smyth line (SL) chicken is the only animal model for autoimmune vitiligo that spontaneously displays all clinical and biological manifestations of the human disorder. To understand the genetic components underlying the susceptibility to develop SL vitiligo (SLV), whole genome resequencing analysis was performed in SLV chickens compared with non-vitiliginous parental Brown line (BL) chickens, which maintain a very low incidence rate of vitiligo. Results Illumina sequencing technology and reference based assembly on Red Jungle Fowl genome sequences were used. Results of genome resequencing of pooled DNA of each 10 BL and SL chickens reached 5.1x and 7.0x coverage, respectively. The total number of SNPs was 4.8 and 5.5 million in BL and SL genome, respectively. Through a series of filtering processes, a total of ~1 million unique SNPs were found in the SL alone. Eventually of the 156 reliable marker SNPs, which can induce non-synonymous-, frameshift-, nonsense-, and no-start mutations in amino acid sequences in proteins, 139 genes were chosen for further analysis. Of these, 14 randomly chosen SNPs were examined for SNP verification by PCR and Sanger sequencing to detect SNP positions in 20 BL and 70 SL chickens. The results of the analysis of the 14 SNPs clearly showed differential frequencies of nucleotide bases in the SNP positions between BL and SL chickens. Bioinformatic analysis showed that the 156 most reliable marker SNPs included genes involved in dermatological diseases/conditions such as ADAMTS13, ASPM, ATP6V0A2, BRCA2, COL12A1, GRM5, LRP2, OBSCN, PLAU, RNF168, STAB2, and XIRP1. Intermolecular gene network analysis revealed that candidate genes identified in SLV play a role in networks centered on protein kinases (MAPK, ERK1/2, PKC, PRKDC), phosphatase (PPP1CA), ubiquitinylation (UBC) and amyloid production (APP). Conclusions Various potential genetic markers showing amino acid changes and potential roles in vitiligo development were identified in the SLV chicken through genome resequencing. The genetic markers and bioinformatic interpretations of amino acid mutations found in SLV chickens may provide insight into the genetic component responsible for the onset and the progression of autoimmune vitiligo and serve as valuable markers to develop diagnostic tools to detect vitiligo susceptibility. Electronic supplementary material The online version of this article (doi:10.1186/1471-2164-15-707) contains supplementary material, which is available to authorized users.


Background
Vitiligo is an acquired hypomelanotic disorder characterized by circumscribed depigmented macules in the skin resulting from the loss of melanocytes. Autoimmunity has been identified as the major etiological factor in vitiligo, although many other factors including infections, stress, neural abnormalities, aberrant melanocyte function, and genetic susceptibility have been implicated [1]. The Smyth line (SL) chicken is the only animal model for autoimmune vitiligo that spontaneously displays all clinical and biological manifestations of the human disorder [2,3]. Like other autoimmune diseases, SL vitiligo (SLV) is multi-factorial in nature and involves the interplay of genetic, immune system, and environmental-factors. SLV susceptibility is manifested in part in an inherent melanocyte defect and loss of melanocytes is due to melanocyte-specific cell-mediated and humoral immune activities [2,3].
Recent genome wide association studies (GWAS) in humans to understand the role of genetic components in a variety of autoimmune diseases including vitiligo have identified hundreds of loci harboring risk alleles [4]. Several GWAS results identified vitiligo susceptible loci in human populations [5][6][7][8][9][10]. However, most susceptible loci identified by GWAS results were found in regulatory regions of gene expression, therefore the identified associations were not sufficient to identify the causal gene or deduce alterations caused by risk variants, which generally do not induce profound changes to genes (e.g. coding sequence changes, deletions, or duplications). Recently, the encyclopedia of DNA elements (ENCODE) of mammalian species suggested that~90% of disease associated genetic variations in human lie in noncoding regions, while only~10% of variations in coding regions were causative mutations associated with human disease [11,12]. Nevertheless, the identification of potential coding mutations that alter protein functionalities is a prerequisite process to understand disease etiology. Moreover, the functional study of candidate genetic risk factors is almost impossible without appropriate model systems.
The SL chicken is an excellent model to conduct a functional verification study of candidate genes that underlie genetic susceptibility for vitiligo due to the tractable, definite phenotype, the high vitiligo incidence in the population (80-90%), the feasibility of in vivo characterization and the relatively short generation time. Recently, microarray analysis examined global gene expression during SLV development and provided comprehensive information at the transcriptome level that supported the multifactorial etiology of vitiligo [13]. In this study, whole genome resequencing analysis using an Illumina platform was performed to more deeply investigate the genetic aspects of SLV expression in comparison with the parental Brown line (BL) of chickens from which the SL was originally derived. BL chickens retain vitiligo susceptibility but with a very low (0 -2%) incidence rate of vitiligo development [2,3], although none of the BL chickens used in this study had vitiligo. Millions of single nucleotide polymorphisms (SNP) were identified by genome resequencing and only potentially causal genes containing non-synonymous mutations that can induce amino acid changes in proteins were focused on in this study.

Results and discussion
Genome resequencing for BL and SL chickens Genome sequencing of pooled DNA from 10 nonvitiliginous BL and SL chickens each with confirmed SLV produced~63 and 89 million sequence reads of 200 bp, respectively (Table 1). Of those,~80% of the reads were used for sequence alignment, while 20% of sequence reads were not aligned. Therefore, genome coverage for BL and SL reached 5.1x and 7.0x, respectively, of the Red Jungle Fowl chicken genome. The total number of SNPs was 4.8 and 5.5 million (~0.5% of template genome) for the BL and SL genome, respectively. The large number of SNPs per examined chicken line was based on data of at least 2 read coverage depths (number of read counts per nucleotide location). Most SNPs were found on the larger chromosomes (Chr), including Chr 1 through 5 and Z (sex chromosome) (data not shown). To identify genetic biomarkers that are responsible for the incidence of SLV, unique SNPs that are found in SL only were selected by removing SNPs that overlapped with those found in the parental BL. Then, mutations with ≥75% SNP rates were chosen as reliable marker SNPs. Since the objective of this study was to identify mutational SNPs uniquely found in SL compared to the parental BL, the filtering process used did not involve a typical SNP calling and filtering method based on quality score (Q call column in Additional file 1: Table S1) [14]. Instead, SNP filtering was conducted by removing overlapping SNPs found in both BL and SL, and applying fixed %SNP rates (≥75%) as described in Methods. As a result, a total of~1 million unique SNPs were identified throughout the SL chicken genome ( Figure 1). Over 100,000 SNPs were found in Chr 1, 2, 3, and Z while Chr 32 did not contain any unique SNPs for SL. When SNPs were grouped by the feature types of chromosome regions,~50% of SNPs were in the intergenic (heterochromatic) regions and 13,710 SNPs were found in CDS sequences (protein coding regions) ( Figure 2). Most genes containing SNPs in regulatory regions, not in CDS regions, identified by human GWAS studies were also observed to contain unique SNPs in the current SL study (data not shown). Around 60% of SL SNPs in protein coding CDS regions were synonymous mutations that did not induce amino acid changes. To identify potentially causal mutations that induce protein coding alterations, SNP analysis focused on SNPs leading to changes in amino acid sequences. Using this approach, a total of 3518 SNPs were identified that could induce non-synonymous-, frameshift-, nonsense-, and no-start-changes in the CDS region ( Figure 2 and Additional file 1: Table S1), suggesting that the 3518 SNPs are part of the genetic components in functional protein coding regions that may drive the high incidence rate of SLV. Of the 3518 candidate SNPs that are associated with amino acid changes, SNPs showing ≥10 read depths (considered to be more reliable candidate genetic markers) were chosen for the further analysis. Using this approach, 195 SNPs remained (data not shown). To reduce false positives due to possible errors in the assembly process, re-scanning of each SNP position for the 195 potentially more reliable protein coding SNPs was conducted using the Seqman-Pro viewer program. This process yielded 156 more reliable SNPs that were chosen as candidate marker SNPs for further analysis (Table 2).

SNP validation using PCR and Sanger sequencing
Since pooled DNA samples of 10 chickens for each line were used for genome sequencing, individual SNPs were          subjected to the verification process with larger bird populations. For this, 14 SNPs were randomly chosen from the 156 candidate marker SNPs showing ≥10 read depths and were subjected to SNP verification analysis using PCR and Sanger sequencing to detect SNP positions with larger numbers of birds; specifically 20 nonvitiliginous BL chickens and 70 SL chickens exhibiting vitiligo. The results clearly showed differential frequencies of nucleotide bases in the 14 SNP positions between BL and SL chickens (Table 3). Thus, the 156 SNPs known to induce amino acid changes can become potential genetic biomarkers for vitiligo in SL chickens.

Bioinformatic analyses of genes containing amino acid change SNPs
Amino acid changes may have impacts on the functional interpretations for vitiligo induction in SL chickens. The Ingenuity Pathway Analysis (IPA) program generated bioinformatics data sets including functional groups (gene ontology; GO) and gene networks for genes containing amino acid changes in SL chicken.
The 156 SNPs were found in 139 genes encompassing known-and unknown functions, chromosomal open reading frames, and hypothetical proteins (Additional file 2: Table S2).

Functional roles
Genes were categorized in 76 functional groups (Additional file 3: Table S3). Of these, six functional groups are of particular interest to autoimmune vitiligo development, including dermatological diseases/conditions, inflammatory response, inflammatory disease, immunological disease, immune cell trafficking, and infectious disease ( Table 4). The functional group of genes for dermatological diseases/conditions contained the following genes: ADAMTS13 (ADAM metallopeptidase with thrombospondin type 1 motif 13); ASPM [asp (abnormal spindle) homolog, microcephaly associated; Drosophila)]; ATP6V0A2     -Is a RhoGEF protein that interacts with cytoskeletal calmodulin and titin and is part of the giant sarcomeric signaling protein family of myosin light chain kinases -Mutant human OBSCN protein (E4574K) is associated with melanoma in human.
-Somatic missense homozygous mutant human STAB2 gene (c.3862 T > G translating to p.S1288A) is associated with melanoma in skin from human leg (observed in 2 of 2 samples) LRP2 (or megalin) [17,18] -Is a member of the low density lipoprotein (LDL) receptor gene family essential for brain development.
-Somatic missense heterozygous mutant human LRP2 gene (c.6284G > A translating to p.R2095Q) is associated with melanoma in skin from human leg (observed in 2 of 2 samples).
ASPM [17,19] -This gene is the human ortholog of the Drosophila melanogaster 'abnormal spindle' gene (asp), which is essential for normal mitotic spindle function in embryonic neuroblasts.
-Somatic nonsense heterozygous mutant human ASPM gene (c.7174C > T translating to p.Q2392*) is associated with melanoma in skin from human leg (observed in 2 of 2 samples).

GRM5 [20]
-Is a member of G protein-coupled receptor that are widely expressed in the brain and modulate many diverse signaling pathways -In mouse melanocytes, transgenic rat GRM5 protein (S901A mutant) affects development of melanoma in mouse.
BRCA2 [21,22] -At the cellular level, loss of BRCA2 function results in sensitivity to cross-linking agents, a decrease in homologydirected repair of double-stranded DNA breaks (DSBs), and defects in replication and checkpoint control -Inherited mutations in BRCA1 and this gene, BRCA2, confer increased lifetime risk of developing breast or ovarian cancer.
-Mutant human BRCA2 gene is associated with malignant melanoma in Homo sapiens.
XIRP1 [17,23] -Its function is unknown, but it is upregulated in wounded skeletal muscle cells in zebrafish -Somatic nonsense heterozygous mutant human XIRP1 gene (c.2838G > A translating to p.W946*) is associated with melanoma in skin from human leg (observed in 2 of 2 samples).

RIDDLE syndrome RNF168 [24]
-Is an E3 ubiquitin ligase -Mutant human RNF168 gene (deletion c.1323_1326del of ACAA and DNA duplication mutation) is associated with RIDDLE (radiosensitivity, immunodeficiency, dysmorphic features, and learning difficulties) syndrome, which a novel human immunodeficiency disorder associated with defective double strand break repair.
-Is associated with the development of thrombotic thrombocytopenic purpura (TTP), known as the Schulman-Upshaw syndrome.
Autosomal recessive cutis laxa type IIA ATP6V0A2 [26] -Is a subunit of the vacuolar ATPase (v-ATPase), a heteromultimeric enzyme that is essential for the acidification of diverse cellular components.
Interestingly, a recent report by Nikolaev et al. (2012) [17] indicated that amino acid changes found in ASPM,  -Is a serine protease involved in degradation of the extracellular matrix and possibly tumor cell migration and proliferation.
-Heterozygous-and heterozygous mutant mouse Plg gene in mouse affects development of blister in mouse subepidermal skin that is altered by transgenic uPAR (PLAUR) protein and transgenic uPA (product of PLAU) protein and development of blister.
LRP2, STAB2, and XIRP1 proteins were associated with human melanoma by exome sequencing. Melanocytes in vitiligo also exhibit morphological and biological melanocyte defects/alterations compared to melanocytes from individuals with normal pigmentation [29]. While these alterations may be different from those observed in melanoma, e.g. slower growth, and higher sensitivity to oxidative stress of cultured melanocytes [30,31], alterations in amino acid sequences found in homolog proteins but different residues may result in opposite phenotypes of dermatological diseases/conditions [32]. In addition to these molecules, BRCA2, GRM5, MKI67, and OBSCN associated with SLV are also known to be associated with human melanoma. The relationship between candidate genes and other dermatological diseases including melanoma is summarized in Table 5.

Gene networks
Gene network analysis, which represents the intermolecular connections among interacting genes based on functional knowledge inputs, was performed on genes with amino acid changes in SLV chickens using the IPA program. The gene network analysis was carried out using the simplest setting of 35 focus molecules to facilitate and summarize the intermolecular connections (Table 6 and Figures 3, 4, 5, 6, 7, 8, and 9). A discussion of the 7 gene networks is provided below and gene information for focus molecules in each network is listed in Additional file 4: Table S4. Candidate genes in Network #1 are associated with signaling pathways of the mitogen activated protein kinase (MAPK; also ERK1/2) and protein kinase C (Pkc) connected to VEGF (vascular endothelial growth factor) and PDGF (platelet derived growth factor) with PLAU in the center (Figure 3). The top functions related to network #1 are cardiovascular disease, hematological disease, and cardiac infarction. Interestingly, molecules including LRP2, PLAU, ADAMTS13, and GRM5 that are part of Network #1 were also identified as functional factors for melanoma as described above [17]. Additionally, mutations in the amino acid sequence of LRP2, altered function of MAP2K1 and MAP2K2 induced by genetic mutations in melanoma patients [17], and mutations in GRM5 in mouse melanoma models [20] were also reported. The connections in Network #1 therefore suggest genetic mutations that generated amino acid changes in LRP2, PLAU, ADAMTS13, and GRM5 may influence dermatological diseases, including vitiligo, through MAPK and ERK1/2 signaling pathways.
The top functions of Network #2 include Cell Cycle, DNA Replication, Recombination and Repair, and Developmental Disorder (Figure 4) and Network #3 is involved in Developmental Disorder, Endocrine System Disorders, Hereditary Disorder ( Figure 5). Most molecules in Networks #2 and #3 directly bind to UBC (ubiquitin C). Ubiquitinylation, the covalent attachment of ubiquitin to proteins, regulates numerous cellular processes such as protein degradation and signal transduction. Recently, many ubiquitinylated proteins and their lysine ubiquitinylation site were identified using proteomic technologies in mammalian species [33][34][35][36][37]. Indeed, in SLV, the amino acid changes in lysine residues for CEP192 (centrosomal protein 192 kDa; pK169R) and API5 (apoptosis inhibitor 5; pK1219Q) were identified ( Table 2), suggesting that the various cellular functions including protein degradation by altered ubiquitinylation properties of proteins may play a significant role in the induction of vitiligo.
Molecules in Network #4 are involved in Cell Morphology, Cellular Function and Maintenance, Hair and Skin Development and Function ( Figure 6). Molecules in network #4 mainly interact with IFNG (interferon gamma), IL4 (interleukin 4), MAPK8, and calcium signaling pathways in addition to protein phosphatase (PPP1CA; protein phosphatase 1 catalytic subunit alpha isozyme) functions. MYO16 (myocin 16), KIF18A (kinesin family member 18A), and CASC1 (cancer susceptibility candidate 1) are known to directly bind to PPP1CA [38][39][40], suggesting that amino acid changes in those proteins may induce alterations (hyper vs hypo) in protein phophorylation states in SL chickens, resulting in vitiligo development. Molecules interacting with IFNG, IL4, MAPK8 and calcium signaling pathways showed an indirect relationship, not a direct relationship with each other making it difficult to explain how amino acid changes in these molecules affect vitiligo induction in SL chicken.
Molecules in Network #5 also mainly bind to UBC as discussed in Networks #2 and #3 and the functions include Lipid Metabolism, Small Molecule Biochemistry, Digestive System Development and Function (Figure 7).
Network #6 contains molecules involved in Hereditary Disorder, Ophthalmic Disease, Neurological Disease ( Figure 8). In this network, ZC2HC1A (zinc finger, C2HC-type containing 1A) and RUFY3 (RUN and FYVE domain containing 3) directly bind to APP [amyloide beta (A4) precursor protein]. APP is a precursor protein for beta-amyloide, which is the main constituent of amyloid plaques in the brains of Alzheimer disease patients [41]. RUFY3 (also known as single axon-related; singar1), which is a brain specific protein, regulates neuronal polarity by suppressing formation of surplus axons [42]. Though the binding of ZC2HC1A and RUFY3 to APP was found during the progression of Alzheimer disease [41], the functional roles for this binding in the progression of this disease has not been characterized. Similarly, amino acid changes found in ZC2HC1A and RUFY3 are implicated in SLV development possibly as a result of altered APP binding properties. USH2A [Usher syndrome 2A (autosomal recessive, mild)] was included in network #6. Various mutations in USH2A have been identified in patients of Usher syndrome type II, which is characterized by moderate to severe sensorineural hearing loss and progressive retinitis pigmentosa [43]. Vitiliginous SL chickens may also develop severe visual impairment and blindness due to autoimmune activity directed against choroidal melanocytes and subsequent damage to the retinal pigment epithelium. [44,45]. Taken together, the amino acid change in USH2A also may affect vitiligo progression and retinal depigmentation.
Molecules in Network #7 mainly function in Cellular Response to Therapeutics, Cellular Assembly and Organization, DNA Replication, Recombination, and Repair. PRKDC (protein kinase, DNA-activated, catalytic polypeptide) and its interacting protein kinases are mainly involved in this network ( Figure 9). In addition to knock-out and inactive mutations, alteration of autophosphorylation capability by single amino acid change of PRKDC has been known to influence rejoining of DNA double stranded breaks in mammalian cells [46], suggesting that vitiligo development may be affected by aberrant PRKDC kinase activity due to an observed amino acid change.
Vitiligo susceptible loci in human populations identified by several GWAS [5][6][7][8][9][10]  , and GZMB (granzyme B). One similar mutation in UBASH3A protein coding region from human [6] was also found in the SLV chicken model (Table 2). Additionally, when the long list of 3518 candidate amino acid altering SNPs (read depth <10) was considered, several genes, including IFIH1 (interferon-induced helicase C domaincontaining protein 1), CD44 antigen and DNAH5 (dynein, axonemal, heavy chain 5), matched those identified by human GWAS studies [Additional file 1: Table S1 and [5,6]], although the amino acid position and alterations did not match. UBASH3A is one of the two family members belonging to the T cell ubiquitin ligand (TULA) family and can negatively regulate T cell signaling [47]. Together with the UBC molecule discussed elsewhere in this paper, functions for UBASH3A related to ubiquitinylation and T cell signaling pathway may be important for SLV development. IFIH1 encodes an interferoninduced RNA helicase involved in antiviral innate immune responses, associated with type 1 diabetes, Graves' disease, Figure 6 Gene network #4. Molecular interaction and symbols are the same as the description in Figure 3. multiple sclerosis, psoriasis, and perhaps lupus [48][49][50][51][52][53]. CD44 encodes a cell surface glycoprotein with various functions, including a role in T cell development [54], and is associated with lupus [55]. DNAH5 gene mutation is found in patients with primary ciliary dyskinesia (PCD) [56], a rare disease transmitted as an autosomal recessive trait and characterized by recurrent airway infections due to abnormal ciliary structure and function. Primary defects in the structure and function of sensory and motile cilia result in multiple ciliopathies [57].

Conclusions
In this study, various potential genetic markers showing amino acid changes were identified in the SLV model through genome re-sequencing. When considering functionality based on the interpretation of factors involved, development of vitiligo appeared to be associated with the interactions among cytoskeletal factors (OBSCN, ASPM, XIRP1, ADAMTS13), protein kinases (MAPK, ERK1/2, PKC, PRKDC), phosphatase (PPP1CA), ubiquitinylation (UBC) and amyloid (APP) production. Further functional validation study, such as allele specific expression of the candidate genes with candidate SNPs at the target tissues involved in SLV development will be carried out using the SL chicken model for spontaneous autoimmune vitiligo.

Animals and Illumina sequencing
Adult SL chickens with vitiligo and parental nonvitiliginous BL chickens, maintained by G. Erf at the University of Arkansas (Fayetteville, AR), were selected from the breeder populations. Blood (3 ml) was collected from 12 birds each following an animal use protocol approved by the University of Arkansas Institutional Animal Care and Use Committee (IACUC; approval number: 11019). Genomic DNA was isolated from each whole blood sample using QiaAmp DNA mini kit (Qiagen, Hilden, Germany) following manufacturer's instructions. DNA quality was determined by agarose gel electrophoresis and 10 samples having the highest quality in each line were pooled to represent each line. Library preparation and Illumina genome sequencing for the pooled DNA samples were performed by the National Center for Genome Resources (NCGR; Santa Fe, NM). Illumina HiSeq system 2x100 bp paired end read technology was used for genome sequencing.

Genome sequence assembly and data analysis
Illumina sequencing data received from NCGR was aligned to the chicken reference genome sequence for Red Jungle Fowl (GBK 4.0) that was retrieved from NCBI. For the reference based genome alignment, the NGen genome sequence assembly program of the Lasergene software package (DNAStar, Madison, WI) was used. Assembly parameters were as follows: file format, BAM; mer Size, 21; mer skip query, 2; minimum match percentage, 93; maximum gap size, 6; minimum aligned length, 35; match score, 10; mismatch penalty, 20; gap penalty, 30; SNP calculation method, diploid bayesian; minimum SNP percentage, 5; SNP confidence threshold, 10; minimum SNP count, 2; minimum base quality score, 5. After assembly, the SeqMan Pro program of the Lasergene package was used for further analyses including SNP data.

SNP detection and analysis
JMP genomics (SAS Institute, Inc., Cary, NC) program was used for filtering unique SNPs for vitiligo SL chickens. SNPs occurring in both SL and BL lines were filtered out, leaving behind unique SNPs for each line. To identify highly fixed and homozygous SNPs, the SNPs were filtered based on SNP percentages (SNP%). SNPs with a SNP% of ≥75 (%) (for example, number of SNP = 3 of read depth = 4) were chosen. The 75% cutoff for SNP selection was set by considering potential sequencing errors that can be generated by the massively parallel sequencing method. Potential causal SL SNPs that induce non-synonymous changes in CDS regions were chosen for further analysis. Since the read depth of many SL SNPs was low, unique SNPs showing ≥10 read depths were considered as reliable SNPs. Reliable and causal SNPs, which were chosen by criteria described above were confirmed by double-checking the raw assembly data with alignment view to reduce false positives.

SNP validation using PCR and Sanger sequencing
Fourteen randomly chosen SNPs, which induce amino acid changes in the CDS region, were subjected to validation using PCR and Sanger sequencing with larger numbers of SL and BL chickens. Twenty BL and 70 SL chickens that were verified phenotypically to be nonvitiliginous and vitiliginous, respectively, were used for blood sampling. Approximately 100 μL of blood was collected from each bird by wing vein puncture into tubes containing citrate (anticoagulant). Genomic DNA was isolated from whole blood using the Wizard SV 96 Genomic DNA Purification System (Promega; Madison, WI) following manufacturer's instructions with modifications. Isolated DNA was quantified using a Nanodrop 1000 spectrophotometer (Thermo Fisher Scientific Inc., Waltham, MA) and a dilution of 1 ng/ μL was prepared in 96 well PCR format for all samples. For PCR reaction, forward and reverse primers were designed based on the RJF genome sequence (GenBank assembly ID: GCA_000002315.2) using Primer 3 online software ( Table 7). The sequencing primers were designed to anneal at least 50 bp upstream of the SNP position and forward/reverse primers were chosen at the flanking regions of the seq primer and the SNP position. All primers were commercially synthesized by Integrated DNA Technology (Ames, IA). PCR was carried out as 25 μL reaction volumes in 96 well plates with cyle conditions as follows: denaturation 95°C for 1 min, 40 cycles of amplification (95°C for 30 s, 55-63°C for 1 min, 72°C for 1 min), final extension 72°C for 10 min. Verification of PCR was performed by 1% agarose gel electrophoresis. PCR products were purified using the Wizard SV 96 PCR Clean-Up System (Promega; Madison, WI) following manufacturer's instructions. Briefly, four plates (four different PCR products) were pooled into one plate and were subjected to PCR clean-up. Cross-specificity of seq-primers to the four pooled PCR products was examined by BLAST function (NCBI) and only products that were not cross- specific with other seq primers were pooled. Purified PCR products were subjected to Sanger sequencing performed by the University of Arkansas DNA Resources Center (Fayetteville, AR). Results were analyzed using ABI Sequence scanner software (Life Technologies, Carlsbad, CA). Ratios of bases occurring at SNP locations were recorded.

Bioinformatics
Functional interpretation of 139 genes showing ≥75 SNP%, ≥10 read depths and non-synonymous changes was analyzed in the context of gene ontology and molecular networks using Ingenuity Pathways Analysis (IPA; Ingenuity Systems®; www.ingenuity.com). Since IPA is based on human and mouse bioinformatics, functionalities for differentially expressed genes in the chicken were interpreted based primarily on mammalian biological mechanisms. The limit of number of molecules in the network was set to 35, leaving only the most important molecules based on the number of connections for each focus gene (a subset of uploaded significant genes having direct interactions with other genes in the database) to other significant genes [58].