Hamzić et al.

Genet Sel Evol (2015) 47:91

DOI 10.1186/s12711-015-0170-0
Ge n e t i c s
Se l e c t i o n
Ev o l u t i o n

RESEARCH ARTICLE Open Access

Genome‑wide association study

and biological pathway analysis of the Eimeria
maxima response in broilers
Edin Hamzić1,2,3, Bart Buitenhuis3, Frédéric Hérault4, Rachel Hawken5, Mitchel S. Abrahamsen5,
Bertrand Servin6, Jean‑Michel Elsen6, Marie‑Hélène Pinard ‑ van der Laan1,2 and Bertrand Bed’Hom1,2*

Abstract
Background: Coccidiosis is the most common and costly disease in the poultry industry and is caused by protozo‑
ans of the Eimeria genus. The current control of coccidiosis, based on the use of anticoccidial drugs and vaccination,
faces serious obstacles such as drug resistance and the high costs for the development of efficient vaccines, respec‑
tively. Therefore, the current control programs must be expanded with complementary approaches such as the use
of genetics to improve the host response to Eimeria infections. Recently, we have performed a large-scale challenge
study on Cobb500 broilers using E. maxima for which we investigated variability among animals in response to the
challenge. As a follow-up to this challenge study, we performed a genome-wide association study (GWAS) to identify
genomic regions underlying variability of the measured traits in the response to Eimeria maxima in broilers. Further‑
more, we conducted a post-GWAS functional analysis to increase our biological understanding of the underlying
response to Eimeria maxima challenge.
Results: In total, we identified 22 single nucleotide polymorphisms (SNPs) with q value <0.1 distributed across
five chromosomes. The highly significant SNPs were associated with body weight gain (three SNPs on GGA5, one
SNP on GGA1 and one SNP on GGA3), plasma coloration measured as optical density at wavelengths in the range
465–510 nm (10 SNPs and all on GGA10) and the percentage of β2-globulin in blood plasma (15 SNPs on GGA1 and
one SNP on GGA2). Biological pathways related to metabolic processes, cell proliferation, and primary innate immune
processes were among the most frequent significantly enriched biological pathways. Furthermore, the network-based
analysis produced two networks of high confidence, with one centered on large tumor suppressor kinase 1 (LATS1)
and 2 (LATS2) and the second involving the myosin heavy chain 6 (MYH6).
Conclusions: We identified several strong candidate genes and genomic regions associated with traits measured
in response to Eimeria maxima in broilers. Furthermore, the post-GWAS functional analysis indicates that biological
pathways and networks involved in tissue proliferation and repair along with the primary innate immune response
may play the most important role during the early stage of Eimeria maxima infection in broilers.

Background and cause coccidiosis: E. brunetti, E. necatrix, E. tenella,

Coccidiosis is an animal parasitic disease caused by E. acervulina, E. maxima, E. mitis, and E. praecox. E.
protozoans belonging to the Coccidia subclass. In the maxima is the most immunogenic of the seven species
chicken, seven species of the genus Eimeria are infectious [1] and mostly infects the lining of the jejunum, causing
mucoid enteritis [2]. Chicken coccidiosis is one of the
most common and costly diseases currently affecting the
poultry industry, with worldwide costs caused by produc-
*Correspondence: [email protected] tion losses as well as by prevention and treatment actions
2
UMR1313 Animal Genetics and Integrative Biology Unit, INRA, Domaine
de Vilvert, 78350 Jouy‑en‑Josas, France that are estimated to exceed USD 3 billion per year [3, 4].
Full list of author information is available at the end of the article Current coccidiosis management is based on the use of

© 2015 Hamzić et al. This article is distributed under the terms of the Creative Commons Attribution 4.0 International License
(http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium,
provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license,
and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/
publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 2 of 17

anticoccidial drugs and vaccination [5]. The first anticoc- was performed to further understand the biology of the
cidial drugs, sulfonamides, began to be used in the early underlying response to the E. maxima challenge. The
1940s, and then over the years, several different classes of functional analysis comprised the following two inde-
drugs were developed and extensively used for the con- pendent approaches: a biological pathway analysis based
trol of coccidiosis in broiler production [6]. However, on the publicly available Kyoto encyclopedia of genes and
the future use of anticoccidial drugs has caused concern genomes (KEGG) pathways and a network-based analysis
due to the tendency of Eimeria species to rapidly develop that was performed using the ingenuity pathway analysis
resistance to drugs [7] as well as public dissatisfaction (IPA) software.
regarding the presence of chemical residues in food. In
addition, the development of efficient multiple-species Methods
live vaccines is primarily limited by their high economic Ethics statements
costs [5]. All procedures were conducted under License No.
Due to these issues, the current control programs A176661 from the Veterinary Services, Charente Mari-
must be expanded using a complementary approach, time, France and in accordance with guidelines for the
including the application of genetics to improve the Care and Use of Animals in Agricultural Research and
host response to Eimeria infection. Genetic approaches Teaching (French Agricultural Agency and Scientific
and manipulations have been shown to have a small Research Agency) (http://www.gouvernement.fr/en/
effect at each generation, but also a cumulative effect culture-education-and-research).
and is a long-term, cost-effective and environmentally
friendly way of keeping livestock animals healthy [8]. Experimental population and phenotyping
The first studies indicating that genetic diversity may In this study, we used phenotype data collected during a
account for differences in susceptibility to coccidiosis large-scale Eimeria maxima challenge study on Cobb500
were published in the 1940s and 1950s [9, 10]. Consid- broilers [19]. The challenge study was performed using
erable variation in coccidiosis susceptibility has been 2024 Cobb500 broilers randomly distributed in 44 (chal-
observed between different chicken breeds [11]. The lenge) and two (control) litter pens (3 m × 1 m) each
first successful performance of divergent selection for containing 44 birds. For the GWAS, we used only data
resistance and susceptibility to acute cecal coccidiosis collected from the challenged animals, with control ani-
was performed by Johnson and Edgar [12]. Likewise, mals excluded from the analysis. Traits were measured
marked differences in the response to Eimeria infection at two levels: “global phenotyping”, which was performed
were observed between inbred and outbred lines [13, on all animals and “detailed phenotyping”, which was per-
14]. The availability of genome-wide dense markers has formed on a subset of 176 animals. The experimental lay-
enabled the identification of several quantitative trait out of the challenge is presented in Fig. 1.
loci (QTL) regions associated with resistance to Eimeria BWG, HEMA, BT and PC were measured on all ani-
tenella and Eimeria maxima in experimental popula- mals. BWG was calculated as BWG = (BW at day 22-BW
tions [15–17]. In addition, several genes that are located at day 15)/BW at day 15. We used relative BWG with
in highly significant previously detected QTL [16] were respect to day 15 to focus solely on the BWG in response
also found to be differentially expressed in a follow-up to the challenge. HEMA levels and PC were measured
transcriptome study [18]. from the blood samples that were collected on days 16
Recently, we conducted a large-scale challenge study and 23.
with Cobb500 broilers using E. maxima as the infective The subset of 176 animals was chosen by selecting two
agent, and high variability was observed in the meas- birds among those with the lowest and two birds among
ured traits among challenged animals [19]. The measured those with the largest BWG from each challenged animal
traits included a wide range of physiological, immuno- pen. Animal selection within each pen was based on the
logical and disease resistance parameters. This study is ranking per pen. The subset of 176 animals with extreme
a follow-up on the aforementioned large-scale challenge BWG values was used for further detailed phenotyping.
study in which we assessed the effects of the E. maxima Detailed phenotyping included LS (duodenum and jeju-
challenge on the measured traits and also evaluated the num), OC, BC and PPP. PPP was assessed using capillary
level of variability of the measured traits. Taking advan- electrophoresis, which separates the protein components
tage of the size and structure of this large-scale challenge into five major fractions by size and electrical charge:
study, we performed a genome-wide association study prealbumins, albumins, α-1 globulins, α-2 globulins, α-3
(GWAS) to identify genomic regions underlying the globulins, β-1 globulins, β-2 globulins and γ globulins.
broiler response to E. maxima. In addition, based on the BC included blood cell count (BCC) and red blood indi-
results of the GWAS, a post-GWAS functional analysis ces (RBI). A detailed description of the methodologies
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 3 of 17

LARGE-SCALE STUDY
Before challenge After challenge

Blood sampling
Challenge (50000 oocysts) Blood sampling
Animal delivery 44 x 44 infected animals Fecal sampling for OC
Identification 2 x 44 controls Slaughter

Day 1 Day 8 Day 15 Day 16 Day 22 (Day 6 p.i.) Day 23 (Day 7 p.i.) After experiment

BW BW BW BW BT, LS BC, BWG, HEMA,

OC, PC, PPP
Global phenotyping and detailed phenotyping
Fig. 1 Experimental layout of the large-scale challenge study. In total, 2024 1-day-old broilers were used in the experiment with 88 control animals
and 1936 challenged animals. The challenge was performed on day 16 of the experiment by inoculating 50,000 Eimeria maxima oocysts. Traits were
measured at two levels: “global phenotyping” and “detailed phenotyping”. Global phenotyping was performed on all animals and included body
weight (BW), plasma coloration (PC), body temperature (BT) and hematocrit (HEMA) levels. Detailed phenotyping was performed on the subset
of 176 animals and included lesion score (LS), oocyst count (OC), plasma protein profiles (PPP), and blood composition (BC). Body weight gain
(BWG) was calculated using the following formula (BW at day 22-BW at day 15)/BW at day 15. PC is optical density of blood plasma measured for 44
wavelengths (every 5 nm, 380–600 nm) on days 15 and 22. BT was measured at day 23, and HEMA was measured from blood samples obtained on
days 16 and 23. BC, PPP, OC and LS were assessed using samples obtained on day 23. Please refer to Hamzic et al. [19] for a detailed description of
the trait measurements

used for measuring all traits was reported by Hamzic (LGE22C19W28_E50C23 and LGE64) and chromosome
et al. [19]. Z. Mean, median and standard deviation for base pair
distances between neighboring SNPs across chromo-
Genotype data somes are in Table S2 (See Additional file 2: Table S2).
DNA was extracted from blood samples obtained from

typing was performed using a 580 K Affymetrix® Axiom®

1972 animal samples at day 16 of the experiment. Geno- Genome‑wide association analysis
The genome-wide association analyses were performed
HD genotyping array (Affymetrix, Santa Clara, USA) [20] using the Genome-wide Efficient Mixed Model Associa-
at a commercial laboratory (GeneSeek, Lincoln, USA). tion (GEMMA) algorithm [25] for the univariate linear
Quality control of the genotype data was performed mixed model (LMM). In this study, we used Cobb500
using PLINK 1.9 [21–23] and included the sample call broilers, which are the final products of a four-way cross-
rate (>98 %), SNP call rate (>98 %), minor allele fre- breeding scheme, indicating the presence of a strong
quency (>2 %) and removing SNPs with extreme F-sta- population structure. The linear mixed model approach
tistic values. A total of 138,568 SNPs out of 580,961 were was used since this method has been proven to suc-
removed during the quality control (See Additional file 1: cessfully account for population structure in associa-
Table S1). SNP thresholds and sample call rates were set tion mapping studies [26]. In the context of the LMM,
at 98 % and all SNPs and individuals with a call rate less the correction for population structure is performed by
than 98 % were excluded from the analysis. We did not creating the genomic relationship matrix (K) that mod-
perform imputation of the residual missing genotypes. els the structure present in the analyzed population by
The F-statistics distribution for the SNPs that remained estimating the contribution of genetic relatedness to the
after the quality control steps was assessed, and all SNPs phenotypic variance. K presents a pairwise relationship
outside the adjusted interquartile range were excluded between individuals, and its structure is also influenced
from the analysis. The interquartile range for skewed dis- by population structure, family structure and cryptic
tribution of F-statistics was calculated as suggested by relatedness.
Hubert and Vandervieren [24]. Sex was assessed using The model used for the analysis is presented in expres-
PLINK 1.9 by analyzing SNPs on the Z and W chromo- sion (1):
somes. Detected females were included in the analysis
with sex considered as a covariate (See Additional file 1:
y = Wα + xβ + u + ǫ, (1)
Table S1). Finally, after quality control, we obtained a where y is an n-vector of observations (or trait meas-
dataset that included 443,587 SNPs distributed along urements) for n individuals, W is an n × c matrix of
28 autosomal chromosomes, two linkage groups covariates which contains information about pen and
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 4 of 17

sex including a column of 1s for the general mean; α is package [30], and the curated genes were assigned to spe-
a c-vector of the corresponding coefficients, including cific pathways. The biological KEGG pathways represent
the intercept, x is an n-vector of marker genotypes, β is a collection of biological pathway maps that integrate
the effect size of the marker, u is an n-vector of polygenic many units, including genes, proteins, RNAs, chemical
effects, and ε is an n-vector of residual effects. compounds, and chemical reactions, as well as disease
For polygenic effects, u follows a multivariate normal genes and drug targets, which are stored as individual
distribution (MVNn) (u ∼ MVNn (0, τ −1 K)), with λ entries in other KEGG databases. The chicken genome
indicating the ratio between the genetic variance (more KEGG PATHWAY contains 162 pathways associated
precisely the variance explained by the SNPs) and vari- with 4342 genes. The list of genes with assigned SNPs
ance of the residuals, τ−1 is the variance of the residuals, from the GWAS and the list of genes with assigned bio-
and K is a known n × n, which is the genomic relation- logical pathways (KEGG) were combined, ultimately
ship matrix, calculated using a n × p matrix of genotypes resulting in a list of 52,204 SNPs assigned to 162 biologi-
(X). Expression (2) was used to calculate K where xi is cal KEGG pathways (See Additional file 3: Table S3).
the ith column representing genotypes of the ith SNP in
the X matrix, xi is the sample mean, and 1n is an n × 1 Statistical analysis for biological pathway analysis
vector of 1s. Residual effects are presented as an ε vector As described in the previous step, a set of SNPs that
of length n following the multivariate normal distribu- were mapped to genes were further assigned to the
tion (MVNn) (ε ∼ MVNn (0, τ −1 In )) where In is an n × n corresponding biological pathways. For each biologi-
identity matrix. cal pathway that was characterized by a set of SNPs, an
appropriate summary of statistics was constructed as
1
p
described in detail by Jensen et al. [31]. The statistics
K= (xi − 1n xi )(xi − 1n xi )T . (2)
p summary was based on the negative log-transformed
t=1
p values from the association of individual SNPs to the
To correct for multiple hypothesis testing, the false-dis- traits. By summing these negative log-transformed p val-
covery rate (FDR) was calculated for each SNP from the ues, we imitated a genetic model that captures variants
distribution of p-values: SNPs with an FDR less than 0.1 with small to moderate effects [32, 33].
were considered significant [27]. The results are shown The observed summary statistics for a particular set of
in Manhattan plots constructed by the qqman R package SNPs were compared with an empirical distribution for
[28]. the summary statistics of random samples of SNP sets of
the same size using a permutation approach. Considering
Biological pathway analysis that the distribution of summary statistics is affected by a
In general, biological pathway analysis is used to test correlation structure of closely linked SNPs, as the con-
the association between a curated set of genes (biologi- sequence of linkage disequilibrium, the following proce-
cal pathways) and a trait of interest. This approach tests dure was used for statistical testing.
for the cumulative effect across many genes, which ena- The vector of observed SNPs with corresponding test
bles the detection of effects at the biological pathway statistics was ordered according to the physical posi-
level. Biological pathway analysis has been conducted as tion on the genome. SNPs were then mapped to genes
described in the following paragraphs. This analysis was and consequently to a biological pathway as described
divided into two main steps: (1) assigning SNPs to their in the first step. The elements in this vector were num-
corresponding annotated genes and assigning these genes bered 1,2,…,N, and the permutation was performed in
to their corresponding biological pathways (KEGG), and two steps. The first step included randomly picking an
(2) statistical testing for the biological pathway (set of element (ej) from this vector. This jth test statistic was
genes with assigned list of SNP) based on the test statis- the first element in the permuted vector, and the remain-
tics obtained in the genome-wide association analysis. ing elements were ordered ej+1, ej+2, …, eN, e1, e2, …, ej-1
according to their original numbering. Therefore, all
Assigning SNPs to genes and biological pathways elements from the original vector were then shifted to
The SNPs used in the GWAS were mapped to the a new position starting with ej; however, the gene posi-
ICGSC Gallus_gallus-4.0 chicken genome assembly tion was kept fixed with respect to the original one. The
(GCA_000002315.2) using the NCBI2R R package [29]. second step involved the computation of summary statis-
A list of annotated genes corresponding to the SNPs tics for each set of SNPs based on the original position
used in the genome-wide analysis was retrieved. For this of the set of SNPs in the original vector of test statistics.
study, the publicly available biological KEGG pathways The connections between SNPs and genes were broken
were downloaded using the Bioconductor KEGGREST while keeping the correlation structure among the test
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 5 of 17

statistics. Steps 1 and 2 were repeated 1000 times, and In addition to these traits, animals from the subset of
from this empirical distribution of summary test statis- 176 were phenotyped for lesion score (LS), oocyst count
tics for each set of SNPs, a p value was obtained. This (OC), blood composition (BC) and plasma protein pro-
empirical p value corresponds to a one-sided test of the files (PPP) (See Additional file 5: Table S5). BC included
proportion of randomly sampled summary statistics that two sets of measurements: blood cell count (BCC) and
were larger than the observed summary statistic with the red blood indices (RBI). BCC included erythrocyte, leu-
arbitrary significance level set to 0.01. kocyte, lymphocyte, heterophil, and thrombocyte counts
as well as the percentage of lymphocytes and heterophils
Network‑based analysis of the total number of leukocytes. RBI included hemo-
The network modeling was performed using the inge- globin content, mean corpuscular volume (MCV), mean
nuity pathway analysis (IPA) tool as a complementary corpuscular hemoglobin (MCH) and mean corpuscular
approach to the KEGG biological pathway analysis. The hemoglobin concentration (MCHC). PPP were assessed
ingenuity pathways database is the manually curated using protein capillary zone electrophoresis, and the pro-
database of previously published relationships on human files included the following fractions: prealbumin, albu-
and mouse biology [34]. The gene input list included min, α1-globulin, α2-globulin, α3-globulin, β1-globulin,
genes that contained SNPs with p values below the β2-globulin, and γ-globulin. Detailed results from the
inferred genome-wide threshold (P < 10−4) for all traits large-scale challenge study were reported by Hamzic
on which GWAS was performed (See Additional file 4: et al. [19]. Descriptive statistics of the traits measured on
Table S4). In this manner, we restricted the list of putative the challenged animals are in Table S5 (See Additional
candidate genes and exploited their documented interac- file 5: Table S5).
tions in biological pathways related to the study. The lists
of SNPs (P < 10−4) were assigned to their respective genes Genome‑wide association study (GWAS)
using the biomaRt R package [35], producing the gene list In total, 22 SNPs were significantly associated (q value
for each trait used in the GWAS. The obtained Human <0.1) with the measured traits and were distributed over
Genome Gene Nomenclature Committee (HGNC) iden- five chicken chromosomes (See Additional file 6: Table
tifiers were mapped onto networks that are available in S6). The most significant SNPs were associated with
the Ingenuity Pathway repository. In the case of PC, we BWG, PC for wavelengths from 465 to 510 nm and the
merged all gene lists obtained for individual wavelengths percentage of β2-globulin in blood plasma of one of the
into one collective gene list. Furthermore, we also created PPP fractions.
a global list of genes by combining all gene lists together. The GWAS identified five SNPs that were significantly
In summary, we performed IPA using 27 gene lists of associated with BWG (Fig. 2), and the quantile–quantile
which 25 gene lists corresponded to each individual trait, (Q–Q) plot for BWG showed a large deviation from the
a gene list for PC obtained by combining gene lists for all distribution under the null hypothesis, indicating that
wavelengths and the global list of genes. strong associations were observed (See Additional file 7:
Figure S1). The five observed SNPs are located on GGA1,
Results 3 and 5. These SNPs explain 7.5 % of the total variance for
Experimental population and phenotyping BWG between days 15 and 23 of the challenge study. The
In this study, we analyzed the data collected during the significantly associated genomic region on GGA5 was
large-scale challenge study on Cobb500 broilers with E. located in the upstream region of the THBS1 gene, which
maxima as the infective agent [19]. The experimental encodes the multi-domain matrix glycoprotein termed
scheme of the large-scale challenge study is illustrated thrombospondin-1. Similarly, the significantly associated
in Fig. 1. The traits were measured at two levels: “global genomic regions on GGA1 and GGA3 are in the vicinity
phenotyping” performed on all 1936 challenged animals of MGAT4C and KCNK3, respectively.
and “detailed phenotyping” performed on a subset of 176 For PC, the AX-75604378 SNP located on GGA10 was
animals (Fig. 1). To form the subset of 176 animals, two significantly associated with PC values measured for
animals among those with the lowest and two animals wavelengths ranging from 465 to 510 nm (Fig. 3). Figure 3
among those with the highest body weight gain (BWG) shows the Manhattan plot and Q–Q plot for PC meas-
were selected from each pen containing challenged ured at 485 nm. The Q–Q plot shows a strong deviation
animals. from the distribution under the null hypothesis, which
The traits measured on all challenged animals included indicates the presence of a strong association between
BWG, plasma coloration (PC) measured as optical den- the SNP and PC values (See Additional file 8: Figure S2).
sity (OD) of blood plasma in the 380–600 nm range, The associated SNP explains between 2.1 and 2.3 % of the
body temperature (BT) and hematocrit level (HEMA). total variance depending on the measured wavelength
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 6 of 17

Fig. 2 Manhattan plot for body weight gain. Manhattan plot of genome-wide −log10 (p values) for body weight gain. P values were adjusted
using the false-discovery rate (FDR) at a significance level of q value <0.1. SNPs labeled in green have q value <0.1

Fig. 3 Manhattan plot for plasma coloration (485 nm). Manhattan plot of genome-wide −log10 (p values) for plasma coloration (485 nm). P values
were adjusted using the false-discovery rate (FDR) at a significance level of q value <0.1. SNPs labeled in green have q value <0.1

(See Additional file 6: Table S6). The AX-75604378 SNP genome assembly. This list of SNPs and their correspond-
is a non-synonymous polymorphism present in the ing genes was combined with the publicly available bio-
MAN2C1 gene. logical KEGG pathway database. Finally, we obtained
Among the traits measured in the subset of 176 ani- a list containing 52,204 SNPs that were included in the
mals, the genomic region located between 52.31 and GWAS; these were assigned to 162 of the biological path-
52.63 Mb on GGA1 and the SNP AX-76165289 that ways associated with 4342 genes in the chicken genome
mapped to GGA2 were significantly associated with the KEGG PATHWAY repository (see “Methods”). The list
percentage of β2-globulin in the blood plasma (Fig. 4). of KEGG pathways which were significantly (P < 0.05)
The genomic region on GGA1 contains 16 SNPs that enriched with genes in genomic regions associated with
were associated with the percentage of β2-globulin in the measured traits is in Table S7 (See Additional file 10:
blood plasma. Individually, the SNPs on GGA1 explain Table S7). The distributions of the most frequent signifi-
approximately 13.4 % of the total variance, and SNP cant biological pathways differed considerably accord-
AX-76165289 explains 14.2 % of the total variance. The ing to whether all measured traits or all measured traits
Q–Q plot for β2-globulin also shows that strong asso- except PC measurements were considered (See Addi-
ciations were observed (See Additional file 9: Figure S3). tional file 11: Figure S4). Several biological pathways were
SNP AX-76165289 is located in the FHOD3 gene on characteristic of the genomic regions associated with PC
GGA2, while five SNPs on GGA1 (between 52.31 and measurements (See Additional file 11: Figure S4), and
52.63 Mb) are located in the LARGE gene (See Additional due to the multiple wavelength measurements, they were
file 6: Table S6). more frequent, which was not the case when the results
were summarized without considering PC measure-
Biological pathway analysis ments (See Additional file 11: Figure S4). For example,
To conduct the biological pathway analysis, we used all the phenylalanine, tyrosine and tryptophan biosynthesis
SNPs that were included in the GWAS. These SNPs were pathway was one of the most frequent pathways when
mapped to their corresponding genes using the latest considering only the PC measurement results and was
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 7 of 17

Fig. 4 Manhattan plot for the percentage of β-globulin. Manhattan plot of genome-wide −log10 (p values) for the percentage of β-globulin. P
values were adjusted using the false-discovery rate (FDR) at a significance level of q value <0.1. SNPs labeled in green have q value <0.1

not common for other measured traits [(See Additional Considering the long list of measured traits for BC
file 11: Figure S4) and Fig. 5]. Therefore, we present the and PPP, we summarize the most interesting results.
PC results (Fig. 5) and all other measured pathways sepa- Regarding the percentage of lymphocytes and hetero-
rately (Fig. 6). phils, we observed 15 and 13 significant biological path-
For PC measured from 380 to 600 nm, we detected 20 ways, respectively, with 13 pathways being common for
significant biological pathways (See Additional file 10: both traits (Fig. 6). Furthermore, the percentages of lym-
Table S7). The most frequent significant biological path- phocytes and heterophils were associated with the larg-
ways for PC across all measured wavelengths were: est number of significant biological pathways (Fig. 6).
phenylalanine, tyrosine and tryptophan biosynthesis, Regarding the number of thrombocytes and erythro-
endocytosis, purine metabolism, glycosphingolipid bio- cytes, 13 and 10 significant biological pathways were
synthesis—lacto and neolacto series, ether lipid metab- observed, respectively (See Additional file 10: Table
olism, caffeine metabolism, spliceosome, and the p53 S7), and these two traits had four significant pathways
signaling pathway (Fig. 5). Each of the mentioned biologi- in common (Fig. 6). Regarding PPP, both the percent-
cal pathways occurred more than 10 times as significant ages of α1-globulin and α2-globulin were associated
when considering all PC measurements (See Additional with 10 significant biological pathways (See Additional
file 10: Table S7). file 10: Table S7). Percentages of prealbumin and albumin
In total, seven biological pathways were significantly shared the following significant pathways: amino sugar
associated with BWG (See Additional file 10: Table S7). and nucleotide sugar metabolism, glycosaminoglycan
Glycosaminoglycan biosynthesis—heparan sulfate/ biosynthesis—chondroitin sulfate/dermatan sulfate and
heparin, glycerolipid metabolism, primary bile acid bio- glycosylphosphatidylinositol (GPI)-anchor biosynthesis
synthesis, melanogenesis, proteasome and regulation (Fig. 6).
of autophagy were the most frequent significantly asso- The three most frequent significant pathways, consid-
ciated pathways across all measured traits except PC ering all measured traits, including PC, were phenylala-
(Fig. 6). In total, 13 and 12 biological pathways were sig- nine, tyrosine and tryptophan biosynthesis, endocytosis
nificantly associated with BT and HEMA, respectively and purine metabolism (See Additional file 11: Figure
(See Additional file 10: Table S7). For BT, the strongest S4). The most frequent significant pathways, considering
associations were observed with metabolic pathways all traits except PC, included purine metabolism, glycosa-
(P = 0.002) and the ErbB signaling pathway (P = 0.003). minoglycan biosynthesis—heparan sulfate/heparin, the
For HEMA, the strongest associations were observed pentose phosphate pathway, and peroxisome and vascu-
with vascular smooth muscle contraction (P < 0.001) lar smooth muscle contraction (Fig. 6) and (See Addi-
and pantothenate and CoA biosynthesis (P = 0.004) (See tional file 11: Figure S4).
Additional file 10: Table S7). For duodenal and jejunal LS,
five and four biological pathways were significant, respec- Network‑based analysis
tively (See Additional file 10: Table S7). Among these, The network-based analysis was performed using the
vascular smooth muscle contraction and porphyrin and ingenuity pathway analysis (IPA) tool. Gene lists used
chlorophyll metabolism were significantly enriched for as input files for the IPA tool were obtained by extract-
both duodenal and jejunal LS (See Additional file 10: ing genes containing SNPs, with an inferred genome-
Table S7). wide significance level (P < 10−4) for all QTL. The IPA
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 8 of 17

Distribution of the biological KEGG pathways

for plasma coloration
Wnt signaling pathway

Steroid biosynthesis

Spliceosome

RNA transport

Purine metabolism

Phenylalanine, tyrosine and tryptophan biosynthesis

Peroxisome

p53 signaling pathway

Non−homologous end−joining

MAPK signaling pathway

Lysine degradation

Hedgehog signaling pathway

Glyoxylate and dicarboxylate metabolism

Glycosphingolipid biosynthesis − lacto and neolacto series

Glycosphingolipid biosynthesis − globo series

FoxO signaling pathway

Ether lipid metabolism

Endocytosis

Caffeine metabolism

Adherens junction

380 395 410 425 440 455 470 485 500 515 530 545 560 575 590

Plasma coloration, nm
Fig. 5 Distribution of the biological pathways that were significantly enriched with genes in genomic regions associated with plasma coloration
(PC). Distribution of the biological pathways significantly (P < 0.05) that were enriched with genes in genomic regions associated with plasma
measurement for all 45 wavelengths. Coloration of blood plasma was measured as the level of absorbance for 45 wavelengths in the range from
380 to 600 nm. This figure illustrates significant KEGG pathways for PC measurements in the range from 380 to 600 nm

tool produced 76 putative networks using 27 gene lists. IPA score 41 (Fig. 7) and the second network with myo-
The majority of networks were related to general molec- sin heavy chain 6 (MYH6) in the center with IPA score 51
ular and cellular processes such as cell to cell signaling, (Fig. 8).
nucleic acid metabolism and replication, cell cycle, with
several of the networks involved in more specific path- Discussion
ways such as gastrointestinal diseases and organismal Identifying genomic regions that underlie the response to
injury and abnormalities. In the context of the large-scale Eimeria maxima in broilers allows us to understand the
challenge study, the most informative networks were the associated molecular mechanisms and provides us with
network of interacting molecules grouped around large candidate genes and genomic regions that can be used
tumor suppressor kinases 1 (LATS1) and 2 (LATS2) with in breeding for improved resistance to coccidiosis. The
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 9 of 17

Fig. 6 Distribution of the biological pathways that were significantly enriched with genes in genomic regions associated to all measured traits
except plasma coloration (PC). Distribution of significantly (P < 0.05) that were enriched with genes in genomic regions associated with all meas‑
ured traits except plasma coloration (PC). This figure summarizes significant biological pathways that occurred more than three times for measured
traits

host response to Eimeria is a complex trait controlled by for which selection would be performed within the pure
a wide range of biological processes, which are in turn lines. This approach may potentially improve general
controlled by many genes with a small effect and a small innate immunity as well as resistance to specific patho-
number of genes with a moderate or large effect. Several gens such as Eimeria species.
QTL regions that are significantly associated with traits However, a more detailed understanding of the genetic
measured after the response to Eimeria challenge were mechanisms that control the response to Eimeria infec-
detected in F2 crosses obtained from experimental lines tion, as in the case of all other complex traits, requires
with different degrees of susceptibility to coccidiosis [16, a sufficiently large sample size, dense SNP coverage that
17, 36]. In addition, similar approaches have been used can exploit the linkage disequilibrium and informative
for the identification of genomic regions associated with phenotypes [39]. Therefore, we performed the large-

mals, which were genotyped using the 580K Affymetrix®

innate and adaptive immunity in laying hens [37] as well scale challenge study on 1936 commercial Cobb500 ani-

Axiom® high-density genotyping array, providing con-

as survival rate and aimed at improving general robust-
ness, especially in laying hens [38]. In addition, infor-
mation regarding the genomic regions that are strongly siderable statistical power for GWAS for traits measured
associated with desirable traits can be incorporated in on all 1936 animals. Taking advantage of the size and
commercial poultry breeding programs. Regarding broil- structure of the challenge experiment [19], we conducted
ers, haplotypes that are associated with desirable traits a GWAS to obtain more information on the underlying
identified in the final product of the four-crossway breed- genetic determinism of the response to E. maxima in
ing scheme can be traced back in pedigreed populations broilers. In addition, a post-GWAS functional analysis
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 10 of 17

Fig. 7 Network of interactions between GWAS candidate genes using ingenuity pathway analysis (IPA). The network shows molecular interac‑
tions between the products of the candidate genes enriched for significantly associated SNPs (P < 10−4) with large tumor suppressor kinases 1
(LATS1) and 2 (LATS2) as the most interesting candidates. Relationships were determined using information contained in the IPA repository. The blue
label indicates the genes that were enriched for significantly associated SNPs. The network was obtained using the global gene list determined by
combining gene lists for all traits

was performed to further understand the biology of the The KCNK3 gene is the closest gene to the GGA3
response to E. maxima in broilers. genomic region, which was significantly associated with
BWG. The KCNK3 gene encodes a member of the potas-
GWAS sium channel superfamily, which has been associated
Previous QTL studies have reported that several genomic with pulmonary hypertension in humans [43]. Broil-
regions are associated with BWG in response to coccidi- ers are known to suffer from cardiovascular disorders
osis [15, 20]. Comparing our results with these previous [44, 45], which may be related to the KCNK3 gene. Fur-
studies, only one highly significant SNP that overlaps thermore, this gene may be associated with the general
with a QTL on GGA3 (between 263 and 282cM; 98.1 and robustness of broilers and their ability to cope with stress
107.0 Mb), detected by Pinard-van der Laan et al. [16] induced by the E. maxima challenge.
and Bacciu et al. [17], was observed. However, these dif- The THBS1 gene is located upstream of the signifi-
ferent results were not completely unexpected consider- cantly associated region GGA5 [between 28.95 and
ing that the previously conducted challenge study was 29.11 Mb, (See Additional file 6: Table S6)]. In addition,
performed with animals at a different age and that origi- we also analyzed putative GENSCAN Gene Predictions
nated from very different experimental lines. that are supported by a few spliced EST (See Additional
The MGAT4C gene is close to the significantly asso- file 12: Figure S5); however, no similarity was observed
ciated genomic region on GGA1. This gene encodes a with known proteins, which indicates that THBS1 was
glycosyltransferase that is involved in the transfer of a better candidate. Human THBS1 is involved in the
N-acetylglucosamine (GlcNAc) to the core mannose resi- regulation of angiogenesis and tumorigenesis in healthy
dues of N-linked glycans, also known as N-linked glyco- tissues and cell adhesion [46, 47]. Furthermore, in the
sylation. N-linked glycosylation has been shown to be study by Heams et al. [18], THBS1 was among the top
essential to HIV-1 pathogenesis [40]. Furthermore, there 10 of 1473 significantly differentially expressed genes in
is a wide range of well-described disorders that affect pri- the caecum between Fayoumi (resistant) and Leghorn
marily N-glycan assembly, with several including gastro- (susceptible) animals infected with E. tenella. Moreo-
intestinal disorders [41]. In addition, a study on humans ver, a porcine transcriptome study showed that THBS1
showed a strong association between the apolipoprotein is strongly repressed in the in vitro stimulation of por-
B level and SNP variants in the MGAT4C gene [42]. cine peripheral-blood mononuclear cells (PBMC) with
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 11 of 17

Fig. 8 Network of interactions between GWAS candidate genes using ingenuity pathway analysis (IPA). The network shows molecular interactions
between the products of the candidate genes enriched for significantly associated SNPs (P < 10−4) with myosin-6 (MYH6) as the most interesting
candidate. Relationships were determined using information contained in the IPA repository. The blue label indicates the genes that were enriched
for significantly associated SNPs. The network was obtained using the global gene list determined by combining gene lists for all traits

tetradecanoyl phorbol acetate (TPA)/ionomycin [48]. and clearance, receptor activation, signal transduction,
Finally, a recent study showed that human THBS1 plays and endocytosis [53]. In contrast, changes in the PC
an important role in the innate immune response dur- reflect the status of intestinal absorption, changes in the
ing respiratory bacterial infection [49], which may be production of protein carriers and their antioxidant effect
of interest regarding Eimeria infection. Based on these in response to Eimeria infection [54] and reflect several
observations, THBS1 is a good candidate gene for further of the processes that MAN2C1 may impact. Further-
functional studies. more, several significant glycol-related pathways were
Regarding PC, our GWAS results show little overlap significantly enriched in the biological pathway analysis
with previous QTL mapping studies [16]. We observed for PC measurements (Fig. 5).
the strongest signal with the SNP in the MAN2C1 gene Regarding the percentage of β-globulins, we identified
(See Additional file 6: Table S6) in association with PC FHOD3 and LARGE as potential candidate genes (See
measured in the range from 465 to 510 nm. The associ- Additional file 6: Table S6). FHOD3 encodes a protein
ated SNP is described as a non-synonymous polymor- that is a member of a formin subfamily and is involved
phism (http://www.ncbi.nlm.nih.gov/projects/SNP/ in the regulation of cell actin dynamics [55]. How-
snp_ref.cgi?rs=314637018). The MAN2C1 gene encodes ever, how FHOD3 may be involved in the regulation of
the α-mannosidase class 2C enzyme, which is one of plasma β-globulin levels is difficult to discern because
the key enzymes involved in N-glycan degradation [50]. the current knowledge regarding FHOD3 is scarce and
Recent studies have indicated that MAN2C1 expression restricted to human and mouse functional studies [56,
is crucial for maintaining efficient protein N-glycosyla- 57]. The LARGE gene encodes an enzyme glycosyltrans-
tion [51] as well as cell–cell adhesion [52]. Glycans are ferase that is involved in alpha-dystroglycan glycosyla-
important molecules in numerous essential biological tion and is capable of synthesizing glycoprotein and
processes, including cell adhesion, molecular trafficking glycosphingolipid sugar chains [58]. The exact function
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 12 of 17

of LARGE is not fully known; however, mutations in the Biological pathway analysis
human LARGE gene have been described to cause con- The biological pathway studies based on the GWAS
genital muscular dystrophy type 1D (MDC1D) [59]. results can potentially extend the knowledge obtained
Taking all measured traits into account, we identified from GWAS studies by identifying the cumulative effect
22 highly associated SNPs; however, GWAS was per- of gene sets [61]. Furthermore, understanding the biolog-
formed on two sets of traits: the first set of traits was ical pathways increases the power to detect statistically
measured on all 2024 animals, and the second set was significant associations because fewer statistical tests are
measured on the subset of 176 animals, for which the sta- performed as a consequence of assigning individual SNPs
tistical power to detect potential candidate genes differed to the respective biological pathways. Therefore, the
due to the different sample sizes. Regarding BWG and number of tests is decreased from over 400,000 (num-
PC that were measured on all animals, two functionally ber of SNPs) to approximately 160 (number of pathways)
well-supported candidate genes (THBS1 and MAN2C1) using a priori biological information. For this purpose, we
were detected. However, we did not identify any can- assigned genes, containing SNPs used in GWAS, to the
didate regions for BT and HEMA, potentially because biological KEGG pathways. Therefore, biological KEGG
these traits have not been as affected in challenged ani- pathways are considered as a set of genes that have been
mals compared with PC and BWG [19]. Furthermore, used for further analysis. However, we have to be aware
an infection caused by E. maxima is often characterized that the biological pathway analysis should still be pri-
by intestinal malabsorption and not by severe bloody marily viewed as an exploratory technique because the
diarrhoea, which is the case with E. tenella. This find- current statistical methodologies used for gene set/path-
ing further indicates that HEMA is not very informative way analysis need further development [61]. Moreover,
when measuring the response to E. maxima infection the available biological pathway databases necessary for
[19]. Similarly, BT is difficult to interpret with respect to this kind of analysis are not completely annotated and do
Eimeria infection because this trait may be influenced not contain all the genes present in the chicken genome.
by other factors as discussed by Hamzic et al. 2015 [19]. In this study, we used a statistical modeling methodol-
However, we have not identified many candidate genes ogy that assesses the cumulative effect of sets of SNPs on
or genomic regions in the GWAS performed by using the measured traits as presented by Jensen et al. [31] and
traits measured in the 176 animals, which may be pri- Buitenhuis et al. [62]. For this purpose, we succeeded in
marily due to the sample size, which was 10 times smaller assigning 52,204 SNPs from GWAS to 4342 annotated
in comparison to the traits measured on all animals, and genes in the chicken genome.
to the complexity of the genetic parameters that con- For PC, the top four most frequently affected pathways
trol these traits. In addition, this absence of identifiable across all wavelengths include phenylalanine, tyrosine
candidate genes may be partially due to the precision of and tryptophan biosynthesis, endocytosis, purine metab-
measured traits and rather stringent significance thresh- olism, and glycosphingolipid biosynthesis—lacto and
olds used in GWAS. Therefore, we performed biological neolacto series (Fig. 5).
pathway analysis, which enabled an increase in the statis- Phenylalanine, tyrosine and tryptophan biosynthesis
tical power to detect significant association. An increase is the most commonly affected pathway when only PC is
in the statistical power is possible because we decreased considered (Fig. 5) and (See Additional file 11: Figure S4).
the number of statistical tests performed by compressing Phenylalanine, tyrosine and tryptophan have important
individual SNPs into biological pathways [60]. roles in the regulation of the immune response [63]. Phe-
The final aim is to transfer the acquired knowledge to nylalanine is indirectly involved in the regulation of nitric
the poultry breeding industry. The identified candidate oxide (NO) synthesis [64], and NO is known to have mul-
genes and genomic regions can be used in breeding for tiple roles related to the immune response such as signal-
improved resistance to coccidiosis. The primary obstacle ing properties, regulating cytokine production and killing
in achieving this goal is relating the knowledge regard- pathogens [65]. Tyrosine is used as a precursor for the
ing the candidate genes identified in the final product production of dopamine, catecholamines and melanin.
(Cobb500) with the grand-parent pure lines where the Dopamine is a neurotransmitter known to be involved
actual selection is performed. This task becomes rather in the regulation of immune response, and melanin
complex, considering that the poultry populations and has antioxidant properties [63]. In addition, interferon
the crossing design are not in the public domain, and due gamma (IFN-γ) suppresses the growth of Toxoplasma
to the proprietary nature of information regarding the gondii through the intracellular depletion of tryptophan
grandparent lines. Future studies will identify the best [66]. Deprivation of tryptophan produces a deleterious
approach to trace the identified regions to the grandpar- effect on Toxoplasma gondii replication. Toxoplasma
ent lines. gondii belongs to the same order (Eucoccidiorida) of
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 13 of 17

intracellular single-cell parasites as E. maxima. Fur- in resilience and immunity to infectious diseases [73].
thermore, Laurent et al. [67] reported a strong increase Moreover, the peroxisome pathway was significantly
in IFN-γ mRNA expression in chickens infected with associated with blood components such as number of
Eimeria spp. Based on this finding and the results from erythrocytes and percentage of heterophils and lympho-
the biological pathway analysis, tryptophan depletion cytes (Fig. 6).
may also be involved in the innate immune response dur- These findings demonstrate that the response to
ing E. maxima infection. Eimeria infection is characterized by a strong effect on
The second most commonly affected pathway when essential metabolic pathways as well as innate immune
considering PC is the purine metabolism pathway (Fig. 5) response-related pathways. Among the essential meta-
and (See Additional file 11: Figure S4). This pathway reg- bolic pathways, the most frequently affected are pheny-
ulates nucleotide metabolism and is important for suc- lalanine, tyrosine and tryptophan biosynthesis, purine
cessful cell division. In addition, we observed that the metabolism and the pentose phosphate pathway. In addi-
purine metabolism pathway is significantly associated tion, the most frequent innate immune response-related
with erythrocyte number, mean corpuscular hemoglobin pathways are glycol-related pathways, the peroxisome
(MCH), mean corpuscular hemoglobin concentration pathway and the endocytosis pathway.
(MCHC) and mean cellular volume (MCV) (Fig. 6). This
association may be explained by an increased demand Network‑based analysis
for cell division of blood cell progenitors and the regen- The network-based analysis was performed as a com-
eration of the intestinal epithelium due to the effects of plementary approach to the biological pathway analysis
the infection because the E. maxima infection has been to build gene networks associated with responses to the
characterized by severe lesions of the intestinal lining in E. maxima challenge using bibliography-based proven
broilers [19]. relationships that are available through the Ingenuity
Endocytosis has been shown to play an important role Pathway repository. The network-based analysis was per-
in both innate and adaptive immune responses [68], formed independently from the biological pathway analy-
which may explain why endocytosis may be among the sis and was based on a list of genes enriched for SNPs that
most commonly affected pathways when PC measure- are associated (p < 10−4) with the traits measured during
ments are considered (Fig. 5). Finally, glycosphingolipid the E. maxima challenge (See Additional file 4: Table S4).
biosynthesis—lacto and neolacto series, like other glycol- The network-based analysis approach, implemented in
related pathways (Fig. 5), is involved in the production the IPA tool, assumes that genes used as an input interact
of glycoconjugate receptors, which are used by microbes with each other, and these interactions are reconstructed
to enter the host cell and are of critical importance in based on the relationships shown in the literature [34].
the early stage of the innate immune response [69, 70]. The network-based analysis aims at exploring the cumu-
Moreover, similar paths of the host cell invasion have lative effect of sets of genes that individually explain a
been previously described in the cases of Eimeria and moderate part of the variation for a measured trait and
Toxoplasma [71, 72]. that cannot be identified during GWAS when the strict
We also summarized the frequency of the most com- significance threshold was applied.
mon significant biological pathways excluding PC (Fig. 6). Figure 7 illustrates the network formed by several of the
In this case, glycosaminoglycan biosynthesis involving molecules grouped around LATS1 and LATS2 with mul-
heparan sulfate/heparin, the pentose phosphate pathway tiple direct connections with other molecules enriched
and the peroxisome are the most frequent significantly with significantly associated SNPs. LATS1 and LATS2
enriched biological pathways. The pentose phosphate are known to be involved in the regulation of intesti-
pathway is a metabolic pathway that regulates the pro- nal epithelium renewal [74], which may be explained by
duction of nicotinamide adenine dinucleotide phos- intensified tissue repair upon E. maxima challenge. In
phate (NADPH) and pentose, which are essential for the addition, phenylalanine, tyrosine and tryptophan bio-
synthesis of nucleic and ribonucleic acids, respectively. synthesis, purine metabolism and the pentose phosphate
The pentose phosphate pathway seems to play a role in pathway are the most frequent significant pathways asso-
plasma protein component and heterophil and lympho- ciated with all measured traits in the KEGG biological
cyte production, which may be explained as a response pathway analysis. All three of these KEGG pathways are
to the increased production of these components during associated with increased DNA replication, cell metabo-
the primary response to the challenge (Fig. 6) [73]. The lism and protein degradation, which are essential during
peroxisome pathway controls the metabolism of reac- the tissue repair process.
tive oxygen components, which are known to be toxic to The second network has the MYH6 as a key molecule
bacteria and several parasites and play a significant role connected with several direct significant relationships to
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 14 of 17

the other genes that are enriched for significantly asso- The post-GWAS functional analysis, which combined
ciated SNPs (Fig. 8). The MYH6 gene encodes the alpha two independent approaches (the biological pathway
heavy chain subunit of cardiac myosin. In mice, inactiva- analysis and network-based analysis), indicated that the
tion of the specific mutant MYH6 transcript suppresses genes and biological pathways involved in tissue repair,
hypertrophic cardiomyopathy [75]. In addition, we also general robustness as well as the primary immune
identified KCKN3, which is associated with pulmonary response may play an important role during the primary
hypertension, as one of the candidate genes because stage of E. maxima infection in broilers.
heart failure and ascites have been well documented in Studies that focus on the transfer of the acquired
broiler chickens [44, 45]. The primary reason for these knowledge to poultry breeding considering the specifici-
problems can be attributed to an intensive selection in ties of the broiler breeding scheme are currently under
poultry breeding during the last 60 years [76]. Therefore, way. Finally, a follow-up transcriptome study is ongoing,
we can potentially indicate which animals are able to which aims at integrating results from the GWAS study
maintain a normal function of the cardiovascular system and further investigation of the genetic mechanisms that
and have an advantage in the face of Eimeria infection. control the response to Eimeria infection in broilers.
Based on the post-GWAS functional analysis, the
broiler response to E. maxima is centered on tissue Availability of supporting data
repair and recovery, general robustness and mainte- No new SNPs were discovered in the preparation of this

the chicken 580 K Affymetrix® Axiom® HD genotyping

nance of tissue integrity, restoring intestinal homeostasis manuscript. The SNPs used in this manuscript are from
after the challenge. The described processes, which may
bring a comparative advantage in the broiler’s ability to array: http://media.affymetrix.com/support/technical/
cope with the challenge, can be described as resilience datasheets/axiom_chicken_array_plate_datasheet.pdf.
to acute Eimeria infection. In contrast to the previously SNP names and location can be found at http://www.
conducted studies [16, 17], which reported associations affymetrix.com/catalog/prod670010/AFFY/Axiom%26%
with genes involved in the immune response, we primar- 23174%3B+Genome%26%2345%3BWide+Chicken+Ge
ily observed associations with genes, biological pathways notyping+Array#1_3.
and gene networks that are involved in tissue repair and
recovery and tissue integrity maintenance. However, pre-
vious studies [16, 17] were conducted on experimental Additional files
layer populations challenged with E. tenella at 28 days of
age. In regard to this, broilers are also able to establish Additional file 1: Table S1. Title: Quality control of genotype data.
complete immunity 16 days after being challenged with Description: The genotype data were obtained from 1972 animals

produced using 580 K Affymetrix® Axiom® HD genotyping array. The call

that were alive at day 16 of the challenge study. Genotype data were
E. maxima [77]. Therefore, one would assume that this
challenge study would identify genetic variants associ- rate threshold was set at 98 %, and all SNPs and samples with values less
ated with processes related to immune responses, which than 98 % were excluded. The minor allele frequency threshold was set at
1 % and animals with values less than 1 % were excluded. Heterozygote
did not occur in this study and may be due to the effect excess was assessed based on the distribution of F-statistics adjusted for
of the infection doses (50,000 oocysts), which were opti- skewness using adjbox function from the R package robustbase. Control
mized to produce severe clinical signs as reported by animals were also excluded from the analysis. A total of 138,568 SNPs out
of 580,961 were removed during the quality control.
Hamzic et al. [19]. We argue that the more resilient ani-
mals are able to maintain their biological homeostasis Additional file 2: Table S2. Title: Descriptive statistics for the distances
between SNPs across chromosomes. Description: The table contains infor‑
and manage the consequences of the infection, which, in mation on the number of SNPs per chromosome and the mean, median
the context of this challenge, exceeded the importance of and standard deviation (SD) for the average distances between neighbor‑
building an adequate immune response. ing SNPs. Results in the table refer to genotype data after quality control.
Additional file 3: Table S3. Title: List of annotated SNPs used for
Biological Pathway Analysis. Description: The table contains SNPs used
Conclusions in the genome-wide association study with p-values < 10−4 with the
We identified 22 SNPs significantly associated with four corresponding information on the genes in which they are located.
different traits at q value <0.1. Two candidate genes, Abbreviations used for column names: RS ID: Reference SNP ID number,
SYMBOL: HGNC symbol for gene, NCBI GENE ID: NCBI gene ID, CHR: chro‑
MAN2C1 and FHOD3, were significantly associated with mosome, CHR POS: chromosome position, GENE REGION: gene region of
PC measured in the range from 465 to 510 nm and the the corresponding gene, SS ID: submitted SNP ID, AFFYMETRIX: Affymetrix
percentage of β2-globulin in blood plasma, respectively. SNP ID, ALLELE: SNP alleles. The remaining columns are p-values for each
trait with column names corresponding to trait abbreviations used in the
Moreover, we identified three genomic regions on GGA1 manuscript.
(MGAT4C), GGA3 (KCNK3) and GGA5 (THBS1) that Additional file 4: Table S4. Title: List of annotated SNPs and their cor‑
are significantly associated with body weight gain and the responding genes used for Ingenuity Pathway Analysis. Description: The
percentage of β2-globulin. table contains SNPs with the corresponding information on the genes
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 15 of 17

FH was involved in the early stages of data analysis and helped draft the
in which they are located and together with p-values obtained in the manuscript. BS and JME helped with the GWAS methodology selection. MHP
genome-wide association study. designed the study. BBH participated in the study design, supervised the
Additional file 5: Table S5. Title: Descriptive statistics of the measured analysis with emphasis on the functional section and helped draft the manu‑
traits. Description: Descriptive statistics of the traits measured on the script. All authors read and approved the final manuscript.
challenged animals during the large-scale challenge study. Plasma colora‑
tion was measured as absorbance (optical density) of blood plasma at Author details
1
different wavelengths. Oocyst count is expressed as the raw number of UMR1313 Animal Genetics and Integrative Biology Unit, AgroParisTech,
oocysts per gram of feces. Erythrocyte count is expressed in millions per 16 rue Claude Bernard, 75005 Paris, France. 2 UMR1313 Animal Genetics
mm3 (M∙mm−3). White blood cell and thrombocyte count is expressed in and Integrative Biology Unit, INRA, Domaine de Vilvert, 78350 Jouy‑en‑Josas,
absolute number per mm3 (count∙mm−3. Lymphocytes and heterophils France. 3 Department of Molecular Biology and Genetics, Center for Quantita‑
are also expressed in the percentage of total number of white blood tive Genetics and Genomics, Aarhus University, Blichers Allé 20, P.O. Box 50,
cells. MCV, MCH and MCHC correspond to mean cellular volume, mean 8830 Tjele, Denmark. 4 UMR1348 Physiology, Environment and Genetics for the
corpuscular hemoglobin and mean corpuscular hemoglobin concentra‑ Animal and Livestock Systems Unit, INRA, Domaine de la Prise, 35590 Saint
tion, respectively. Gilles, France. 5 Cobb-Vantress Inc., Siloam Springs, AR 72761, USA. 6 UMR1388
Genetics, Physiology and Breeding Systems, INRA, 24 chemin de Borde‑Rouge,
Additional file 6: Table S6. Title: Significant SNPs associated with traits
31326 Castanet‑Tolosan, France.
measured on all challenged animals. Description: Significant SNPs associ‑
ated with traits measured on all challenged animals. A genome-wide
Acknowledgements
association analysis was performed using Cobb500 broilers for all the
EH benefited from a joint grant from the European Commission and Cobb-
measured traits during the large-scale study. SNPs were considered signifi‑
Vantress Inc. within the framework of the Erasmus-Mundus joint doctorate
cant at q-value  < 0.10 using False Discovery Rate adjustment.
“EGS-ABG”. BB was supported by the “Development of genetic selection
Additional file 7: Figure S1. Title: Quantile–quantile plot for body technology for polyvalent resistance to infectious diseases” (POLY-REID) (Grant
weight gain. Description: Quantile–quantile plot for the body weight gain Number 10-093534) project granted by the Danish Council for Strategic
test statistics. Research, the Danish Poultry Council, The Hatchery Hellevad, and Cobb-
Vantress Inc.
Additional file 8: Figure S2. Title: Quantile–quantile plot for plasma
coloration (485 nm). Description: Quantile–quantile plot for plasma
Competing interests
coloration (485 nm) test statistics.
The authors declare that they have no competing interests.
Additional file 9: Figure S3. Title: Quantile–quantile plot for the
percentage of β2-globulin. Description: Quantile–quantile plot for Received: 23 May 2015 Accepted: 5 November 2015
β2-globulin test statistics.
Additional file 10: Table S7. Title: List of biological pathways signifi‑
cantly enriched with genes in genomic regions associated with measured
traits. Description: Results for the biological pathway significantly
(P < 0.05) enriched with genes in genomic regions associated with meas‑
References
ured traits during the large-scale challenge study. The p-values provided
1. Smith AL, Hesketh P, Archer A, Shirley MW. Antigenic diversity in Eimeria
in the table represent values for the one-sided test for overlapping SNPs
maxima and the influence of host genetics and immunization schedule
after 1000 permutations.
on cross-protective immunity. Infect Immun. 2002;70:2472–9.
Additional file 11: Figure S4. Title: Distribution of the most frequent 2. Conway DP, McKenwie EM. Poultry coccidiosis diagnostic and testing
significant biological pathways. Description: Distribution of the most procedures. Ames: Blackwell Publishing Professional; 2007.
frequent biological pathways being significantly enriched with genes in 3. Dalloul RA, Lillehoj HS. Poultry coccidiosis: recent advancements in
genomic regions associated with measured traits during the large-scale control measures and vaccine development. Expert Rev Vaccines.
challenge study. The distribution is considerably different when all 2006;5:143–63.
measured traits are considered together and when PC measurements are 4. Williams RB. A compartmentalised model for the estimation of the cost
excluded. of coccidiosis to the world’s chicken production industry. Int J Parasitol.
1999;29:1209–29.
Additional file 12: Figure S5. Title: Illustration of the genomic region
5. Blake DP, Tomley FM. Securing poultry production from the ever-present
between 28.95 and 29.11 Mb on GGA5. Description: The genomic region
Eimeria challenge. Trends Parasitol. 2014;30:12–9.
between 28.95 and 29.11 Mb contains three SNPs that were significantly
6. Jeffers TK. Anticoccidials: past. Feedstuffs. 2011;11:16.
associated with body weight gain (BWG) and also contains several
7. Jeffers TK. Eimeria acervulina and E. maxima: incidence and anticoccidial
chicken spliced mRNA and EST, which may indicate the presence of a
drug resistance of isolants in major broiler-producing areas. Avian Dis.
non-annotated chicken gene. Illustration was obtained from the UCSC
1974;18:331–42.
Genome Browser (ICGSC Gallus_gallus-4.0/galGal4).
8. Lamont SJ. Impact of genetics on disease resistance. Poult Sci.
1998;77:1111–8.
Abbreviations 9. Rosenberg MM. A study of the inheritance of resistance to Eimeria tenella
BWG: body weight gain; HEMA: hematocrit; LS: lesion score; OC: oocyst in the domestic fowl. Poult Sci. 1941;20:472.
count; PC: plasma coloration; GWAS: genome-wide association study; QTL: 10. Edgar SA. Control of avian coccidiosis through breeding or immunization.
quantitative trait loci; PPP: plasma protein profiles; BC: blood composition; NO: Poult Sci. 1951;30:911.
nitric oxide; IFN-γ: interferon gamma; EST: expressed sequence tags; PBMC: 11. Long PL. The effect of breed of chickens on resistance to Eimeria infec‑
peripheral blood mononuclear cell; IPA: Ingenuity Pathway Analysis; HGNC: tions. Br Poult Sci. 1968;9:71–8.
Human Genome Gene Nomenclature Committee; TPA: tetradecanoyl phorbol 12. Johnson LW, Edgar SA. Responses to prolonged selection for resistance
acetate; MDC1D: congenital muscular dystrophy type 1D; MCH: mean corpus‑ and susceptibility to acute cecal coccidiosis in the Auburn strain single
cular hemoglobin; MCV: mean corpuscular volume; MCHC: mean corpuscular comb white Leghorn. Poult Sci. 1982;61:2344–55.
hemoglobin concentration; HEMA: hematocrit level; BT: body temperature; 13. Bumstead N, Millard B. Genetics of resistance to coccidiosis: response of
BW: body weight. inbred chicken lines to infection by Eimeria tenella and Eimeria maxima.
Br Poult Sci. 1987;28:705–15.
Authors’ contributions 14. Pinard-van der Laan MH, Monvoisin JL, Pery P, Hamet N, Thomas M.
EH performed the data analysis and wrote the manuscript. BB supervised the Comparison of outbred lines of chickens for resistance to experimental
biological pathway and GWA analysis and helped draft of the manuscript. infection with coccidiosis (Eimeria tenella). Poult Sci. 1998;77:185–91.
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 16 of 17

15. Zhu JJ, Lillehoj HS, Allen PC, Van Tassell CP, Sonstegard TS, Cheng HH, 39. Robinson MR, Wray NR, Visscher PM. Explaining additional genetic varia‑
et al. Mapping quantitative trait loci associated with resistance to coccidi‑ tion in complex traits. Trends Genet. 2014;30:124–32.
osis and growth phenotypic measurement of disease susceptibility. Poult 40. Montefiori DC, Robinson WE, Mitchell WM. Role of protein N-glycosyla‑
Sci. 2003;82:9–16. tion in pathogenesis of human immunodeficiency virus type 1. Proc Natl
16. Pinard-van der Laan MH, Bed’Hom B, Coville JL, Pitel F, Feve K, Leroux S, Acad Sci USA. 1988;85:9248–52.
et al. Microsatellite mapping of QTLs affecting resistance to coccidiosis 41. Leroy JG. Congenital disorders of N-glycosylation including diseases
(Eimeria tenella) in a Fayoumi × White Leghorn cross. BMC Genomics. associated with O- as well as N-glycosylation defects. Pediatr Res.
2009;10:31. 2006;60:643–56.
17. Bacciu N, Bed’Hom B, Filangi O, Romé H, Gourichon D, Répérant JM, et al. 42. Kathiresan S, Manning AK, Demissie S, D’Agostino RB, Surti A, Guiducci C,
QTL detection for coccidiosis (Eimeria tenella) resistance in a Fayoumi × et al. A genome-wide association study for blood lipid phenotypes in the
Leghorn F2 cross, using a medium-density SNP panel. Genet Sel Evol. Framingham Heart Study. BMC Med Genet. 2007;8:S17.
2014;46:14. 43. Girerd B, Perros F, Antigny F, Humbert M, Montani D. KCNK3: new gene
18. Heams T, Bed’Hom B, Rebours E, Jaffrezic F, Pinard-van der Laan M-H. target for pulmonary hypertension? Expert Rev Respir Med. 2014;8:385–7.
Insights into gene expression profiling of natural resistance to coccidiosis 44. Baghbanzadeh A, Decuypere E. Ascites syndrome in broilers: physiologi‑
in contrasting chicken lines. BMC Proc. 2011;5:S26. cal and nutritional perspectives. Avian Pathol. 2008;37:117–26.
19. Hamzic E, Bed’Hom B, Juin H, Hawken R, Abrahamsen MS, Elsen JM, 45. Olkowski AA. Pathophysiology of heart failure in broiler chickens: structural,
et al. Large-scale investigation of the parameters in response to Eimeria biochemical, and molecular characteristics. Poult Sci. 2007;86:999–1005.
maxima challenge in broilers. J Anim Sci. 2015;93:1830–40. 46. Ernens I, Bousquenaud M, Lenoir B, Devaux Y, Wagner DR. Adenosine
20. Kranis A, Gheyas AA, Boschiero C, Turner F, Yu L, Smith S, et al. Develop‑ stimulates angiogenesis by up-regulating production of thrombospon‑
ment of a high density 600 K SNP genotyping array for chicken. BMC din-1 by macrophages. J Leukoc Biol. 2015;97:9–18.
Genomics. 2013;14:59. 47. Leclair P, Lim CJ. CD47-independent effects mediated by the TSP-derived
21. Chang CC, Chow CC, Tellier LC, Vattikuti S, Purcell SM, Lee JJ. Second- 4N1K peptide. PLoS One. 2014;9:e98358.
generation PLINK: rising to the challenge of larger and richer datasets. 48. Gao Y, Flori L, Lecardonnel J, Esquerré D, Hu ZL, Teillaud A, et al. Transcrip‑
Gigascience. 2015;4:7. tome analysis of porcine PBMCs after in vitro stimulation by LPS or PMA/
22. Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. ionomycin using an expression array targeting the pig immune response.
PLINK: a tool set for whole-genome association and population-based BMC Genomics. 2010;11:292.
linkage analyses. Am J Hum Genet. 2007;81:559–75. 49. Zhao Y, Olonisakin TF, Xiong Z, Hulver M, Sayeed S, Yu MT, et al. Throm‑
23. PLINK Version 1.90b3. 2015. https://www.cog-genomics.org/plink2. bospondin-1 restrains neutrophil granule serine protease function and
Accessed 23 April 2015. regulates the innate immune response during Klebsiella pneumoniae
24. Hubert M, Vandervieren E. An adjusted boxplot for skewed distributions. infection. Mucosal Immunol. 2015;8:896–905.
Comput Stat Data Anal. 2008;52:5186–201. 50. Kuokkanen E, Smith W, Mäkinen M, Tuominen H, Puhka M, Jokitalo E, et al.
25. Zhou X, Stephens M. Efficient algorithms for multivariate linear mixed Characterization and subcellular localization of human neutral class II
model algorithms in genome-wide association studies. Nat Methods. alpha-mannosidase [corrected]. Glycobiology. 2007;17:1084–93.
2014;11:407–9. 51. Bernon C, Carré Y, Kuokkanen E, Slomianny MC, Mir AM, Krzewinski F,
26. Yu J, Pressoir G, Briggs WH, Vroh Bi I, Yamasaki M, Doebley JF, et al. A et al. Overexpression of Man2C1 leads to protein underglycosylation and
unified mixed-model method for association mapping that accounts for upregulation of endoplasmic reticulum-associated degradation pathway.
multiple levels of relatedness. Nat Genet. 2006;38:203–8. Glycobiology. 2011;21:363–75.
27. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a 52. Qu L, Ju JY, Chen SL, Shi Y, Xiang ZG, Zhou YQ, et al. Inhibition of the
practical and powerful approach to multiple testing. J R Stat Soc B. alpha-mannosidase Man2c1 gene expression enhances adhesion of
1995;57:289–300. Jurkat cells. Cell Res. 2006;16:622–31.
28. Turner SD. qqman: an R package for visualizing GWAS results using Q-Q 53. Ohtsubo K, Marth JD. Glycosylation in cellular mechanisms of health and
and Manhattan plots. bioRxiv. 2014. doi:10.1101/005165. disease. Cell. 2006;126:855–67.
29. Melville S. NCBI2R-An R package to navigate and annotate genes and 54. Yvore P, Mancassola R, Naciri M, Bessay M. Serum coloration as a criterion of the
SNPs. R package version 1.4.7. 2015. http://cran.rproject.org/src/contrib/ severity of experimental coccidiosis in the chicken. Vet Res. 1993;24:286–90.
Archive/NCBI2R/. Accessed 23 April 2015. 55. Taniguchi K, Takeya R, Suetsugu S, Kan-O M, Narusawa M, Shiose A,
30. Tenenbaum D. KEGGREST: Client-side REST access to KEGG. R package et al. Mammalian formin fhod3 regulates actin assembly and sarcomere
version 1.8.0. 2014. http://bioconductor.org/packages/release/bioc/html/ organization in striated muscles. J Biol Chem. 2009;284:29873–81.
KEGGREST.html. Accessed 23 April 2015. 56. Wooten EC, Hebl VB, Wolf MJ, Greytak SR, Orr NM, Draper I, et al. Formin
31. Rohde PD, McKinnon-Edwards S, Sarup PM, Sorensen P. Gene based homology 2 domain containing 3 variants associated with hypertrophic
association approach identify genes across stress traits in fruit flies. In cardiomyopathy. Circ Cardiovasc Genet. 2013;6:10–8.
10th world congress on genetics applied to livestock production: 17–22 57. Kan-O M, Takeya R, Abe T, Kitajima N, Nishida M, Tominaga R, et al. Mam‑
August 2014; Vancouver. https://asas.org/docs/default-source/wcgalp- malian formin Fhod3 plays an essential role in cardiogenesis by organ‑
posters/673_paper_9223_manuscript_493_0.pdf. izing myofibrillogenesis. Biol Open. 2012;1:889–96.
32. Newton MA, Quintana FA, den Boon JA, Sengupta S, Ahlquist P. Random- 58. Peyrard M, Seroussi E, Sandberg-Nordqvist AC, Xie YG, Han FY, Fransson
set methods identify distinct aspects of the enrichment signal in gene- I, et al. The human LARGE gene from 22q12.3-q13.1 is a new, distinct
set analysis. Ann Appl Stat. 2007;1:85–106. member of the glycosyltransferase gene family. Proc Natl Acad Sci USA.
33. Jiang Z, Gentleman R. Extensions to gene set enrichment. Bioinformatics. 1999;96:598–603.
2007;23:306–13. 59. Longman C, Brockington M, Torelli S, Jimenez-Mallebrera C, Kennedy C,
34. Krämer A, Green J, Pollard J Jr, Tugendreich S. Causal analysis approaches Khalil N, et al. Mutations in the human LARGE gene cause MDC1D, a novel
in ingenuity pathway analysis. Bioinformatics. 2014;30:523–30. form of congenital muscular dystrophy with severe mental retardation
35. Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the and abnormal glycosylation of alpha-dystroglycan. Hum Mol Genet.
integration of genomic datasets with the R/Bioconductor package 2003;12:2853–61.
biomaRt. Nat Protoc. 2009;4:1184–91. 60. Fridley BL, Biernacka JM. Gene set analysis of SNP data: benefits, chal‑
36. Kim ES, Hong YH, Min W, Lillehoj HS. Fine-mapping of coccidia-resistant lenges, and future directions. Eur J Hum Genet. 2011;19:837–43.
quantitative trait loci in chickens. Poult Sci. 2006;85:2028–30. 61. Mooney MA, Nigg JT, McWeeney SK, Wilmot B. Functional and
37. Biscarini F, Bovenhuis H, van Arendonk JAM, Parmentier HK, Jungerius AP, genomic context in pathway analysis of GWAS data. Trends Genet.
van der Poel JJ. Across-line SNP association study of innate and adaptive 2014;30:390–400.
immune response in laying hens. Anim Genet. 2010;41:26–38. 62. Buitenhuis B, Janss LLG, Poulsen NA, Larsen LB, Larsen MK, Sørensen P.
38. Ellen ED, Visscher J, van Arendonk JAM, Bijma P. Survival of laying hens: Genome-wide association and biological pathway analysis for milk-fat
genetic parameters for direct and associative effects in three purebred composition in Danish Holstein and Danish Jersey cattle. BMC Genomics.
layer lines. Poult Sci. 2008;87:233–9. 2014;15:1112.
Hamzić et al. Genet Sel Evol (2015) 47:91 Page 17 of 17

63. Li P, Yin YL, Li D, Kim SW, Wu G. Amino acids and immune function. Br J 71. Lai L, Bumstead J, Liu Y, Garnett J, Campanero-Rhodes MA, Blake DP, et al.
Nutr. 2007;98:237–52. The role of sialyl glycan recognition in host tissue tropism of the avian
64. Shi W, Meininger CJ, Haynes TE, Hatakeyama K, Wu G. Regulation of parasite Eimeria tenella. PLoS Pathog. 2011;7:e1002296.
tetrahydrobiopterin synthesis and bioavailability in endothelial cells. Cell 72. Friedrich N, Santos JM, Liu Y, Palma AS, Leon E, Saouros S, et al. Members
Biochem Biophys. 2004;41:415–34. of a novel protein family containing microneme adhesive repeat domains
65. Bogdan C. Nitric oxide synthase in innate and adaptive immunity: an act as sialic acid-binding lectins during host cell invasion by apicompl‑
update. Trends Immunol. 2015;36:161–78. exan parasites. J Biol Chem. 2010;285:2064–76.
66. Pfefferkorn ER. Interferon gamma blocks the growth of Toxoplasma gondii 73. Allen PC, Fetterer RH. Recent advances in biology and immunobiology
in human fibroblasts by inducing the host cells to degrade tryptophan. of Eimeria species and in diagnosis and control of infection with these
Proc Natl Acad Sci USA. 1984;81:908–12. coccidian parasites of poultry. Clin Microbiol Rev. 2002;15:58–65.
67. Laurent F, Mancassola R, Lacroix S, Menezes R, Naciri M. Analysis of 74. Imajo M, Ebisuya M, Nishida E. Dual role of YAP and TAZ in renewal of the
chicken mucosal immune response to Eimeria tenella and Eimeria max- intestinal epithelium. Nat Cell Biol. 2015;17:7–19.
ima infection by quantitative reverse transcription-PCR. Infect Immun. 75. Jiang J, Wakimoto H, Seidman JG, Seidman CE. Allele-specific silencing of
2001;69:2527–34. mutant Myh6 transcripts in mice suppresses hypertrophic cardiomyopa‑
68. Gleeson PA. The role of endosomes in innate and adaptive immunity. thy. Science. 2013;342:111–4.
Semin Cell Dev Biol. 2014;31:64–72. 76. Havenstein G, Ferket P, Qureshi M. Growth, livability, and feed conversion
69. Svensson M, Platt FM, Svanborg C. Glycolipid receptor depletion as of 1957 versus 2001 broilers when fed representative 1957 and 2001
an approach to specific antimicrobial therapy. FEMS Microbiol Lett. broiler diets. Poult Sci. 2003;82:1500–8.
2006;258:1–8. 77. Stiff MI, Bafundo KW. Development of immunity in broilers continuously
70. Lingwood CA. Glycosphingolipid functions. Cold Spring Harb Perspect exposed to Eimeria sp. Avian Dis. 1993;37:295–301.
Biol. 2011;3:a004788.

Submit your next manuscript to BioMed Central

and take full advantage of:

• Convenient online submission

• Thorough peer review
• No space constraints or color figure charges
• Immediate publication on acceptance
• Inclusion in PubMed, CAS, Scopus and Google Scholar
• Research which is freely available for redistribution

Submit your manuscript at

www.biomedcentral.com/submit