Genetic Analyses of Tanzanian Local Chicken Ecotypes Challenged with Newcastle Disease Virus

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.


Introduction
Newcastle disease (ND) is a major threat to poultry globally, and severe outbreaks of the velogenic strains of Newcastle disease virus (NDV) in African local chicken populations often have devastating economic impacts on households. In Africa, local chicken production is characterized by low input production systems [1,2] and serves as a major source of high quality protein (eggs and meat) and financial assets, particularly for children and women [2][3][4]. Because local chickens are managed in backyards or under a free-range scavenging system [2], households often lose entire chicken flocks to poultry diseases such as NDV [5,6]. Farmers usually control NDV through vaccination, but this is an inadequate control method in many small-holder farms in rural Africa because of limited husbandry 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 log 10 . 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.

Genotyping and Quality Control
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].

Population Stratification
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.

Genetic Parameters and Correlations
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.

Genome-Wide Association 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].

Multiple Test Correction
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.

Bioinformatics Analyses
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/).

Population Stratification and Phenotypic Data
Population stratification using Admixture and Multi-Dimensional Scaling showed overlaps among the three chicken ecotypes (Figures 1 and 2). The clustering analyses indicated that ecotypes assigned to birds showed shared ancestry of genotypes among birds across ecotypes ( Figure 1). Two main clusters were identified, with Ching and MoroMid belonging to one discrete cluster and Kuchi to the other cluster. Admixture analyses based on identity by state also showed clear overlaps between the three ecotypes. Birds mainly belonged to two populations, with Ching and MoroMid having higher average proportions of population one (0.78 and 0.75 for Ching and MoroMid, respectively), while Kuchi birds had a higher (0.67) average proportion of population two.
Genes 2019, 10 5 assigned to birds showed shared ancestry of genotypes among birds across ecotypes ( Figure 1). Two main clusters were identified, with Ching and MoroMid belonging to one discrete cluster and Kuchi to the other cluster. Admixture analyses based on identity by state also showed clear overlaps between the three ecotypes. Birds mainly belonged to two populations, with Ching and MoroMid having higher average proportions of population one (0.78 and 0.75 for Ching and MoroMid, respectively), while Kuchi birds had a higher (0.67) average proportion of population two.  In total, 1399 chicks were challenged with NDV. The mean ± standard error for growth rate was 5.12 ± 1.31 and 6.85 ± 2.82 g/d for pre-and post-infection growth rates, respectively ( Table 2). The mean log10 anti-NDV antibody level was 3.45 ± 0.45. The mean log10 viral copy number was 4.72 ± 1.03 and 4.25 ± 1.18 for 2 and 6 dpi, respectively. Mean viral clearance from 26 dpi was 6%. Genes 2019, 10 5 assigned to birds showed shared ancestry of genotypes among birds across ecotypes ( Figure 1). Two main clusters were identified, with Ching and MoroMid belonging to one discrete cluster and Kuchi to the other cluster. Admixture analyses based on identity by state also showed clear overlaps between the three ecotypes. Birds mainly belonged to two populations, with Ching and MoroMid having higher average proportions of population one (0.78 and 0.75 for Ching and MoroMid, respectively), while Kuchi birds had a higher (0.67) average proportion of population two.  In total, 1399 chicks were challenged with NDV. The mean ± standard error for growth rate was 5.12 ± 1.31 and 6.85 ± 2.82 g/d for pre-and post-infection growth rates, respectively ( Table 2). The mean log10 anti-NDV antibody level was 3.45 ± 0.45. The mean log10 viral copy number was 4.72 ± 1.03 and 4.25 ± 1.18 for 2 and 6 dpi, respectively. Mean viral clearance from 26 dpi was 6%. In total, 1399 chicks were challenged with NDV. The mean ± standard error for growth rate was 5.12 ± 1.31 and 6.85 ± 2.82 g/d for pre-and post-infection growth rates, respectively ( Table 2). The mean log 10 anti-NDV antibody level was 3.45 ± 0.45. The mean log 10 viral copy number was 4.72 ± 1.03 and 4.25 ± 1.18 for 2 and 6 dpi, respectively. Mean viral clearance from 26 dpi was 6%.

Genetic Parameter Estimates
Estimates of heritability from single-trait analyses are shown in Table 2. Heritability estimates from single-trait and bivariate analyses were similar. Estimates of heritability for pre-and post-infection growth rates were moderate, 0.35 and 0.21, respectively. Viral load was moderately heritable, with estimates of 0.18 and 0.35 at 2 and 6 dpi, respectively, but the estimate of heritability for viral clearance was not different from zero. The heritability estimate for anti-NDV antibody level (10 dpi) was 0.22.
Estimates of phenotypic and genetic correlations among the traits are shown in Table 3. For phenotypic correlations, pre-and post-infection growth rates were highly correlated (0.54 ± 0.02) and both were positively correlated with antibody levels and negatively correlated with viral load at 2 and 6 dpi. Anti-NDV antibody levels were positively correlated with viral load at 2 and 6 dpi. Viral clearance was positively and negatively correlated with viral load at 2 and 6 dpi, with estimates of 0.18 ± 0.03 and -0.29 ± 0.03, respectively.
Pre-and post-infection growth rates were genetically highly correlated (0.74 ± 0.08) ( Table 3). Viral clearance was genetically positively correlated with pre-and post-infection growth rates and with antibody level; genetically, birds with higher antibody levels at 2 dpi had higher viral clearance. Viral clearance was negatively correlated with viral load at 6 dpi but not significantly (−0.11 ± 0.21). Viral load at 6 dpi was negatively correlated with pre and post-infection growth rate with estimates of −0.23 and −0.13, respectively. Table 3. Estimates (± SE) of genetic (above diagonal) and phenotypic (below diagonal) correlations based on bivariate analyses.

Genome-Wide Association Studies
A total of 1399 animals and 396,055 SNPs remained after quality control and were utilized for GWAS analyses. Because the chicken ecotypes were highly heterogeneous and admixed, analyses were performed Genes 2019, 10, 546 7 of 16 across ecotypes by fitting population proportion (C) to account for the structure of the populations. Using principle component analysis, 71,374 components accounted for 95% of the variance of SNP genotypes across the genome. Bonferroni corrected thresholds of 10 (p-value = 2.05 × 10 −6 ) and 20% (p-value = 7.035 × 10 −6 ) were used to declare suggestive associations in the single SNP GenABEL analyses. Markers significantly associated with pre-and post-infection growth rates, antibody, and viral load at 2 and 6 dpi from single-SNP GenABEL analyses are presented in Table 4. Gensel results (≥0.5% genetic variance in a 1-Mb window) obtained using the Bayes-B method are presented in Table 5.   There were 10, 1-Mb windows on eight chromosomes that together explained 6.5% of genetic variance for pre-infection growth rate. The location of two suggestive significant SNPs in two QTL on chromosomes 3 and 22 ( Figure 3) corresponded with two 1-Mb windows identified for pre-infection growth rate. For post-infection growth rate, 2 windows on two chromosomes explained 1.7% of genetic variance. The window that explained the most genetic variance (1.2%) corresponded with the location for a suggestive significant SNP on chromosome 19 ( Figure 4). One SNP on chromosome 2 ( Figure 5) and one window explaining > 1% genetic variance on chromosome 9 were associated with anti-NDV antibody levels. None of the 1-Mb windows corresponded in location with the significant SNP on chromosome 2 from single-SNP analyses. Four windows on three chromosomes explained a total of 3.8% genetic variance for viral load at 2 dpi. One significant SNP on chromosome 5 ( Figure 6) corresponded in location with a 1-Mb window explaining 2% of genetic variance. Three windows on three chromosomes together explained 13.7% of genetic variance for viral load at 6 dpi. Three SNPs associated with viral load at 6 dpi ( Figure 7B) corresponded in location with a 1-Mb window that explained 12.4% of genetic variance in the Bayes-B analysis (Table 4). Although heritability was low for viral clearance (Table 2), five windows explaining >1% genetic variance were identified. Four significant SNP locations were in correspondence with these windows (Figure 8).

Population Stratification
The admixture results for the three chicken ecotypes (Kuchi, Ching and MoroMid) indicated a common genetic background of the Tanzanian local chicken ecotypes that were evaluated. This admixture can be attributed to inter-mating (breeding) among chickens across different areas of the country and movement of chickens between country districts through trading and market chains. The MoroMid and Ching ecotypes were sampled from districts (regions) that neighbor each other, compared to Kuchi sample districts. The complimentary Admixture ( Figure 1) and MDS (Figure 2) plots support this geography, where Ching and MoroMid birds had similar admixture patterns and clustered together in the MDS plot compared to Kuchi.

Genetic Parameters
Estimates of heritability were moderate for most traits, ranging from 0.18 for viral load at 2 dpi to 0.35 for pre-infection growth rate. Viral clearance had the lowest heritability estimate, at 0.04. These results agree with previous findings from commercial egg laying chickens that were challenged with the LaSota NDV strain [18,19]. Heritability estimates for antibody level at 10 dpi were similar to those for antibody level in two Tanzanian (Kuchi and MoroMid chicks) ecotypes measured at 2 weeks post-infection [20]. To the best of our knowledge, this study is the first to report heritability estimates for growth rates and viral load in local Tanzanian chicken ecotypes. The moderate estimates of heritability suggest that all NDV response traits measured could be improved through selective breeding of chickens to enhance NDV resilience.
Growth rates were negatively and positively genetically correlated with viral load at 2 and 6 dpi, and viral clearance, respectively. Although these correlations were low, this indicates that selection for decreased viral load at 6 dpi and increased viral clearance could increase pre-and post-infection growth rates. In addition, the positive genetic correlations between anti-NDV antibody levels and growth rates indicate that selection for increased antibody levels would also increase growth rate. Our results contrast with the findings of Rowland et al. [18] in commercial layers, which showed unfavorable genetic correlations of viral load at 6 dpi and antibody levels with pre-and post-infection growth rates. The estimates of genetic correlations of anti-NDV antibody levels with viral load at 2 and 6 dpi were not different from zero. This suggests that 10 dpi anti-NDV antibody may not be a good indicator trait for the ability of the virus to replicate within the host. Although we identified moderate genetic correlations between some of the traits, these estimates should be considered with care because of the high standard errors.

Genome-Wide Association Analyses
Two analysis platforms, Gensel and GenABEL, were used for GWAS analyses. We utilized the Bayes B approach in the Gensel software to compute genetic variance accounted for by a 1-Mb window for the NDV response traits. The Gensel software does not allow for fitting random effects in the model, except SNP effects. An R package, GenABEL, was used to identify individual SNPs associated with the NDV response traits. GenABEL allows fitting single SNPs one-by-one during the association analyses.
One suggestive QTL region was identified for post-infection growth rate based on a significantly associated SNP on chromosome 19 at 1.6 Mb. The location of the QTL identified by GenABEL analyses corresponded with a 1-Mb window ( Table 4) that explained 1.2% of genetic variance based on the Gensel analyses. This SNP is within the intron of the AUTS2 activator of transcription and developmental regulator gene. The role of AUST2 is not well known in chickens but previous studies in other species found AUST2 to be associated with various neurological diseases [34,35]. A study conducted in zebrafish reported an increase in cell proliferation in the brain when AUST2 was knocked down [34]. Therefore, this gene could be vital for growth under NDV challenge conditions. Another gene, SBDS ribosome maturation factor, was 793,675 bp downstream of the significant SNP. In mammalian cells, SBDS has been implicated in several pathways, including cell motility [36], regulation of reactive oxygen species, and ribosome biogenesis [37]. The SBDS gene was differentially expressed at 14 dpi in White Leghorn chickens inoculated with Salmonella enterica serovar Enteritidis. The SBDS gene could be important in the regulation of cellular processes during disease challenge in chickens.
Single-SNP analyses revealed one significant suggestive SNP associated with antibody level. A window explaining >1% of genetic variance for antibody level contained 16 genes (Table 6). However, the location of this window did not correspond with the significant SNP from single SNP analyses. One of these genes from this window, HES1, is a critical mediator of canonical notch signaling in lymphocyte development and transformation in mice [38]. HES1 is expressed in B, T, and STEM cell blood lineages, and is important in transmitting Notch functions [39,40].  One significant suggestive QTL for viral load at 2 dpi was identified on chromosome 5 by the GenABEL analyses. This SNP corresponds in location with a 1-Mb window that explained 2% of genetic variance based on the Gensel analyses. The significant SNP was within the intron of the PLEKHH1 gene, which interacts with the MYC transcription factor to activate the transcription of growth-related genes [41]. These results indicate that this SNP could be the first to associate the PLEKHH1 gene with viral load at 2 dpi. Other genes located in this 1-Mb window (Table 6) are reported to be important in immune response, including ZFP36L1, SLC39A9, and ACTN1. ZFP36L1 is an important regulator of innate immune response and may modulate Mkp-1 mRNA basal levels to control p38 MAPK activity during lipopolysaccharide stimulation [42]. SLC39A9 induces an increase in extracellular zinc levels that leads to Akt and Erk activation/phosphorylation in response to B cell-receptor activation [43]. SLC39A9 may be an important gene in the early stages of NDV infection and may play an important role in the chickens' adaptive immune response. Another gene, ACTN1 has been connected to phagocytosis and the immune system [44].
The QTL on chromosome 24 associated with viral load at 6 dpi had the most significant and largest number of SNPs identified in our study. The approximate 1-Mb window that contained these SNPs also had the highest genetic variance explained and contained 36 genes ( Table 6). Some of these genes, including TIRAP, ETS1, KIRREL3, and ST3GAL4, are potential candidates for this QTL. TIRAP had three significant SNPs downstream of it that were associated with viral load at 6 dpi (Table 5). TIRAP is a Toll-interleukin 1 receptor [45] that recognizes pathogens within a host and is part of the microbial pathogen recognition Toll-like receptor system [46]. TIRAP is an adaptor protein that activates TLR4 signaling and is involved in modulating early innate immune response through detection of viral envelope glycoprotein [47]. Another immune response gene in the identified QTL region was ETS1, which is a transcription factor that regulates cytokines and chemokine gene expression [48] and is required for development of natural killer cells [49]. ETS1 may be an important gene in the early stages of NDV infection and could activate the chickens' adaptive immune response. A previous study [19] conducted in commercial layer chickens identified significant SNPs in the same region as the current study for viral load at 6 dpi. That study investigated resistance to NDV of Hy-Line Brown chickens under the effects of heat stress, and authors utilized the GenABEL analysis R package for their association analyses. The two analysis platforms (GenABEL and Gensel) utilized in the current study revealed the same potential candidate genes on chromosome 24 associated with viral load at 6 dpi as identified by Saelao et al. [19]. A parallel study to Saelao et al. [19] investigated resistance to NDV of the same bird population, but without heat stress [18]. That study did not identify any significant SNPs associated to viral load at 6 dpi. A possible reason why our results for chromosome 24 agree with Saelao et al. [19] but not with Rowland et al. [18] could be because Tanzanian local chicken ecotypes in our study were exposed to natural heat stress. Experiments were conducted during the hottest months of January to May, and NDV challenge during this period may have influenced our results.

Conclusions
Our results revealed that heritability estimates for most NDV response traits were moderately high. Growth rates were positively genetically correlated with anti-NDV antibody level and were negatively correlated with viral load at 6 dpi. The genetic correlation between anti-NDV antibody levels and viral load at 2 and 6 dpi were not different from zero. Six suggestive QTL for NDV response traits were identified, some of which were consistent with QTL identified in previous independent studies. The strongest QTL region was identified on chromosome 24 and contained several candidate genes, including EST1, TIRAP, and KIRREL3, that could be vital in the chickens' response to NDV infection. The moderate estimates of heritability and the identified QTL for NDV response traits, suggest that NDV response traits can be improved through selective breeding of chickens to enhance NDV resistance and vaccine efficacy of chickens in regions where NDV is endemic. The variants and important genomic regions identified in this study warrant further investigations to comprehensively understand underlying molecular mechanisms of NDV challenge. Funding: This study was funded by the US Agency for International Development Feed the Future Innovation Lab for Genomics to Improve Poultry AID-OAA-A-13-00080 (H. Zhou). The data supporting the findings of this study will be made available.

Conflicts of Interest:
The authors have no conflict of interest.