Newcastle Disease (ND) is a continuing global threat to domestic poultry, especially in developing countries, where severe outbreaks of velogenic ND virus (NDV) often cause major economic losses to households. Local chickens are of great importance to rural family livelihoods through provision of high-quality protein. To investigate the genetic basis of host response to NDV, three popular Tanzanian chicken ecotypes (regional populations) were challenged with a lentogenic (vaccine) strain of NDV at 28 days of age. Various host response phenotypes, including anti-NDV antibody levels (pre-infection and 10 days post-infection, dpi), and viral load (2 and 6 dpi) were measured, in addition to growth rate. We estimated genetic parameters and conducted genome-wide association study analyses by genotyping 1399 chickens using the Affymetrix 600K chicken SNP chip. Estimates of heritability of the evaluated traits were moderate (0.18–0.35). Five quantitative trait loci (QTL) associated with growth and/or response to NDV were identified by single-SNP analyses, with some regions explaining ≥1% of genetic variance based on the Bayes-B method. Immune related genes, such as ETS1, TIRAP, and KIRREL3, were located in regions associated with viral load at 6 dpi. The moderate estimates of heritability and identified QTL indicate that NDV response traits may be improved through selective breeding of chickens to enhance increased NDV resistance and vaccine efficacy in Tanzanian local ecotypes.
All animal procedures were approved by the Institutional Animal Care and Use Committee of University of California, Davis (#:17853) Breeder chickens consisted of 65 sires and 324 dams from three popular Tanzanian chicken ecotypes, Ching’wekwe (Ching), Kuchi, and Morogoro Medium (MoroMid) that were randomly collected from the five regions of Tanga, Shinyanga, Morogoro, Singida, and Mwaza, representing the central, coastal, and lake zones of Tanzania [22]. Sires and dams were vaccinated using the recommended vaccination schedule by project veterinary personnel. Challenge experiments were conducted for a total of 1399 chicks (Ching (477), Kuchi (315), and MoroMid (607)) from hatch to 38 days of age (doa) across five replicates. All birds were raised under similar conditions with ad libitum access to feed and water. Blood samples were collected from all birds at 27 doa and IDEXX NDV ELISA was used to quantify maternal antibody levels (IDEXX Laboratories, Inc., Westbrook, ME, USA). At 28 doa, all birds were challenged via the oculo-nasal route with 107, 50% embryonic infectious dose (107EID50) of a live attenuated type B1 LaSota lentogenic NDV strain. This lentogenic vaccine strain induces an immune response without causing severe disease and death, allowing for observation of phenotypic responses to the virus. To quantify viral load, viral RNA was isolated from lachrymal fluid samples at 2 and 6 days post infection (dpi) and quantified using qPCR, as described by Rowland et al. [18]. The mean viral RNA was computed per sample and transformed to log10 for downstream analyses. Viral clearance was computed as the difference between viral loads at the two time points divided by viral load at 2 dpi. Blood samples were collected and ELISA was used to measure the anti-NDV antibody levels at 10 dpi, which is the time required for birds to generate an acquired immune response [18,23]. The anti-NDV antibody level data were also transformed to log10. Body weights were measured at hatch and at 7, 14, 21, 28 (0 dpi), 34 (6 dpi), and 38 (10 dpi) doa. Pre- and post-infection growth rates were calculated as grams per day from these weights using linear regression of weight on age. Outliers greater than three standard deviations from the mean were removed for all response traits. Blood samples were collected from chicks before challenge using Whatman FTA cards (Sigma-Aldrich, St. Louis, MO, United States). Genomic DNA was isolated from the FTA cards and genotyping was conducted using the Affymetrix Axiom® 600k Array at GeneSeek (Lincoln, NE, USA). Genotyping Array annotation files were based on the chicken Gallus gallus genome version 5 (Thermo Fisher Scientific Inc., Calsbad, CA, USA). Genotype data quality filtering was performed with Axiom™ Analysis Suite 3.1 (Applied Biosystems, Thermo Fisher Scientific Inc., Calsbad, CA, USA) and included single nucleotide polymorphism (SNP) call rate ≥99% and minor allele frequency ≥0.05. Other Affymetrix genotype metrics used for filtering, with their corresponding cutoffs, are listed in Table 1. After filtering, a total of 396,055 SNPs remained. Imputation of missing genotypes was performed using Fimpute [24]. Genotype quality metrics provided by Affymetrix and the requirements used in quality control filtering. Population structure was examined by constructing a Multi-dimensional Scaling (MDS) plot in two dimensions using the cluster algorithm in the PLINK v1.9 software [25]. Shared ancestry of birds was explored using the Admixture software [26], with the number of subpopulations ranging from 1 to 4. The optimal number of subpopulations was determined based on the lowest cross-validation error rate and was determined to be 2. The generated population proportions for each individual were used in downstream GWAS analyses. Estimation of variance components and heritabilities was done using ASReml 4 [27] based on the following univariate mixed linear animal model: where Y is the dependent phenotype variable: pre- and post-infection growth rate, antibody at 10 dpi, viral load at 2 dpi, and viral load at 6 dpi. Fixed effects included death prior to 10 dpi (D = 0/1), trial (replicate, R = 1–5), and sex (S = male/female). Only one covariate, population proportion (C), obtained as described above, was fitted. Random effects included animal genetic effects (A), using the genomic relationship matrix computed based on procedures reported by [28], X to account for maternal effects, and residuals (e). The dam effect (X) was removed for some traits because it was not significant. For viral load at 2 and 6 dpi, and antibody at 10 dpi, qPCR plate (55 and 60 plates at 2 and 6 dpi, respectively) and replicate plate (46), respectively, were added as fixed effects. Phenotypic variance was obtained by summation of variance due to animal genetic, dam, and residuals. Heritability was computed as a ratio of the estimates of animal genetic to phenotypic variance. Bivariate animal models were used to estimate pairwise genetic and phenotypic correlations between traits, with the same fixed and random effects as specified for the univariate analyses. Two approaches were used for whole genome association analyses. First, a Bayesian approach called Bayes B [29], as implemented in the Gensel software [30], was used to compute the genetic variance accounted for by every-one mega base (Mb) window of SNPs. Model (1) above was used for the analyses but with A was replaced by SNP genotype effects. The genotypes at all SNPs, coded as −10/0/10 for AA/AB/BB, respectively, were fit in a Bayes B approach using a prior probability of the SNP having no effect (pi) of 0.999. In total, 41,000 iterations of a Markov chain were run for each trait, with 1000 iterations burn-in and 100 iterations as the output frequency. We used Bayes B instead of Bayes C approach because it performs better in detecting QTL windows compared to Bayes C. For the second approach, single SNP association analyses were conducted using the R package GenABEL [31], with a hierarchical generalized linear model [32]. The same fixed (class) and covariate effects described in the Gensel analyses were used in the GenABEL analyses. A polygenic model was fitted using the “polygenic_hglm” function, with a genomic kinship matrix that was created using the ibs() function. To test for association between a trait and genotypes at a SNP for related individuals, the “mmscore” function was used and residuals were obtained from the polygenic_hglm analysis. The mmscore function was performed using the formula described by Rowland et al. [18]. To determine significance thresholds for the GenABEL analyses, multiple test correction was performed. Suggestive genome-wide significance thresholds of 10 and 20% were computed using a modified Bonferroni correction as 0.1 or 0.2 divided by the number of independent tests. To determine the number of independent tests, the SNP genotype data were separated by chromosome and then further divided to form chromosomal segments such that the number of SNPs were equal to half the number of animals, as described by Waide et al. [33]. The number of independent tests in each segment was determined by the minimum number of principle components required to account for 95% of the variance among genotypes in each segment. The total number of independent tests was the sum across segments. Gene annotation for 1-Mb windows that explained more than 1% of genetic variance was obtained using NCBI’s Genome Data Viewer on the chicken genome version Gallus gallus 5 (https://www.ncbi.nlm.nih.gov/genome/gdv/)
N/A