The Identification of Loci for Immune Traits in Chickens Using a Genome-Wide Association Study

The genetic improvement of disease resistance in poultry continues to be a challenge. To identify candidate genes and loci responsible for these traits, genome-wide association studies using the chicken 60k high density single nucleotide polymorphism (SNP) array for six immune traits, total serum immunoglobulin Y (IgY) level, numbers of, and the ratio of heterophils and lymphocytes, and antibody responses against Avian Influenza Virus (AIV) and Sheep Red Blood Cell (SRBC), were performed. RT-qPCR was used to quantify the relative expression of the identified candidate genes. Nine significantly associated SNPs (P < 2.81E-06) and 30 SNPs reaching the suggestively significant level (P < 5.62E-05) were identified. Five of the 10 SNPs that were suggestively associated with the antibody response to SRBC were located within or close to previously reported QTL regions. Fifteen SNPs reached a suggestive significance level for AIV antibody titer and seven were found on the sex chromosome Z. Seven suggestive markers involving five different SNPs were identified for the numbers of heterophils and lymphocytes, and the heterophil/lymphocyte ratio. Nine significant SNPs, all on chromosome 16, were significantly associated with serum total IgY concentration, and the five most significant were located within a narrow region spanning 6.4kb to 253.4kb (P = 1.20E-14 to 5.33E-08). After testing expression of five candidate genes (IL4I1, CD1b, GNB2L1, TRIM27 and ZNF692) located in this region, changes in IL4I1, CD1b transcripts were consistent with the concentrations of IgY, while abundances of TRIM27 and ZNF692 showed reciprocal changes to those of IgY concentrations. This study has revealed 39 SNPs associated with six immune traits (total serum IgY level, numbers of, and the ratio of heterophils and lymphocytes, and antibody responses against AIV and SRBC) in Beijing-You chickens. The narrow region spanning 247kb on chromosome 16 is an important QTL for serum total IgY concentration. Five candidate genes related to IgY level validated here are novel and may play critical roles in the modulation of immune responses. Potentially useful candidate SNPs for marker-assisted selection for disease resistance are identified. It is highly likely that these candidate genes play roles in various aspects of the immune response in chickens.


Introduction
Great effort has been expended globally to understand and genetically improve disease resistance in domestic animals [1][2][3][4]. Immune capacity associated with specific diseases may be useful indicators for indirect selection for general disease resistance, because such traits can be evaluated and quantified in live animals [5][6][7][8][9].
More recently, genome wide association studies (GWASs) have become one of the most commonly used strategies for identifying genes for complex traits in humans, as well as in animals [22]. In chickens, some major loci associated with resistance to Marek's disease [23] and immune response to Newcastle disease virus and infectious bronchitis virus were identified by GWASs [9,22]. Despite these studies, there is still limited information about the multiple immune traits that underly the full immune response at the genome-wide level in chickens.
This study aimed to identify major genomic regions (loci) and candidate genes associated with the immune response using GWAS for an array of important immune traits, including total serum concentrations of immunoglobulin Y (IgY), numbers of, and the ratio of heterophils and lymphocytes (H/L), and antibody responses against Avian Influenza Virus (AIV) and Sheep Red Blood Cells (SRBC) in chickens. This approach may offer valuable information for understanding the molecular mechanisms of the regulation of immune traits and facilitate the application of marker-assisted selection in breeding programs for disease resistance in chickens.

Results and Discussion
Genome-wide association analysis A total of 7,175 independent SNP markers were obtained with multidimensional scaling (MDS) analysis of these SNPs using the first two principal components ( Fig. 1) indicating that chickens within each half-sibling family were clustered together. No obvious population substructure was detected. The analytical method was as described previously [24] and resulted in the elimination of approximately 13,506 SNPs from further analysis. The distribution of SNPs after this quality control is summarized in Table 1. All raw genotypes are available from the Dryad Repository (doi:10.5061/dryad.8m4r5).
The descriptive statistics for six immune traits in the 720 Beijing-You chickens used for the present GWAS are shown in Table 2. The numbers of genome-wide significant SNPs detected for the six immune traits, and the details of these significant SNPs, including their positions in the genome, the nearest known genes and the P-values, are given in Table 3 and Table 4. The global view of P-values (in terms of-log 10 (P-value)) for all SNP markers of each trait is represented by a Manhattan plot shown in Fig. 2. The genome-wide association analysis in this study identified 39 SNPs, consisting of nine significantly associated SNPs and 30 suggestively associated SNPs.

Serum Concentrations of IgY at d 80
Five SNPs, associated with the serum total IgY level, were clustered within a 0.26 Mb region on chromosome 16 (S1 Table). One SNP, rs14738106, was located 2.9 kb downstream of the B-G gene (within the MHC region). Four of the SNPs (rs16057130, rs14738106, rs15788101 and rs15788237) were located within CD1b and CD1c genes, an intron of GNB2L1 gene, an intron of IL4I1 gene and an intron of BMA1 gene, respectively. A QTL for titres of Natural antibodies (NA) was detected at the SNP rs15788216 on chromosome 16 in chickens [14][15], near the SNP rs15788237 found here.
The chicken MHC, which is composed of the B and Y complexes [25] and plays a central role in the immune response, has also been mapped to this NOR (Nucleolus Organizing Region, cluster of genes for ribosomal RNA 18S, 5.8S and 28S)-bearing microchromosome 16 [26].
CD1b and CD1c encode two members of the CD1 family, which are structurally similar to the MHC class I antigen and have been shown to present lipid antigens for recognition by T lymphocytes [27]. The expression of the IL4I1 gene can be induced by interleukin 4 in B cells. The IL4I1 protein has L-phenylalanine oxidase activity, and oxidative deamination of phenylalanine by IL4I1 produces H 2 O 2 , and functions in the innate immune defenses of vertebrates. Previous studies have revealed that the IL4I1 protein inhibits T-cell proliferation with a preference for memory T cells, and H 2 O 2 mediates this inhibition [28][29].
In the current study, significant SNPs on GGA1, 3, 5 and 12 were also identified that related to serum total IgY level. Those were not located within the previously known QTLs for immune traits, for example, five QTLs associated with humoral innate immune response in chickens (NA titers binding LPS, located on GGA5, GGA9, GGA14, GGA18, and GGAZ) were recently validated in a cross between Green-Legged Partridge-like and White Leghorn [19,21]. The reason for the discrepancy might result from the different chicken populations being studied.

Heterophils and Lymphocytes
The heterophil to lymphocyte (H/L) ratio in the peripheral blood of chickens has been widely accepted as a reliable and accurate physiological indicator of the stress response in chickens [30][31]. The number of lymphocytes decrease and the number of heterophils increase in response to stressors such as high temperature, hypothermy, excessive NH 3 , bacterial infection, inter alia [32][33]. The ratio, and individual numbers of heterophils and lymphocytes are highly heritable traits, with heritability estimates exceeding 0.5 [34], suggesting that good responses to selection should be possible. There were three SNPs detected as being suggestively associated with the H/L ratio (S2 Table), two of which were simultaneously associated with the number of heterophils (S3 Table). One SNP (rs14562554) was mapped to 559kb downstream of the growth hormone inducible transmembrane protein (GHITM) gene on chromosome 6. The GHITM protein is part of the interferon signaling system that includes Janus kinases and their downstream target Gene name alone indicates the SNP is within the gene. c D indicates that the SNP is downstream of the gene. d U indicates that the SNP is upstream of the gene. These definitions are also used in Table 4. STAT proteins [35]. One SNP at 7.2Mb on chromosome 13 is located 192kb downstream of gamma-aminobutyric acid receptor subunit alpha-6 (GABRA6). This gene is associated with increased production of cortisol and increased blood pressure in response to psychological stress [36]. The rs16743572 SNP is located downstream of DnaJ homolog, subfamily C, member 16 (DNAJC16) in chromosome 21; the gene is a member of the DnaJ/Hsp40 (heat shock protein 40) family. These proteins are molecular chaperones induced in response to cellular stress [37] and have been preserved throughout evolution, functioning in protein translation, folding, unfolding, translocation, and degradation, primarily by stimulating the ATPase activity of Hsp70s chaperone proteins [38]. Human Hsp40s have an immune-modulatory effect on T cells and cytokine secretion in rheumatoid arthritis patients [39].
There were two SNPs, rs15950814 and GGaluGA134321, that were suggestively associated with the number of lymphocytes (S4 Table). The former is located within phosphotyrosine interaction domain containing 1 (PID1), also called NYGGF4, on chromosome 9, which encodes a protein possessing a phosphotyrosine interaction domain (PID). Mouse and human NYGGF4 are involved in the regulation of proliferation of adipocyte precursors [40]. The gene is more widely expressed in chicken tissues than in those of mice and humans [41]. Proteins with PID are widely involved in physiological processes, such as neural development, immune function, homeostasis and cell growth [42], which indicates that the chicken NYGGF4 homolog possibly possesses multiple functions. The GGaluGA134321 SNP, is located 11kb downstream of the zinc finger E-box binding homeobox 1 (ZEB1) gene on chromosome 2. In human and mouse studies, ZEB1 is a transcriptional repressor and plays an important role in regulating lymphoid differentiation [43]. It has been implicated in the regulation of other important genes in the immune system such as the Ig heavy chain enhancer, CD4, GATA-3, α7-integrin and IL-2 [44][45]. It has been shown that ZEB1-deficient mice have poorly developed thymuses and significant decreases in their T cell population [46].

Antibody responses to SRBC
The direct response to selection on antibody titer to SRBC challenge is accompanied by correlated responses in production and disease-related traits [13,16]. The heritability of anti-SRBC titer in chickens is moderate (0.20 to 0.25), and the trait responds to divergent selection [47].
In this study, of special interest, five SNPs in a narrow region (9.16Mb-9.31Mb) on chromosome 19 were suggestively associated with antibody titer against SRBC (S5 Table). The rs10725789 and GGaluGA128939 are located within and 48 kb downstream of the inducible nitric oxide synthase (iNOS) gene. The expression of iNOS is regulated by cytokines such as IFN-γ, IL-1β, IL-6, TNF-α, and iNOS-derived nitric oxide (NO) is well known to play important roles in infection and immunity [48]. The rs15050742, rs13576070 and GGaluGA128957 SNPs are located 21kb downstream, 1kb and 8kb upstream of the tumor necrosis factor, alphainduced protein 1 (TNFAIP1). The TNFAIP1 protein is a TNF-α and IL-6 inducible protein.
Research in humans and rats indicate that the protein may have roles in DNA synthesis and repair and in apoptosis [49]. Such results refine the detail of this genomic region and lead to more conclusive identification of the contributing loci involved and the putative causal mutations.
In addition, three SNPs on chromosome 1 were suggestively associated with anti-SRBC titer (S5 Table), and rs13905265 is 68 kb downstream of G protein-coupled receptor 89B (GPR89B), but no known genes were found within 500 kb of the other two SNPs. A QTL has been detected in the vicinity of this region in a F2 cross of two lines divergently selected for SRBC antibody responses [16]. Two SNPs of suggestive significance at 13.9 Mb and 14.6 Mb are consistent with a genomic region identified in a previous study in chickens [50]. There are several other genomic regions previously identified as QTLs for SRBC antibody response in chickens, including those on GGA2, GGA3, GGA5, GGA11, GGA16, GGA23 and GGA24 [13, 16-19, 21, 50-51].

Antibody responses to AIV
Two SNPs on chromosome 1 were suggestively associated with AIV antibody titer (S6 Table). One SNP, rs13937877 is located 234 kb downstream of solute carrier family 25, member 6 (SLC25A6). The present study found that the region at about 100 Mb from the proximal end of chicken chromosome 1, including the ROBO1 (Roundabout, axon guidance receptor, homolog 1) and ROBO2 (homolog 2) genes, has a strong effect on the antibody response to the NDV in chickens [22]. One SNP rs14527240 was located in the SPARC-related modular calciumbinding 1 (SMOC1) gene on GGA5 playing an important role in MD tumorigenesis [23]. GGA1 also had the most significantly associated SNPs on the antibody level against infectious bronchitis virus (IBV), but a genomic region containing a cluster of 13 beta-defensin (GAL1-13) and interleukin-17F genes on GGA3 probably plays an important role in the immune response against IBV [9]. The SLC25A6 gene encodes adenine nucleotide translocase 3 (ANT3), an enzyme which exchanges ATP and ADP through the mitochondrial membrane. Mitochondrial ANT3 is regulated by IL-4 and IFN-γ via STAT-dependent pathways [52]. Considering the critical role of mitochondrial ANTs in energy metabolism and apoptosis, ANT3 regulation by IL-4 and IFN-γ may have a functional implication in cytokine-mediated T cell survival [53].
One SNP (rs15995401) is located 28 kb downstream of hematopoietic cell-specific Lyn substrate 1 (HCLS1). HCLS1 is an intracellular protein expressed mainly in hematopoietic cells. It was reported to be an important actin regulator at the T-cell immunologic synapse and also to influence numerous functions of natural killer cells [54].
Eight SNPs on chromosome Z were suggestively novel and specific associated with AIV antibody titer in this study (S6 Table). One SNP (rs14778781) is located within multiple PDZ domain protein (MPDZ) which is involved in a diverse range of pathways that maintain cell integrity and function [55]. No known genes were found within 500kb of the remaining seven SNPs. These studies concerning immune responses against NDV, IBV and AIV add to basic knowledge of host attributes potentially relevant to improved methods for controlling outbreaks of viral diseases.

Relative expression of candidate genes in spleen
The present GWAS verified that the important immune-related locus was located on chromosome 16 [9,22]. In this study, five of the significant SNPs for total serum IgY level were located within a narrow 247 kb region on chromosome 16. On this basis, the genes within this region were selected for further expression analysis, including CD1b, IL4I1 (Interleukin 4-induced 1), TRIM27 (tripartite motif-containing 27), GNB2L1 (guanine nucleotide binding protein, beta polypeptide 2-like 1), and ZNF692 (zinc finger protein 692).
Compared with SRBC, Salmonella enteritidis is more potent in inducing production of serum antibodies (IgY) in the infected chickens. For just that reason, it was used here to validate the association between total serum IgY level and the corresponding candidate genes identified in this study.
The relative expression of those five candidate genes near associated signals from the GWAS analysis was determined by RT-qPCR in chickens with and without challenge with Salmonella enteritidis. In this study, the changes across days in expression of IL4I1 and CD1b followed a similar pattern to the concentrations of IgY, while changes in TRIM27 (alternatively named RET finger protein, RFP) and ZNF692 expression were the reverse (mirror-image) of concentrations of IgY (Fig. 3).
Chickens possess two CD1 genes designated chCD1-1 and chCD1-2 that are closely linked to the MHC locus on chromosome 16. Both genes are expressed in lymphoid tissues as evidenced by Northern analysis and PCR [56]. The spleen was examined here because it is the major lymphoid organ in birds and has previously been shown to have high level expression of CD1 [57].
IL4I1 modulates the number of antitumor cytotoxic T cells in vivo and affects their capacity to secrete IFN-γ. The expression of IL4I1 is restricted to lymphoid tissues, with the highest levels found in lymph nodes and spleen [58]. In the present study, the expression of CD1b and IL4I1 was up-regulated in spleens of infected birds.
TRIM27 is a member of the TRIM superfamily and exhibits the classical PRY-SPRY domain C-terminal of the RBCC motif. TRIM27 is highly expressed in mouse spleen and thymus [59] and in hematopoietic cells [60]. Roles for TRIM27 in transcriptional repression [60], the negative regulation of NF-κB and in IFN-signalling pathways, apoptosis and in cell cycle regulation have been reported [61], indicating its involvement the control of multiple cellular processes. These findings account for the present demonstration that the abundance of TRIM27 showed reciprocal changes to those of IgY concentrations.
The genes ZNF692 and GNB2L1 relate to the MHC genes, but little is known of their precise function in chickens. The changes in relative expression and relationship to changes in IgY after infection, support their having some role, probably negative regulatory, in the magnitude of the immune response.

Conclusions
This study has revealed 39 SNPs associated with six immune traits (total serum IgY level, numbers of, and the ratio of heterophils and lymphocytes, and antibody responses against AIV and SRBC) in Beijing-You chickens. Five of the significant SNPs for total serum IgY level were located within a narrow 247 kb region on chromosome 16. Some SNPs were located within or close to the previously reported QTL regions. Five candidate genes related to IgY level by GWAS, and validated here by experimental challenge, are novel and may play critical roles in the modulation of immune responses. It is highly likely that these candidate genes play an important role in the immune response. Findings . Serum concentrations of total IgY of SE-infected and uninfected birds were determined by chicken ELISA kits. Results were expressed as the mean fold-change in gene expression at each time point, from triplicate analyses, using the non-infected control sample at 0.5 dpi as the calibrator (assigned an expression level of 1). Each point-in-time represents the mean ± S.D. of 6 birds. There were barely detectable or no changes in relative expression of the five candidate genes and total IgY concentration in un-challenged controls (data not shown). herein lay an important foundation for further identifying the immunity-related genes and provide useful information for the molecular mechanisms underlying disease-resistance traits in chickens.

Ethics Statement
All of the animal experiments were conducted in accordance with the Guidelines for Experimental Animals established by the Ministry of Science and Technology (Beijing, China). Animal experiments were approved by the Science Research Department (in charge of animal welfare issue) of the Institute of Animal Sciences, CAAS (Beijing, China).

Animal resource
Beijing-You chickens (a Chinese indigenous breed) from the national conservation farm, Institute of Animal Sciences, CAAS, were used. A total of 728 male offspring from 50 sire families were produced in one hatch. The birds were reared in stair-step cages under standard conditions of humidity, temperature, and ventilation, and had access to water and feed ad libitum. All chickens received the conventional inoculation regimen (S7 Table), including vaccination against H5N1 avian influenza virus (AIV) on d 17.

Phenotypic measurements
All blood samples were collected from the brachial vein by venipuncture. At d 49, chickens were injected intramuscularly in the leg with 1 ml of a 0.25% suspension of SRBC. Sera were isolated from blood collected on d 55 (6d post-injection) for immediate determination of anti-SRBC titers by hemagglutination [62] using 96-well plates and the same preparation of SRBC as used for the immunization. Antibody titers were expressed as the log2 of the reciprocal of the highest dilution giving complete agglutination. On d 62, blood smears were prepared from each chicken, stained with May-Grunwald and Giemsa stains and differentially counted [63]. One hundred leukocytes, including granular (heterophils, eosinophils, and basophils) and nongranular (lymphocytes and monocytes), were counted on one slide from each bird, and the heterophil to lymphocyte ratio (H/L) was calculated. Concentrations of IgY and anti-AIV titers were measured in serum collected on d 80. A double antibody sandwich enzyme-linked immunosorbent assay (ELISA) was used for IgY (R&D Systerm, Minneapolis, USA) and a hemagglutination-inhibition diagnostic test was used for titering anti-AIV. The descriptive statistics of phenotypic values of all traits were analyzed using the MEANS procedure in SAS 9.2 software (SAS Institute Inc., Cary, NC).

Genotyping and quality control
Genomic DNA (gDNA) was extracted from blood samples using phenol-chloroform. The DNA samples were genotyped using the Illumina 60K chicken SNP Beadchip (DNA Land-Marks Inc., Saint-Jean-sur-Richelieu, Quebec, Canada). The intensity files generated were processed with the GenomeStudio GT module with default settings [64].
Quality control was performed by PLINK 1.07 [65] and eight chickens were excluded because of high genotyping missing rate (> 95%). SNPs failing to meet any one of the following filtering criteria were excluded: call rate > 90%; minimum minor allele frequency (MAF) > 0.03; P-value for Hardy-Weinberg equilibrium test > 1×10 -6 .

Statistical analysis
To evaluate if population substructure existed among the individuals, classical MDS, as implemented by PLINK 1.07 [65], was performed. The procedures were as follows: first, to ensure uncorrected LD did not distort the analysis, all autosomal SNPs passing the quality control were subjected to LD-based pruning, thus the remaining SNPs within a given window size of 25 had pairwise r 2 < 0.2; second, the pairwise IBS distance among the 720 individuals was calculated using the remaining 7,175 SNPs; finally, the first two MDS dimensions were extracted via the "MDS-plot" command and visualized in R v3.0.1 (www.r-project.org).
Association analyses were conducted to test the additive effects of allele dosage by PLINK, using the linear regression model. The fitted equation was: Where: Y ij is the phenotypic value of a given immune trait for individual i in the j th sire family, b 0 is the intercept, b k-1 , b 50 is the regression coefficient, f is the dummy variable for the sire family, k is an integer ranging from 2 to 50, X ij denote the genotype at the given SNP, and ε ij is independent, normally distributed random residual, with mean 0 and constant variance. For the family variable with 50 categories, 49 binary dummy variables f (coded 0/1), created by the dummy-coding option within PLINK, were fitted in the model. Given the correlation between SNPs in linkage disequilibrium, the strict Bonferroni adjustment appears to be overly conservative [66]. The sum of independent LD blocks plus singleton markers was used to define the number of independent statistical comparisons [67], as described in [68]. Finally, 17,794 independent tests were used to determine the P-value thresholds, including 5% genome-wide significance (2.81E-6, 0.05/17,794) and suggestive association (5.62E-5, 1/17,794). A global view of the p-value results for the six traits from the association analysis was created by the "gap" package [66] in R v3.0.1 (www.r-project.org).

Relative expression of candidate genes identified by the GWAS
Birds and Treatment. Twelve-day-old Beijing-You chickens were from the same conservation stock and were reared as described earlier, in separate cages. Birds were confirmed to be free of Salmonella by culturing fecal samples in buffered peptone water (BPW) overnight with shaking at 150 rpm and spreading on brilliant green agar containing 100 mg/ml nalidixic acid (37°C, 18-24 h). Chickens were randomly divided into 2 groups of 30 each (Salmonella enteritidis challenge and un-challenged controls). Birds in the challenged group were injected intramuscularly into the leg with 0.5 ml of SE bacterial suspension containing 10 8 CFU and the control group was given 0.5 ml saline. Six chickens from each group were randomly chosen, and sacrificed at 12h, 24h, 3 or 6 days post-infection. Blood was collected, allowed to clot, and serum was stored. The spleen was rapidly removed, snap-frozen in liquid nitrogen, and stored at-80°C for mRNA extraction.
Serological test. Concentrations of total IgY in all sampled birds were determined in triplicate by ELISA, as described earlier.
mRNA abundance of candidate genes by RT-qPCR. Total RNA was isolated from each spleen sample with the RNAsimple Total RNA kit (TIANGEN BIOTECH, Beijing, China) and first-strand cDNA was synthesised from 2 μg total RNA using the Promega Reverse Transcription Kit (Promega, Beijing, China), according to the manufacturer's instructions. Power SYBR Green PCR Master Mix (Applied Biosystems, USA) was used to quantify mRNA abundance of the genes by Reverse Transcription-Quantitative PCR (RT-qPCR) with an ABI 7500 Real-time Detection System (Applied Biosystems, USA). The primers were designed with Primer Premier 5.0 from the GenBank sequences (Table 5). The PCR amplification was performed, as previously described [24]. To determine fold-changes in gene expression, the comparative CT method was used [69], calculated as 2 -ΔΔCT . Results were expressed as the mean fold-change in gene expression at each time, from triplicate analyses, using the non-infected control sample at 0.5dpi as the calibrator (assigned an expression level of 1). All analyses were performed using Student's t-test and data are presented as means ± S.D.