Skip to main content

Detection and replication of QTL underlying resistance to gastrointestinal nematodes in adult sheep using the ovine 50K SNP array



Persistence of gastrointestinal nematode (GIN) infection and the related control methods have major impacts on the sheep industry worldwide. Based on the information generated with the Illumina OvineSNP50 BeadChip (50 K chip), this study aims at confirming quantitative trait loci (QTL) that were previously identified by microsatellite-based genome scans and identifying new QTL and allelic variants that are associated with indicator traits of parasite resistance in adult sheep. We used a commercial half-sib population of 518 Spanish Churra ewes with available data for fecal egg counts (FEC) and serum levels of immunoglobulin A (IgA) to perform different genome scan QTL mapping analyses based on classical linkage analysis (LA), a combined linkage disequilibrium and linkage analysis (LDLA) and a genome-wide association study (GWAS).


For the FEC and IgA traits, we detected a total of three 5 % chromosome-wise significant QTL by LA and 63 significant regions by LDLA, of which 13 reached the 5 % genome-wise significance level. The GWAS also revealed 10 significant SNPs associated with IgAt, although no significant associations were found for LFEC. Some of the significant QTL for LFEC that were detected by LA and LDLA on OAR6 overlapped with a highly significant QTL that was previously detected in a different half-sib population of Churra sheep. In addition, several new QTL and SNP associations were identified, some of which show correspondence with effects that were reported for different populations of young sheep. Other significant associations that did not coincide with previously reported associations could be related to the specific immune response of adult animals.


Our results replicate a FEC-related QTL located on OAR6 that was previously reported in Churra sheep and provide support for future research on the identification of the allelic variant that underlies this QTL. The small proportion of genetic variance explained by the detected QTL and the large number of functional candidate genes identified here are consistent with the hypothesis that GIN resistance/susceptibility is a complex trait that is not determined by individual genes acting alone but rather by complex multi-gene interactions. Future studies that combine genomic variation analysis and functional genomic information may help elucidate the biology of GIN disease resistance in sheep.


Persistence of gastrointestinal nematode (GIN) infection and the related control methods have major impacts on the sheep industry worldwide [1]. The extensive use of anthelmintics has negative consequences, such as the costs of treatments, the emergence of anthelmintic-resistant strains of parasites, and the presence of drug residues in animal products. Among different alternatives to chemical control, the selection of genetically-resistant animals has been suggested to reduce dependence on the use of anthelmintics [2, 3]. Selective breeding for resistance to GIN using fecal egg count (FEC) as an indicator trait has been undertaken for certain sheep breeds [46]. However, classical selection for this complex phenotype is hindered by the time-consuming and costly process of recording information for indicator phenotypes (which may also include serum levels of e.g., immunoglobulin A (IgA), IgE and pepsinogen) and by the requirement for animals to be infected by GIN at sampling. These difficulties suggest that selecting animals resistant to GIN infection would be more efficient if it was based on indirect estimates, such as those generated from molecular marker information. In the last few decades, considerable effort has been made to understand the relationship between host and parasite and the mechanisms that underlie host resistance [7]. Moreover, the recent availability of the Illumina OvineSNP50 BeadChip (Illumina Inc., San Diego, CA) (referred to here as the “50 K chip”) and a high-quality reference genome assembly [8] may allow for a deeper understanding of the genetic architecture of complex traits in sheep. Effective exploitation of this molecular information will increase our chances of developing protocols that will enable efficient selection of animals with increased resistance to GIN infections.

Because GIN are particularly pathogenic to young naïve animals such as growing lambs, gastrointestinal infections constitute a major cost to the sheep meat industry [9]. Accordingly, most of the quantitative trait locus (QTL) studies on GIN resistance traits [10], including those based on microsatellite markers as well as more recent analyses that exploit the ovine 50 K chip, have been conducted primarily on young animals [1126]. Conversely, for the Mediterranean dairy sheep industry, a production system that is based on adult ewes and the sale of suckling lambs fed exclusively on maternal milk, replacement ewes and adult sheep are the only animals subjected to the direct effects of helminth infections [27]. In these animals, the breakdown of the acquired immunity to infection that occurs around the time of parturition [28] and the necessity of anthelmintic treatment determine how severe the economic losses will be [29].

Previously, we performed a genome scan using microsatellite markers to identify QTL that influence indicator traits of parasite resistance in adult Churra dairy sheep, an autochthonous dairy breed of the northwest region of Castilla y León in Spain [20]. The lack of strong coincidence between the QTL that we had identified and those previously detected by using lamb data suggested that aside from differences in host-parasite combinations, these QTL could be related to different mechanisms that underlie resistance between adult sheep and lambs.

Within this context, we undertook a new QTL mapping study based on the use of the ovine 50 K chip to genotype a commercial population of Spanish Churra dairy sheep. To follow on the initial linkage analysis-based genome scan reported by Gutiérrez-Gil et al. [20], our study was designed to replicate some of the QTL that were detected by the microsatellite-based scan and to identify new QTL and allelic variants associated with two previously analyzed indicator traits of parasite resistance: FEC and serum levels of IgA. For this purpose, we performed the new analyses using a different set of half-sib families from the same commercial population of Spanish Churra sheep. Taking advantage of the increased marker density offered by the 50 K chip, in addition to classical linkage analysis (LA), we also implemented combined linkage disequilibrium and linkage analysis (LDLA) and genome-wide association study (GWAS) approaches to provide a more complete picture of the QTL that segregate in this ovine population.


Resource population and sampling

Phenotypic and genotypic information for 518 Churra ewes from the Selection Nucleus of the National Association of Churra Breeders (ANCHE) was analyzed. The animals belonged to 14 half-sib families and were produced by artificial insemination, with an average family size of 37 daughters per sire (ranging from 12 to 89). A single collection of fecal and blood samples was performed for each of the 17 flocks in the Castilla y León region where the animals were raised. The samples were later processed to measure two indicator traits of parasite resistance, FEC and serum IgA levels. The ages of the sheep included in this study ranged from 4 to 11 years. At the time of sampling, all the sheep were undergoing milking and were at least in their third lactation.

Phenotypic records

FEC measurements were determined by floating the feces samples in zinc sulfate (d = 1.33) solution on a McMaster slide and counting the eggs [30]. The detection limit for this technique was 15 eggs per gram (epg). The samples showed a low level of FEC, which was related to the exceptionally small amount of rainfall before and during the sampling period. For each flock, pooled feces were cultured to recover and identify third-stage larvae (L3) using standard parasitological techniques [30]. One hundred L3 were identified per flock to estimate the percentage of each helminth species.

IgA activity in serum was tested against a somatic antigen from the fourth-stage larvae (L4) of Teladorsagia circumcincta by indirect ELISA according to a modified protocol that was previously described by Martinez-Valladares et al. [31]. Briefly, ELISA plates (Sigma) were coated overnight with 100 µL of phosphate buffered saline (PBS) solution containing 2.5 µg/mL of T. circumcincta L4 somatic antigen. On the following day, the ELISA test was performed in four steps. After each step, the content of the plate was removed, the plate was washed, and each well was filled with a specific reagent; the plates were then incubated for 30 min. The following reagents were used for each step: (1) PT-Milk (4 g powdered milk + 100 mL PBS-Tween 20; PBS-Tween 20: 1 L PBS (pH 7.4) + 1 mL Tween 20 (Sigma)); (2) a sheep serum; (3) a rabbit anti-sheep IgA antibody and (4) a peroxidase substrate and tetramethylbenzidine solution to produce a color reaction that was stopped after 30 min by the addition of 50 μL of 2 M H2SO4. The results were measured as optical density (OD) values. Positive and negative controls were included in all the plates; positive controls were obtained from a pool of sera from sheep that were experimentally infected with T. circumcincta and negative controls were obtained from non-infected sheep that were maintained indoors. The results are expressed as optical density ratios (ODR) according to the following formula:

$$ODR = \frac{{\left( {sampleOD - negativeOD} \right)}}{{\left( {positiveOD - negativeOD} \right)}}$$

Statistical analyses

Prior to further analyses, FEC measurements were log-transformed (LFEC) to reduce over-dispersion, since no transformation yielding a normalized FEC dataset was available. However, Box-Cox power transformation was used for the IgA phenotype to obtain a normal distribution of values (IgAt). We used the R ‘car’ library to estimate the power parameter λ and carry out the transformation [32]; the log-transformation was also calculated through a command line in R [33].

To assess the variables that influence the two parasite resistance-related traits under study, an analysis of variance (ANOVA) was performed for LFEC and IgAt using a general linear model (GLM) through the R command line [33], which included the three following fixed effects: flock, age and time point relative to parturition. The ‘flock’ effect was classified into 17 groups. For the ‘age’ effect, two groups were considered i.e. ewes four to six years old and ewes seven or more years old. For the ‘time point relative to parturition’ effect, two categories were also considered i.e. one that included ewes that had a low immune response possibly because they were in the last stage of pregnancy or beginning lactation (animals sampled 2 weeks before giving birth or 30 days after birth) and one that included ewes that were outside that specific period.

Genotypes and physical map

We analyzed the genotypes that were obtained with the 50 K chip for a population of 1696 Churra ewes [34], which included animals with available phenotypic measurements for parasite resistance traits. First, SNP order and genome positions were updated according to the latest available version of the ovine Genome Assembly, Oar_v3.1 [35] by considering a 1 cM–1 Mb conversion rate. Then, quality control (QC) of the genotypes was performed for the entire genotyped population according to the protocol described in [34]. Briefly, QC was performed in seven steps that were applied to raw genotypes using the following criteria: (1) a GenCall score for raw genotypes greater than 0.15; (2) known location of the SNPs on the ovine autosomes; (3) a call rate per individual greater than 0.9; (4) a call rate per SNP greater or equal to 0.95; (5) minor allele frequency (MAF) higher than 0.05; (6) a p value for Hardy–Weinberg equilibrium (HWE) greater than 0.00001; and (7) analysis of the filtered genotypes using the VerifTyp software to check for Mendelian inconsistencies between parents and offspring (Boichard D and Druet T, personal communication). A total of 43,613 SNPs located on the 26 ovine autosomes passed the QC for the population of 1696 Churra ewes. For these 43,613 SNPs, available genotypes for 518 animals with parasite resistance phenotypes were subjected to different QTL mapping analyses.

QTL mapping analyses

Yield deviations (YD) of transformed data were used as dependent variables for statistical analyses to identify genomic regions that influence resistance to GIN infection. For the two traits under study, YD estimates were calculated following a multivariate animal model using the R command line and the ‘lsmeans’ library [36]. LFEC and IgAt were corrected for the fixed effect of ‘flock’, which according to the previously described ANOVA analysis, was the only factor that significantly influenced the studied traits. Then, the following statistical procedures were used for QTL mapping:

(1) Genome scans based on a classical LA and a combined LDLA procedure were performed at 0.1 cM step intervals using the corresponding analysis options (calcul = 4 and calcul = 28) of the QTLMap software [37]. Using this software, we also calculated the significance thresholds at the chromosome-wise significance level through a total of 1000 permutations (at 0.1 cM steps) for LA and 1000 simulations (at 5 cM steps) for LDLA. Genome-wise significance thresholds were based on the chromosome-wise significance threshold by correcting for the total number of chromosomes under analysis. A by-default haplotype size of four SNPs was used for LDLA.

For each QTL identified by the across-family LA scan, linkage-based within-family analyses were performed to identify the corresponding segregating families. For the significant QTL that were detected by LA, likelihood ratio test (LRT) values were converted to logarithm odds ratio (LOD) values [15], and confidence intervals (CI) for the QTL locations were estimated by the widely used 1-LOD drop-off method [38]. The proportion of phenotypic variance that was explained by the QTL detected by LA was calculated based on the corresponding LOD values using the formula \(\sigma_{p} = 1 - 10^{{\frac{ - 2}{n}LOD}}\) [39]. In the LDLA, chromosomal regions that involved consecutive significant haplotype associations within a chromosome (allowing gaps no greater than 5 cM) were grouped as a significant LDLA interval and the remaining ones were considered as isolated significant haplotypes.

For chromosomes with significant effects that were identified by both LA and LDLA genome scans, a linkage disequilibrium analysis (LDA) based on the LDA decay approach of Legarra and Fernando [40] was implemented using the QTLMap software (calcul = 26). The aim of this analysis was to determine whether the significant associations identified by LDLA were exclusively due to linkage pedigree-related information or whether an association with the trait could also be identified at the population level. Similar to the previously described LDLA, LDA was performed at 0.1 cM step intervals using a by-default 4-SNP haplotype size and 1000 (at 5 cM steps) simulations for the chromosome-wise threshold calculation. Significant LDA intervals were defined in the same way as for LDLA.

(2) A GWAS was performed by implementing the following linear mixed model (LMM), which includes the polygenic effect as a random effect and genotypes at single SNPs as fixed effects: (\({\mathbf{y}} = {\rm {\mathbf {Zu}}} + {\rm {\mathbf{X}}}b + e\)) where \({\mathbf{y}}\) is defined as the vector of phenotypes (YD) of the ewes; \({\mathbf{Z}}\) is a matrix associating random additive polygenic effects to individuals; u is a vector containing random polygenic effects; X is a vector with a genotypic indicator (−1, 0, or 1) that associates records to the marker effect; b is the allele substitution effect for the analyzed SNP; and e is the random residual. This association analysis was implemented by the restricted maximum likelihood (REML) method using the DMU package [41], and the SNP effect was tested using a Wald test against a null hypothesis of b = 0.

Bonferroni corrections for multiple-testing were used to estimate the genome-wise and chromosome-wise significant thresholds for the GWAS-based analyses. To account for the existence of linkage disequilibrium (LD) between the analyzed SNPs, rather than performing a conservative Bonferroni correction based on the total number of SNPs analyzed, we implemented the method proposed by Gao et al. [42] to calculate the number of independently analyzed SNPs for each chromosome and for the entire sheep genome. To this end, we used the simpleM test [43], which estimates the actual number of effective tests (Meff) in genome-wide association studies through a principal component analysis (PCA) approach. Using a PCA-cutoff of 0.975, the total number of independently analyzed SNPs across the entire genome was equal to 25,881.

Comparison with previously reported QTL and identification of functional candidate genes

We performed a systematic search for previously reported QTL and associations related to parasite resistance traits in sheep for which a good correspondence was observed with the significant associations that we identified in our study; in addition, we performed a search for positional candidate genes in relation to our results. However, prior to these searches, for each significant QTL and significant SNP association identified, we determined a “target genomic interval” (TGI), which was defined as the genomic region based on the sheep reference genome assembly Oar_v3.1 that corresponded to: (1) the CI that was estimated for the significant QTL detected by LA and for the defined significant LDLA intervals; and (2) a 250 kb-long interval centered on each of the significant isolated haplotypes detected by LDLA and the significant SNPs identified by GWAS.

Once the TGI were defined, they were compared with the Oar_v3.1 intervals that are annotated in the SheepQTL database (SheepQTLdb) [10] for previously reported QTL and that are mainly derived from microsatellite-based genome scans. We also compared these TGI with more recent data from studies based on the 50 K chip that are not included in this database [2126]. For some of these recent data based on the sheep genome assembly Oar_v2.0, when available, the corresponding Oar_v3.1 position of the target marker/interval was considered for comparison. Only regions that mapped within 1 Mb from the defined TGI were considered to coincide with our results. For the QTL that covered a very long region, the position of the QTL peak was prioritized to determine a possible correspondence.

The extraction of positional candidate genes included in the TGI according to the sheep genome assembly (Oar_v3.1) was performed using the BioMart web-based tool [44] based on the Ensembl release 81. Functional candidate genes related to the QTL identified in this study were identified by comparing the complete list of positional candidate genes extracted with BioMart with a database of 5029 immune-related genes. This database was based on the IRIS (1535 genes [45]) and ImmPort (4815 genes) gene lists, both of which are available at [46].



The presence of nematodes was confirmed in all the studied flocks with Trichostrongylus spp. and Teladorsagia spp. being the most prevalent species (49.3 and 48.6 %, respectively) that were identified among the total number of third-stage larvae obtained for the studied population. The prevalence of GIN infection by FEC per flock was 88.2 % (mean = 42.8 epg) and per individual was 45.4 % (mean = 39.4 epg). Faecal egg counts of GIN ranged from 0 to 1290 epg. For individual animals, the mean ODR of the IgA activity was 4.1 and ranged from 0.09 to 32.9.

QTL regions

The LA genome scan identified three 5 % chromosome-wise significant QTL (Table 1); in contrast, the LDLA genome scan identified 63 significant regions at the 5 % chromosome-wise level (Table 2). The LDA, which was performed for the three chromosomes that showed coincident results between the LA and LDLA scans, supported some of the significant signals that were identified previously (See Additional file 1: Table S1, Additional file 2: Figure S1). Although ten significant SNPs associated with IgAt (Table 3) were identified in the GWAS, no significant associations were detected for LFEC. The significant results are described below and those identified by more than one analysis are highlighted. For ease of comparison, Table 4 provides a summarized representation of the results of the three analyses performed across the entire genome (LA, LDLA and GWAS).

Table 1 Significant chromosome-wise QTL detected by linkage analysis (LA)
Table 2 Chromosome-wise significant results (Pc-value <0.01) from the combined linkage disequilibrium and linkage analysis (LDLA)
Table 3 Chromosome-wise SNPs significantly associated with the IgAt trait as identified by the GWAS
Table 4 Summary of the QTL detected by the three analyses performed in this study

LA results

The across-family regression analysis performed for LFEC and IgAt across the ovine autosomes identified three chromosome-wide significant QTL. Two of these QTL that are located on OAR6 (OAR for Ovis aries chromosome) (peak at 88.1 cM) and OAR8 (peak at 2 cM) had an effect on LFEC (Fig. 1a), whereas the other QTL located on OAR22 (peak at 3.4 cM) had effects on IgAt (Fig. 1b).

Fig. 1

Results of linkage analysis (LA; a, b) and combined linkage disequilibrium and linkage analysis (LDLA; c, d) genome scans performed for the two indicator traits of parasite resistance analyzed. Analyzed traits: LFEC Log-transformed faecal egg count, IgA t Box-Cox-transformed optical density ratio (ODR) values of immunoglobulin A activity. Likelihood ratio test (LRT) values obtained across the 26 ovine autosomes are represented. For those chromosomes that harbor significant QTL, the horizontal lines indicate the 5 % chromosome-wise significance threshold for LA (a, b) and the 5 % chromosome-wise significance threshold for LDLA (c, d)

The significant QTL identified by the across-family LA (maximum LRT value and CI estimated by the 1-LOD drop-off method), together with the results of the within-family analyses are in Table 1. The QTL for LFEC on OAR6 and OAR8 segregated in three and two families, respectively, whereas a single family was significant for the QTL for IgAt on OAR22. The CI that were estimated for the individual segregating families were located in the same region as the corresponding across-family CI, except for the peak for the QTL on OAR8 of Family 4, which was located at a more central position (31.2 cM) compared to the across-family peak at the proximal end of OAR8 (2 cM). However, the statistical profile for this family displayed a second peak reaching the 5 % chromosome-wise significance threshold (LRT = 11.76) at 12 cM, which was closer to the across-family QTL peak. The QTL effects estimated for the individual sires ranged from 0.3 (for the QTL for LFEC on OAR6) to 0.78 (for the QTL for LFEC on OAR8) standard deviations (Table 1). The estimated proportions of phenotypic variance explained by the three QTL identified by the LA were very similar and small (0.075, 0.077 and 0.069 % for the QTL on OAR6, 8 and 22, respectively).

LDLA results

Sixty-three significant QTL were detected at the 5 % chromosome-wise significance level by LDLA (30 for LFEC and 33 for IgAt). Among these 63 QTL, 13 (six for LFEC and seven for IgAt) reached the 5 % genome-wise significance level (Table 2; Fig. 1d). For 37 of the significant LDLA associations, nearby significant positions were grouped within a significant LDLA interval (Table 2); the remaining significant QTL identified by LDLA were defined based on isolated significant haplotypes. In addition, the three significant QTL identified by LA (on OAR6, 8 and 22) were supported by the LDLA scan (Table 2) (see Additional file 2: Figure S1). On OAR6, the LDLA results for LFEC revealed two 5 % chromosome-wise significant associations at 36 and 89.9 cM, with the latter being included within the CI of the QTL for LFEC on OAR6 detected by LA (Table 2). This analysis also identified a genome-wise significant association within the interval between 72.3 and 77.2 cM on OAR6.

On OAR8, although the LDLA scan identified a significant association at the proximal end of the chromosome (between 0.3 and 12.8 cM), which corresponded to the across-family CI for the QTL identified by LA, four other significant haplotype associations were identified across the chromosome (Table 2). Coincident with the QTL for IgAt on OAR22 detected by LA (between 0.3 and 5.8 cM), the LDLA scan revealed a chromosome-wise significant haplotype association (maximum LRT at 6.7 cM) at the proximal end of this chromosome.

LDA results

For the three chromosomes for which the QTLMap LDA approach was implemented, several 5 % chromosome-wise significant associations were identified for the same trait for which significant results were observed in the LA and LDLA (See Additional file 1: Table S1). A correspondence was found between the significant LDA association of the 75.8–85.1 cM region on OAR6 with LFEC and the LA and LDLA results. The other significant associations identified by LDA coincided with QTL detected by LDLA.

GWAS results

None of the analyzed SNPs reached significance for LFEC (Table 3; Fig. 2a). For IgAt, the GWAS identified one 5 % genome-wise significant SNP on OAR12 and nine additional 5 % chromosome-wise significant associations that were distributed on six chromosomes (OAR8, 10, 11, 14, 15 and 25) (Table 3; Fig. 2b). The allelic substitution effect of the significant SNPs identified for IgAt ranged from 0.243 to 0.417 phenotypic SD units. Although more than one significant SNP was identified on OAR8 and 10, these SNPs were located at relatively large distances on the chromosome (i.e., 22.8 and 13.9 Mb, respectively). Among the ten significant GWAS associations reported here for IgAt, one located on OAR10 was coincident with a significant QTL identified by LDLA for the same trait (between 21.5 and 27.2 cM), whereas two other associations, located on OAR8, overlapped with QTL for LFEC identified by LDLA.

Fig. 2

Results from the genome-wide association study (GWAS) performed for the two indicator traits of parasite resistance analyzed. Analyzed traits: LFEC Log-transformed faecal egg count, IgA t Box-Cox-transformed optical density ratio (ODR) values of immunoglobulin A activity. The values of the log(1/P-value) are shown for all the 43,613 SNPs that passed the quality control. For the chromosomes that harbor significant SNP associations, the horizontal lines indicate the 5 % chromosome-wise significance threshold obtained by applying a Bonferroni correction considering the number of independent SNPs analyzed for each chromosome. The genome-wise significance threshold, considering the number of independent markers analyzed for the entire genome is also represented

Correspondence of the detected associations with previously reported QTL for parasite resistance traits

The QTL for parasite resistance traits previously reported in sheep that coincide with the TGI reported here and are associated with the significant QTL and SNP associations identified here are summarized in Additional file 3: Table S2. Overall, we found correspondences with other studies for half of the 76 significant QTL identified by the three genome scans performed in this study.

List of functional candidate genes

A total of 905 unique genes were extracted from the TGI that were defined for the significant QTL detected by LA, LDLA and GWAS (416 and 489 unique genes extracted from FEC- and IgAt-associated regions, respectively) (see Additional file 4: Table S3). From the list of 5029 known immune-related genes, we performed a survey for positional candidate genes and identified 205 functional candidate genes (indicated in blue font in Additional file 4: Table S3), which were all extracted from TGI related to significant QTL that were detected by LA or LDLA. Gene symbols of these functional candidate genes are in Tables 1 and 2 based on their genomic locations within the corresponding QTL regions.


The genetic architecture of resistance to internal parasites is a complex trait that is influenced by many loci with small effects [21]. Using two different approaches to correct for sampling errors associated with single-marker regression, Kemper et al. [21] estimated that the largest effects that influence fecal worm egg count for Trichostrongylus colubriformis explained between 0.12 and 0.48 % of the phenotypic variance. These authors suggest that such small effects are shared by many complex traits and are not specific to parasite resistance. The proportions of phenotypic variance explained by the significant LA associations reported here, which were equal to ~0.074 %, are slightly lower than the lower limit of the range reported by Kemper et al. [21], although the estimated effects are within the ranges reported in other related studies [19, 20, 22]. Considering the small size of the targeted genetic effects to be detected, the statistical power of QTL detection for indicators of parasite resistance may be limited in such experiments if the number of sampled individuals is not very large. Based on Weller et al. [47], we estimated that the statistical power of QTL detection for QTL with a substitution effect of 0.2 phenotypic SD units, two alleles with frequencies of 0.25 and 0.75, respectively, and for a trait with a heritability of 0.2 (considering the estimates of Gutiérrez-Gil et al. [48]) was approximately 11 %. This estimate is based on the following assumptions i.e. (1) a type I error rate of 0.05, (2) a 1 % recombination frequency between the QTL and SNP and (3) 37.5 % of the analyzed sires are heterozygous at the QTL.

Our study successfully identified QTL that influence the two indicator traits related to GIN resistance using LA and LDLA, whereas the GWAS analysis only detected significant SNP associations with IgAt. The different analyses performed in this study can detect significant associations with different features. Hence, because classical LA will only detect QTL in our design if several sires are heterozygous at the same QTL (Qq), many marker-trait associations that do not satisfy this assumption but have a genuine association at the population level, will not be detected by LA; however, such associations can be detected by either of the two alternative genome scan analyses performed here i.e. LDLA or GWAS. Therefore, we attempted to present a global picture of the associations that segregate in this commercial sheep population by complementing the limits of classical LA with these alternative LDLA and GWAS approaches, which exploit population information. In our case, the GWAS approach also identified a substantially lower number of associations than LDLA. This may be explained by the fact that modeling both the association (LD) and the transmission (linkage) in a single analysis, LDLA permits to map QTL more accurately than LA while retaining its robustness to spurious associations [40]. In addition, among the different advantages highlighted for the use of LDLA versus GWAS for animal populations, Meuwissen et al. [49] claimed that LDLA is expected to suffer less from multiple-testing, and therefore to have more power to detect the existing QTL.

For the chromosomes that showed coincident significant results identified by LA and LDLA, we performed an exploratory LDA analysis with the QTLMap software (see Additional file 1: Table S1, Additional file 2: Figure S1). This analysis differs from GWAS in that parental haplotypes are pooled in classes that are defined by the identity-by-state (IBS) status of the haplotypes, with each different haplotype class having a specific effect on the quantitative trait [40]. The significant LDA results obtained for OAR6, 8 and 22 supported several of the significant LDLA associations reported for these chromosomes; whereas the LDA result obtained for OAR6 at 85.1 Mb supported the significant QTL that was detected by both LA and LDLA. This observation strengthens the support for the QTL for LFEC identified by LA on OAR6, which suggests that in addition to a family-based linkage information signal, the effect is also due to a genuine association with the trait, although it was not identified in our GWAS (most likely as a consequence of the limited power of the experimental design).

Regarding the LFEC-related results for OAR6 that were obtained by LA, LDLA and LDA, in the current study, we replicated the most significant QTL that was previously identified through a microsatellite-based genome scan using a different set of Churra sheep half-sib families [20]. In the latter study, the peak of the genome-wise significant QTL for LFEC was located in the marker interval BM4621-CSN3 on OAR6, which corresponds to a region between 68 and 85.1 Mb in the current sheep genome assembly (Oar_v3.1). The mentioned flanking interval overlaps with the TGI defined here for LFEC on OAR6 by LA (between 80.8 and 91.4 Mb) (Table 1), LDLA (between 72.3 and 77.2 and between 85 and 90.2 Mb) (Table 2) and LDA (between 75.8 and 77.7 and between 85 and 85.1 Mb) (see Additional file 1: Table S1, Additional file 2: Figure S1). This finding provides support for the design and planning of future fine-mapping studies for this chromosomal region. The higher marker density and information provided by the complementary analyses reported here for this region suggest that the OAR6 region ranging from 68 to 91.4 Mb includes several different QTL that directly influence GIN resistance in Churra sheep. Interestingly, a GWAS on a Red Maasai x Dorper backcross sheep population [26] also suggested the presence of several QTL for FEC in lambs within a region between 55.9 and 78.19 Mb on OAR6. This finding was based on the fact that the most significant SNP association with FEC identified on OAR6 at 74.86 Mb, was proven not to be in LD with nearby clusters of significant markers for the same trait (in intervals between 55.9 and 62.6 Mb, 74.1 and 75.00 Mb, and 78.1 and 78.2 Mb) (see Additional file 3: Table S2). In spite of the remarkable correspondence between these results and our results, the most distal signals that were detected on OAR6 in our study (TGI defined by LA: 80.8 to 91.4 Mb; LDLA: 85 to 90.2 Mb; and LDA: 85 to 85.1 Mb) do not overlap with any previously reported QTL in other populations, but only with those previously reported by Gutiérrez-Gil et al. [20] (see Additional file 3: Table S2) . With the exception of Gutiérrez-Gil et al. [20] work, most studies refer to QTL that are detected for young animals (lambs); thus, the most distal QTL that we identified on OAR6 could be related to specific mechanisms of the immune response that is activated in adult animals. As suggested by Stear et al. [50], the genetic variation in fecal egg counts in lambs is a consequence of genetic variation in worm length and hence worm fecundity; in contrast, mature sheep may be able to regulate both fecundity and worm number. These authors suggested that the lower fecal egg counts observed in adult animals compared to lambs are due to the acquisition of effective immune responses that reduce worm numbers, possibly via immediate hypersensitivity reactions against incoming third-stage larvae [51]. Recent studies have highlighted differences in the pathways involved in innate and acquired resistance [52]. Another correspondence that was observed with the results reported by Gutiérrez-Gil et al. [20] concerned the QTL for LFEC detected by LDLA on OAR10 (TGI: 70.01–71.55 Mb) (see Additional file 3: Table S2). Due to the lack of evidence from the other analyses reported here, this region was not further investigated.

An intriguing finding is that the other two QTL detected by LA in this work did not coincide with QTL that were reported for other sheep populations, whereas three of the ten significant SNP associations identified by GWAS, and 35 of the 63 significant QTL identified by LDLA, overlapped with QTL effects described in other studies (see Additional file 3: Table S2). Indeed, the significant GWAS results coincided with QTL on OAR8 reported by Crawford et al. [13] and Silva et al. [19], on OAR12 by Riggio et al. [24], and on OAR15 by Silva et al. [19] (see Additional file 3: Table S2). In our study, the SNP association on OAR12 at 61.9 Mb was the only one that reached the 5 % genome-wide significance level. Although not mentioned in Additional file 3: Table S2 because there was no complete overlap, Beh et al. [12] used microsatellite markers to identify a QTL in this genomic region (between 63.5 and 71.5 Mb) for FEC-related traits in T. colubriformis infection. It should be noted that we did not find a clear correspondence with the classical regions reported to influence parasite resistance traits, such as those that harbor the ovine IFN-γ gene (OAR3: 151.53 Mb) [11, 14, 17] or the major histocompatibility complex-related genes (OAR20: 7 Mb; 24–26 Mb; 58–60 Mb) [14].

Among the large number of correspondences between our LDLA results and previously reported studies (see Additional file 3: Table S2), those that are based on data from the 50 K chip are of special relevance because of the proximity between the QTL peaks reported here and in other studies. Apart from the correspondences with the findings of Benavides et al. [26] mentioned above for OAR6, those found for the QTL on OAR5 (TGI: 89.68–90.14 Mb) are particularly relevant. This QTL identified by LDLA is located in a region where several significant effects for a wide range of parasite indicator traits were reported by Sallé et al. [22], which suggests the presence of a QTL with pleiotropic effects.

We identified 205 immune-related genes within the TGI defined by the LA and LDLA (Tables 1, 2) but none of these functional candidate genes were found in the significant GWAS-defined TGI. Some of these immune-related genes are involved in the T helper (Th) 2 cell response, which orchestrates the mechanisms of tissue repair as a primary host defense against helminthes [53], whereas others are linked to the Th1 cell response, which is associated with progression to chronic infection [54].

Due to the large number of significant regions identified and the need for additional fine-mapping results to propose reliable promising causal candidate genes, in the following part, we only discuss below the genes that were identified in relation to the QTL for LFEC identified by LA on OAR6 (TGI: 80.9–91.4 Mb), which include the genes extracted for the LDLA-defined TGI between 85 and 90.2 Mb. The fact that this QTL, previously reported by Gutiérrez-Gil et al. [20], was also identified for the population analyzed here and the support provided by the related signals identified by LDLA/LDA, led us to carry out a preliminary assessment of the 20 positional candidate immune-related genes that map to this region (Table 1). Among these genes, several encode chemokines (IL8, CXCL1, CXCL10, CXCL11, CXCL9, PF4, PPBP), a family of small proteins that play important roles in the immune system through leukocyte recruitment, cell communication and cell activation during infection [55, 56]. In particular, IL8 (or CXCL8) and CXCL1 are involved in the recruitment and activation of neutrophils [55]. IL8 also participates in the recruitment of mast cells, which are frequently associated with the Th2 cell response [57]. CXCL9, CXCL10 and CXCL11, which are induced by IFN-γ, are involved in promoting the Th1 immune response. In nematode-infected mice, CXCL10 slows down the intestinal epithelial cell turnover rate and thus, increases worm survival [58]. In addition, both PF4 and PPBP have been suggested to play roles in wound healing [59, 60]. Three genes coding for members of the epidermal growth factor family also map to the considered region on OAR6: AREG (amphiregulin), BTC (betacellulin) and EREG (epiregulin). AREG is expressed by diverse cell types involved in the immune response, such as activated Th2 cells [61], and is a central mediator of epithelial repair [62]. In mice, lack of AREG expression appears to have an effect on the delayed expulsion of GIN [63]. Because wound repair and GIN expulsion are related to the acquired Th2 response [53, 64], genes associated with these mechanisms (e.g., IL8, PF4, PPBP and AREG) could be of interest when searching for candidates to explain an adult-specific QTL, such as the QTL detected on OAR6 between 80.8 and 91.4 Mb.

The large number of QTL identified in this study supports the idea that disease susceptibility is not determined by individual genes acting alone but rather by complex multi-gene interactions [65, 66]. Our results are the first steps towards the identification of allelic variants that directly control the phenotypic variation observed for parasite resistance in adult Churra sheep. The identification of causal variants, or SNPs in strong LD with the casual variants, could contribute to the implementation of these results in breeding schemes for the Churra breed population. Future studies that combine genomic variation analysis and functional genomic information may help to elucidate the biology of resistance to GIN diseases in sheep.


In summary, the 50 K chip was used for a medium marker density scan of the sheep genome to identify regions that influence traits related to resistance to GIN infections in adult animals. By exploiting the information obtained at the within-family level and at the population level, three methods of analysis were used (LA, LDLA and GWAS) to provide a global picture of the QTL that segregate in the commercial population of Churra sheep analyzed. Many of the significant associations reported here overlap with previously reported QTL for different populations of young sheep. These results will contribute to identify target regions that control variation of the complex parasite resistance trait in sheep, independently of the age of the animals. Other significant associations that did not coincide with previously reported QTL could be related to the specific immune response of adult animals. This study also replicated a QTL for FEC on OAR6 that was previously reported in a different subset of animals from the commercial population of Churra sheep. Together, the enhanced marker density provided by the 50 K chip and the complementary analyses reported here suggest that several QTL are present in this genomic region. This replication and the re-definition of these genetic effects in the independent population analyzed here provide support for investing future research efforts aimed at identifying the corresponding causal allelic variants. The combination of high-density SNP genotyping (700 K SNP array) and whole-genome sequencing of segregating trios (composed by a segregating sire carrying the Qq genotype, and two homozygous daughters for alternative haplotype alleles, QQ and qq, and showing extreme divergence for the resistance phenotype) could be a powerful strategy to reach this objective.


  1. 1.

    Kaplan RM, Vidyashankar AN. An inconvenient truth: global worming and anthelmintic resistance. Vet Parasitol. 2012;186:70–8.

    PubMed  Article  Google Scholar 

  2. 2.

    Raadsma HW, Gray GD, Woolaston RR. Genetics of disease resistance and vaccine response. In: Piper L, Ruvinsky A, editors. The Genetics of Sheep. University Press: Cambridge; 1997. p. 199–224.

  3. 3.

    Stear MJ, Doligalska M, Donskow-Schmelter K. Alternatives to anthelmintics for the control of nematodes in livestock. Parasitology. 2007;134:139–51.

    PubMed  CAS  Article  Google Scholar 

  4. 4.

    Morris CA, Wheeler M, Watson TG, Hosking BC, Leathwick DM. Direct and correlated responses to selection for high or low faecal nematode egg count in Perendale sheep. NZJ Agric Res. 2005;48:1–10.

    Article  Google Scholar 

  5. 5.

    Karlsson LJE, Greeff JC. Selection response in fecal worm egg counts in the Rylington Merino parasite resistant flock. Aust J Exp Agric. 2006;46:809–11.

    Article  Google Scholar 

  6. 6.

    Kemper KE, Elwin RL, Bishop SC, Goddard ME, Woolaston RR. Haemonchus contortus and Trichostrongylus colubriformis did not adapt to long-term exposure to sheep that were genetically resistant or susceptible to nematode infections. Int J Parasitol. 2009;39:607–14.

    PubMed  CAS  Article  Google Scholar 

  7. 7.

    Grencis RK. Immunity to helminths: resistance, regulation, and susceptibility to gastrointestinal nematodes. Annu Rev Immunol. 2015;33:201–25.

    PubMed  CAS  Article  Google Scholar 

  8. 8.

    Jiang Y, Xie M, Chen W, Talbot R, Maddox JF, Faraut T, et al. The sheep genome illuminates biology of the rumen and lipid metabolism. Science. 2014;344:1168–73.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  9. 9.

    Hein WR, Pernthaner A, Piedrafita D, Meeusen EN. Immune mechanisms of resistance to gastrointestinal nematode infections in sheep. Parasite Immunol. 2010;32:541–8.

    PubMed  CAS  Google Scholar 

  10. 10.

    Hu ZL, Park CA, Wu XL, Reecy JM. Animal QTLdb: an improved database tool for livestock animal QTL/association data dissemination in the post-genome era. Nucleic Acids Res. 2013;41:D871–9.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  11. 11.

    Coltman DW, Wilson K, Pilkington JG, Stear MJ, Pemberton JM. A microsatellite polymorphism in the gamma interferon gene is associated with resistance to gastrointestinal nematodes in a naturally-parasitized population of Soay sheep. Parasitology. 2001;122:571–82.

    PubMed  CAS  Article  Google Scholar 

  12. 12.

    Beh KJ, Hulme DJ, Callaghan MJ, Leish Z, Lenane I, Windon RG, et al. A genome scan for quantitative trait loci affecting resistance to Trichostrongylus colubriformis in sheep. Anim Genet. 2002;33:97–106.

    PubMed  CAS  Article  Google Scholar 

  13. 13.

    Crawford AM, Paterson KA, Dodds KG, Diez Tascon C, Williamson PA, Roberts Thomson M, et al. Discovery of quantitative trait loci for resistance to parasitic nematode infection in sheep: I. Analysis of outcross pedigrees. BMC Genomics. 2006;7:178.

    PubMed  PubMed Central  Article  Google Scholar 

  14. 14.

    Davies G, Stear MJ, Benothman M, Abuagob O, Kerr A, Mitchell S, et al. Quantitative trait loci associated with parasitic infection in Scottish blackface sheep. Heredity (Edinb). 2006;96:252–8.

    CAS  Article  Google Scholar 

  15. 15.

    Beraldi D, McRae AF, Gratten J, Slate J, Visscher PM, Pemberton JM. Mapping quantitative trait loci underlying fitness-related traits in a free-living sheep population. Evolution. 2007;61:1403–16.

    PubMed  Article  Google Scholar 

  16. 16.

    Marshall K, Maddox JF, Lee SH, Zhang Y, Kahn L, Graser HU, et al. Genetic mapping of quantitative trait loci for resistance to Haemonchus contortus in sheep. Anim Genet. 2009;40:262–72.

    PubMed  CAS  Article  Google Scholar 

  17. 17.

    Dominik S, Hunt PW, McNally J, Murrell A, Hall A, Purvis IW. Detection of quantitative trait loci for internal parasite resistance in sheep. I. Linkage analysis in a Romney × Merino sheep backcross population. Parasitology. 2010;137:1275–82.

    PubMed  CAS  Article  Google Scholar 

  18. 18.

    Matika O, Pong-Wong R, Woolliams JA, Bishop SC. Confirmation of two quantitative trait loci regions for nematode resistance in commercial British terminal sire breeds. Animal. 2011;5:1149–56.

    PubMed  CAS  Article  Google Scholar 

  19. 19.

    Silva MV, Sonstegard TS, Hanotte O, Mugambi JM, Garcia JF, Nagda S, et al. Identification of quantitative trait loci affecting resistance to gastrointestinal parasites in a double backcross population of Red Maasai and Dorper sheep. Anim Genet. 2012;43:63–71.

    PubMed  CAS  Article  Google Scholar 

  20. 20.

    Gutierrez-Gil B, Perez J, Alvarez L, Martinez-Valladares M, de la Fuente LF, Bayon Y, et al. Quantitative trait loci for resistance to trichostrongylid infection in Spanish Churra sheep. Genet Sel Evol. 2009;41:46.

    PubMed  PubMed Central  Article  Google Scholar 

  21. 21.

    Kemper KE, Emery DL, Bishop SC, Oddy H, Hayes BJ, Dominik S, et al. The distribution of SNP marker effects for faecal worm egg count in sheep, and the feasibility of using these markers to predict genetic merit for resistance to worm infections. Genet Res (Camb). 2011;93:203–19.

    PubMed  CAS  Article  Google Scholar 

  22. 22.

    Sallé G, Jacquiet P, Gruner L, Cortet J, Sauvé C, Prévot F, et al. A genome scan for QTL affecting resistance to Haemonchus contortus in sheep. J Anim Sci. 2012;90:4690–705.

    PubMed  Article  Google Scholar 

  23. 23.

    Riggio V, Matika O, Pong-Wong R, Stear MJ, Bishop SC. Genome-wide association and regional heritability mapping to identify loci underlying variation in nematode resistance and body weight in Scottish Blackface lambs. Heredity (Edinb). 2013;110:420–9.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  24. 24.

    Riggio V, Pong-Wong R, Sallé G, Usai MG, Casu S, Moreno CR, et al. A joint analysis to identify loci underlying variation in nematode resistance in three European sheep populations. J Anim Breed Genet. 2014;131:426–36.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  25. 25.

    McRae KM, McEwan JC, Dodds KG, Gemmell NJ. Signatures of selection in sheep bred for resistance or susceptibility to gastrointestinal nematodes. BMC Genomics 2014;15:637.

    PubMed  PubMed Central  Article  Google Scholar 

  26. 26.

    Benavides MV, Sonstegard TS, Kemp S, Mugambi JM, Gibson JP, Baker RL, et al. Identification of novel loci associated with gastrointestinal parasite resistance in a Red Maasai × Dorper backcross population. PLoS One. 2015;10:e0122797.

    PubMed  PubMed Central  Article  Google Scholar 

  27. 27.

    Garcı́a-Pérez AL, Hurtado A, Oregui LM, Juste RA. Effects of a second annual strategic anthelmintic treatment in dairy sheep in Northern Spain. Small Rumin Res. 2002;43:121–6.

  28. 28.

    Houdijk JGM, Kyriazakis I, Coop RL, Jackson F. The expression of immunity to Teladorsagia circumcincta in ewes and its relationship to protein nutrition depend on body protein reserves. Parasitology. 2001;122:661–72.

    PubMed  CAS  Article  Google Scholar 

  29. 29.

    Wolstenholme AJ, Fairweather I, Prichard R, von Samson-Himmelstjerna G, Sangster NC. Drug resistance in veterinary helminths. Trends Parasitol. 2004;20:469–76.

    PubMed  CAS  Article  Google Scholar 

  30. 30.

    Fish and Food Ministry of Agriculture. Manual of veterinary parasitological laboratory techniques. 1986.

  31. 31.

    Martinez-Valladares M. Vara-Del Rio MP, Cruz-Rojo MA, Rojo-Vazquez FA. Effect of a low protein diet on the resistance of Churra sheep to Teladorsagia circumcincta. Parasite Immunol. 2005;27:219–25.

    PubMed  CAS  Article  Google Scholar 

  32. 32.

    Fox J, Weisberg S, Bates D, Fox MJ. Package ‘car’. Vienna: R Foundation for Statistical Computing; 2012.

    Google Scholar 

  33. 33.

    R Core Team. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna. 2014.

  34. 34.

    García-Gámez E, Gutiérrez-Gil B, Sahana G, Sánchez JP, Bayón Y, Arranz JJ. GWA analysis for milk production traits in dairy sheep and genetic support for a QTN influencing milk protein percentage in the LALBA gene. PLoS One. 2012;7:e47782.

    PubMed  PubMed Central  Article  Google Scholar 

  35. 35.

    The CSIRO Animal, Food and Health Sciences livestock genomics web site. CSIRO Australia. The Ovis aries reference genome assembly ( Accessed 21 Jan 2013.

  36. 36.

    Lenth Russell V. lsmeans: least-squares means. R package version 1.06–05. 2013.

  37. 37.

    Filangi O, Moreno C, Gilbert H, Legarra A, Le Roy P, Elsen J. QTLMap, a software for QTL detection in outbred populations. In Proceedings of the 9th World Congress on Genetics Applied to Livestock Production: 1–6 August 2010; Leipzig; 2010.

  38. 38.

    Lander ES, Botstein D. Mapping mendelian factors underlying quantitative traits using RFLP linkage maps. Genetics. 1989;121:185–99.

    PubMed  CAS  PubMed Central  Google Scholar 

  39. 39.

    Broman KW, Sen S. A Guide to QTL Mapping with R/qtl (Vol. 46). New York: Springer; 2009.

  40. 40.

    Legarra A, Fernando RL. Linear models for joint association and linkage QTL mapping. Genet Sel Evol. 2009;41:43.

    PubMed  PubMed Central  Article  Google Scholar 

  41. 41.

    Madsen P, Sørensen P, Su G, Damgaard LH, Thomsen H, Labouriau R. DMU-a package for analyzing multivariate mixed models. In Proceedings of the 8th World Congress on Genetics Applied to Livestock Production: 13–18 August 2006; Belo Horizonte. 2006.

  42. 42.

    Gao X, Becker LC, Becker DM, Starmer JD, Province MA. Avoiding the high Bonferroni penalty in genome-wide association studies. Genet Epidemiol. 2010;34:100–5.

    PubMed  PubMed Central  Google Scholar 

  43. 43.

    simpleM. A multiple testing correction program for correlated SNPs. ( Accessed 15 Nov 2014.

  44. 44.

    Ensembl release 81—August 2015 WTSI/EMBL-EBI ( Accessed 15 Aug 2015.

  45. 45.

    Kelley J, de Bono B, Trowsdale J. IRIS: a database surveying known human immune system genes. Genomics. 2005;85:503–11.

    PubMed  CAS  Article  Google Scholar 

  46. 46.

    Innate Database Gene Lists: ( Accessed 15 Feb 2015.

  47. 47.

    Weller JI, Kashi Y, Soller M. Power of daughter and granddaughter designs for determining linkage between marker loci and quantitative trait loci in dairy cattle. J Dairy Sci. 1990;73:2525–37.

    PubMed  CAS  Article  Google Scholar 

  48. 48.

    Gutiérrez-Gil B, Pérez J, de la Fuente L, Meana A, Martínez-Valladares M, San Primitivo F, et al. Genetic parameters for resistance to trichostrongylid infection in dairy sheep. Animal. 2010;4:505–12.

    PubMed  Article  Google Scholar 

  49. 49.

    Meuwissen, T. Use of whole genome sequence data for QTL mapping and genomic selection. In Proceedings of the 9th World Congress on Genetics Applied to Livestock Production: 1–6 August 2010; Leipzig; 2010.

  50. 50.

    Stear MJ, Strain S, Bishop SC. Mechanisms underlying resistance to nematode infection. Int J Parasitol. 1999;29:51–6.

    PubMed  CAS  Article  Google Scholar 

  51. 51.

    Stear MJ, Bishop SC, Doligalska M, Duncan JL, Holmes PH, Irvine J, et al. Regulation of egg production, worm burden, worm length and worm fecundity by host responses in sheep infected with Ostertagia circumcincta. Parasite Immunol. 1995;17:643–52.

    PubMed  CAS  Article  Google Scholar 

  52. 52.

    Grencis RK, Humphreys NE, Bancroft AJ. Immunity to gastrointestinal nematodes: mechanisms and myths. Immunol Rev. 2014;260:183–205.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  53. 53.

    Allen JE, Wynn TA. Evolution of Th2 immunity: a rapid repair response to tissue destructive pathogens. PLoS Pathog. 2011;7:e1002003.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  54. 54.

    Else KJ, Finkelman FD, Maliszewski CR, Grencis RK. Cytokine-mediated regulation of chronic intestinal helminth infection. J Exp Med. 1994;179:347–51.

    PubMed  CAS  Article  Google Scholar 

  55. 55.

    Schumacher C, Clark-Lewis I, Baggiolini M, Moser B. High- and low-affinity binding of GRO alpha and neutrophil-activating peptide 2 to interleukin 8 receptors on human neutrophils. Proc Natl Acad Sci USA. 1992;89:10542–6.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  56. 56.

    Trotta T, Costantini S, Colonna G. Modelling of the membrane receptor CXCR3 and its complexes with CXCL9, CXCL10 and CXCL11 chemokines: putative target for new drug design. Mol Immunol. 2009;47:332–9.

    PubMed  CAS  Article  Google Scholar 

  57. 57.

    da Silva EZ, Jamur MC, Oliver C. Mast cell function: a new vision of an old cell. J Histochem Cytochem. 2014;62:698–738.

    PubMed  PubMed Central  Article  Google Scholar 

  58. 58.

    Cliffe LJ, Humphreys NE, Lane TE, Potten CS, Booth C, Grencis RK. Accelerated intestinal epithelial cell turnover: a new mechanism of parasite expulsion. Science. 2005;308:1463–5.

    PubMed  CAS  Article  Google Scholar 

  59. 59.

    Senior RM, Griffin GL, Huang JS, Walz DA, Deuel TF. Chemotactic activity of platelet alpha granule proteins for fibroblasts. J Cell Biol. 1983;96:382–5.

    PubMed  CAS  Article  Google Scholar 

  60. 60.

    Dvonch VM, Murphey RJ, Matsuoka J, Grotendorst GR. Changes in growth factor levels in human wound fluid. Surgery. 1992;112:18–23.

    PubMed  CAS  Google Scholar 

  61. 61.

    Zaiss DM, van Loosdregt J, Gorlani A, Bekker CP, Gröne A, Sibilia M, et al. Amphiregulin enhances regulatory T cell-suppressive function via the epidermal growth factor receptor. Immunity. 2013;38:275–84.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  62. 62.

    Monticelli L, Sonnenberg G, Abt M, Osborne L, Wojno ET, Alenghat T, et al. Innate lymphoid cells promote airway epithelial repair through the amphiregulin-EGFR pathway (P3250). J Immunol. 2013;190:136–7.

    Google Scholar 

  63. 63.

    Zaiss DM, Yang L, Shah PR, Kobie JJ, Urban JF, Mosmann TR. Amphiregulin, a TH2 cytokine enhancing resistance to nematodes. Science. 2006;314:1746.

    PubMed  CAS  Article  Google Scholar 

  64. 64.

    Allen JE, Sutherland TE. Host protective roles of type 2 immunity: parasite killing and tissue repair, flip sides of the same coin. Semin Immunol. 2014;26:329–40.

    PubMed  CAS  PubMed Central  Article  Google Scholar 

  65. 65.

    Kadarmideen HN, Watson-Haigh NS, Andronicos NM. Systems biology of ovine intestinal parasite resistance: disease gene modules and biomarkers. Mol BioSyst. 2011;7:235–46.

    PubMed  CAS  Article  Google Scholar 

  66. 66.

    Moore JH. The ubiquitous nature of epistasis in determining susceptibility to common human diseases. Hum Hered. 2003;56:73–82.

    PubMed  Article  Google Scholar 

Download references

Authors’ contributions

JJA and BGG conceived and designed the study and analyses; MA and MMV carried out the data collection and prepared the phenotype; MA, BGG, and JJA analyzed the data; MA, BGG and JJA drafted the manuscript. All authors read and approved the final manuscript.


The authors would like to acknowledge JM Elsen (INRA, France) for providing assistance with interpretation of the QTLMap output. This work was supported by a competitive grant from the Castilla and León regional government (Junta de Castilla y León) (Ref. LE245A12-2) and the AGL2012-34437 project funded by the Spanish Ministry of Economy and Competitiveness (MINECO). M Atlija is a grateful grantee of a Marie Curie fellowship funded by the EC-funded Initial Training Network (ITN) NematodeSystemHealth (FP7-PEOPLE-2010-ITN Ref. 264639). B Gutiérrez-Gil is funded through the Spanish “Ramón y Cajal” Programme (RYC-2012-10230) from the MINECO.

Competing interests

The authors declare that they have no competing interests.

Author information



Corresponding author

Correspondence to Beatriz Gutiérrez-Gil.

Additional files

Additional file 1: Table S1.

Chromosome-wise significant results (Pc-value < 0.05) identified by the linkage disequilibrium analysis (LDA) performed in the present study for chromosomes (OAR) 6, 8 and 22. Characterization of the chromosome-wise significant results (Pc-value < 0.05) identified by the QTLMap linkage disequilibrium analysis (LDA) that was performed for the three chromosomes showing coincident results in the LA and LDLA genome scans presented here for parasite resistance traits.

Additional file 2: Figure S1.

Profiles of the Likelihood Ratio Test (LRT) obtained from the linkage analysis (LA), linkage disequilibrium analysis (LDA) and the combined LDLA performed for chromosomes (OAR) 6 (a; LFEC), 8 (b; LFEC), and 22 (c; IgAt). For the indicated trait, the LRT results of LA (solid line), LDA (dark gray circle), and LDLA (light gray circle) (y-axis) are plotted against the SNP positions analyzed along chromosomes (OAR) 6, 8 and 22 (x-axis). The 5 % chromosome-wise significance thresholds considered for each of three analyses are represented as horizontal lines.

Additional file 3: Table S2.

Summary table of the correspondence between the QTL and SNP associations identified in the present study and other studies previously reported for parasite resistance traits in sheep. This table shows the correspondences found for all the QTL identified in this study (by LA, LDLA and GWAS) (indicated in green cells) with QTL previously reported based on microsatellite-based studies (compiled in the SheepQTLdb; indicated in light orange cells) and SNP chip-based studies (indicated in orange cells).

Additional file 4: Table S3.

Total list of annotated genes extracted from the Sheep Genome Assembly v3.1 using the BioMart web-tool for the significant QTL regions and SNP associations identified for the two parasite resistance traits analyzed in the present study. Among the total list of genes extracted, we identified 205 functional candidate genes involved in the immune response, based on our candidate gene survey, which are indicated in blue font colour. The colour of the rows refer to genes extracted based on the results of the Linkage Analysis (LA; green), Combined Linkage Disequilibrium and Linkage Analysis (LDLA; yellow) and Genome-wise Association Study (GWAS; blue).

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (, 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 ( applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Atlija, M., Arranz, JJ., Martinez-Valladares, M. et al. Detection and replication of QTL underlying resistance to gastrointestinal nematodes in adult sheep using the ovine 50K SNP array. Genet Sel Evol 48, 4 (2016).

Download citation


  • Quantitative Trait Locus
  • Linkage Analysis
  • Significant Quantitative Trait Locus
  • Linkage Disequilibrium Analysis
  • Quantitative Trait Locus Peak