Skip to main content

Genomic regions and pathways associated with gastrointestinal parasites resistance in Santa Inês breed adapted to tropical climate

Abstract

Background

The aim of this study was to estimate variance components and to identify genomic regions and pathways associated with resistance to gastrointestinal parasites, particularly Haemonchus contortus, in a breed of sheep adapted to tropical climate. Phenotypes evaluations were performed to verify resistance to gastrointestinal parasites, and were divided into two categories: i) farm phenotypes, assessing body condition score (BCS), degree of anemia assessed by the famacha chart (FAM), fur score (FS) and feces consistency (FC); and ii) lab phenotypes, comprising blood analyses for hematocrit (HCT), white blood cell count (WBC), red blood cell count (RBC), hemoglobin (HGB), platelets (PLT) and transformed (log10) egg per gram of feces (EPGlog). A total of 576 animals were genotyped with the Ovine SNP12k BeadChip (Illumina, Inc.), that contains 12,785 bialleleic SNP markers. The variance components were estimated using a single trait model by single step genomic BLUP procedure.

Results

The overall linkage disequilibrium (LD) mean between pairs of markers measured by r 2 was 0.23. The overall LD mean between markers considering windows up to 10 Mb was 0.07. The mean LD between adjacent SNPs across autosomes ranged from 0.02 to 0.10. Heritability estimates were low for EPGlog (0.11), moderate for RBC (0.18), PLT (0.17) HCT (0.20), HGB (0.16) and WBC (0.22), and high for FAM (0.35). A total of 22, 21, 23, 20, 26, 25 and 23 windows for EPGlog for FAM, WBC, RBC, PLT, HCT and HGB traits were identified, respectively. Among the associated windows, 10 were shown to be common to HCT and HGB traits on OAR1, OAR2, OAR3, OAR5, OAR8 and OAR15.

Conclusion

The traits indicating gastrointestinal parasites resistance presented an adequate genetic variability to respond to selection in Santa Inês breed, and it is expected a higher genetic gain for FAM trait when compared to the others. The level of LD estimated for markers separated by less than 1 Mb indicated that the Ovine SNP12k BeadChip might be a suitable tool for identifying genomic regions associated with traits related to gastrointestinal parasite resistance. Several candidate genes related to immune system development and activation, inflammatory response, regulation of lymphocytes and leukocytes proliferation were found. These genes may help in the selection of animals with higher resistance to parasites.

Background

Small ruminants, like sheep (Ovis aries) and goats (Capra aegagrus hircus) were the first animals domesticated by human for food supply, being the most important group of ruminants raised in temperate and tropical regions [1] . Sheep are multi-purpose animals, raised for meat, milk, wool, hides and skins, having a huge socioeconomic importance worldwide. However, there are several productive drawbacks associated with gastrointestinal parasites infection in small ruminants, since it represents the type of disease with the highest impact on animal health and productivity [2]. Thus, the losses caused by gastrointestinal parasites and the costs due to excessive use of anthelmintic drugs are a struggling problem that restricts the sheep production in many regions of the world.

The main gastrointestinal nematodes of small ruminants belong to the Trichostrongylidae family. These parasites occur in the gastrointestinal tract and are seen as the main obstacle in the sheep industry, since they lead to significant economic losses due to high mortality and productivity losses. Within the Trichostrongylidae family, the Haemonchus genera has great pathogenic action and is the most prevalent parasite affecting small ruminants, mainly in tropical regions, where environmental conditions are characterized by high temperature and humidity, and abundant rainfall during summer [3,4,5,6].

To reduce the economic losses caused by nematode infections there are several management alternatives to minimize the damage, such as raise breeds or animals that are more resistant to these infections [5]. In this regard, Santa Inês, an American hair breed, showed higher resistance to gastrointestinal nematode infections when compared to European sheep breeds [5, 7,8,9]. Several studies have shown that selection for gastrointestinal parasite resistance is possible in sheep, since genetic progress in research and commercial herds [10,11,12,13,14,15,16].

The identification of genetic markers associated with gastrointestinal parasites resistance could increase the genetic response for this trait by marker-assisted selection [15]. Furthermore, the identification of genomic regions associated with resistance or susceptibility to gastrointestinal parasites will help to deeper understand the biological and physiological processes underlying this trait [17]. Several studies reported genetic markers associated with gastrointestinal parasites resistance close to the Major Histocompatibility Complex (MHC) [18,19,20] and IFN-gama genes [18,19,20,21,22]. Recently, genome-wide association studies have identified genetic variants for gastrointestinal parasite resistance in some sheep breeds [23,24,25,26]. The identification of genomic regions that play a role in gastrointestinal parasite resistance may become an important tool to improve the resistance of Santa Inês or other sheep breeds. Therefore, the aim of this study was to estimate variance components and to identify genomic regions and pathways associated with resistance to gastrointestinal parasites in a Santa Inês population adapted to tropical conditions.

Methods

Data

The phenotypic records were collected from 700 naturally infected animals of Santa Inês breed belonging to four flocks located in the Minas Gerais and São Paulo southeast states and Sergipe northeast state of Brazil. The samples collection and phenotypes evaluations were performed from October to November of 2011. Phenotypes evaluations were achieved to verify resistance to gastrointestinal parasites, and were divided into two categories: i) farm phenotypes, assessing body condition score (BCS), degree of anemia assessed by the famacha card (FAM), fur score (FS) and feces consistency (FC); and ii) lab phenotypes, comprising blood analyses for hematocrit (HCT), white blood cell count (WBC), red blood cell count (RBC), hemoglobin (HGB), platelets (PLT) and the egg counts per gram of feces (EPGlog). The EPG values were transformed to log10 (n + 1) = EPGlog to meet the basic requirements of normality and homogeneity in an attempt to stabilize the variance prior to analyses, where n is the number of EPGlog per sample.

Blood samples were collected by puncture of the jugular vein using vacuum tubes with anticoagulant EDTA and clot activator. Subsequently, samples were submitted to the Veterinary hematology analyzer Sysmex PocH-100iV Diff to perform a complete blood cell count for HCT, WBC, RBC, HGB and PLT. Fecal samples used for EPG evaluation were taken directly from the rectum of the animals and sent to the laboratory for analyses. The EPG count was assessed using the modified McMaster technique [27] and the parasites’ genera were identified using the morphometric key by Van Wyk et al. [28]. The feces were classified through visual inspection according to its consistency appearance (FC) developed by Gordon et al. [29] and modified by the authors considering three out of the original five-value scale: 1 for normal stool, 2 for pasty stool, and 3 for diarrheal feces.

The BCS was assessed after evaluating the prominences of the spinous and transverse bones of the spine, fat coverage, and muscle development between the last rib and the ileum wing, as described by Russel et al. [30]. Coloration of the ocular mucosa was measured by trained people by observing the medial part of the lower conjunctiva and comparing it with the famacha chart (FAM; FAMACHA© System) considering a five-value scale: 1 for red robust, 2 for rosy red, 3 for pink, 4 for pale pink, and 5 for white [31]. The descriptive statistics for the analyzed traits are presented in Table 1.

Table 1 Descriptive statistics for eggs per gram of feces (EPGlog), red blood cells (RBC), famacha (FAM), platelet (PLT), white blood cells (WBC), hematocrit (HCT), and hemoglobin (HGB)

Genotyping of animals

The extraction of the genomic DNA from each animal was performed from blood samples collected with EDTA. For this, red blood cells were lysed using 1 mL of whole blood and 500 μL of lysis solution (0.32 mol/L Sucrose, 12 mmol/L Tris-HCl pH 7.5, 5 mmol/L MgCl2, 1% Triton X-100) followed by centrifugation at 13,000 r/min. The supernatant was discarded to reach the leucocytes. The pellets were slowly washed adding 1 mL of ultrapure water and then the microtube was poured out of the water and held for a few minutes on absorbent paper to completely dry the pellet. This washing step was repeated approximately three times until a clean white pellet was obtained. The white pellet was then kept frozen (−20 °C) until sending to the DEOXI laboratory in Araçatuba-SP, where the DNA was completely extracted. The DNA of each sample was quantified and the degree of purity (ratio of optical densities between 260 and 280 nm) (Table 1) was evaluated. After these processes, the samples were stored at −20 °C until genotyping was performed.

A total of 576 animals were genotyped with the Ovine SNP12k BeadChip (Illumina, Inc.), that contained 12,785 biallelic SNP markers. Quality control consisted of excluding markers with unknown genomic position, located on sex chromosomes, monomorphic, with minor allele frequency (MAF) lower than 0.05, call rate lower than 90%, and with excess heterozygosity. After quality control, there were 11,602 SNPs and 574 samples left for analyses.

Quantitative genetic analyses

The variance components were estimated using a single animal trait model by single step genomic BLUP (ssGBLUP) procedure, under Bayesian inference [32]. For all studied traits the fixed effects considered in the model were contemporary groups (farm and management group), month of sample collection, sex, covariable body condition (linear effect), and age at the collection (linear and quadratic effect). All phenotypes were tested for data consistency and contemporary groups with less than three animals and records out of three standard deviations from the contemporary group mean were discarded, remaining 500 animals with phenotypes records. The ssGBLUP is a modified version of the animal model (BLUP) with additive relationship matrix A −1 replaced by H −1 [33]:

$$ {H}^{-\mathbf{1}}={A}^{-\mathbf{1}}+\left[\begin{array}{cc}\hfill \mathbf{0}\hfill & \hfill \mathbf{0}\hfill \\ {}\hfill \mathbf{0}\hfill & \hfill {G}^{-\mathbf{1}}-{A}_{\mathbf{22}}^{-\mathbf{1}}\hfill \end{array}\right] $$

where A 2 2 is a numerator relationship matrix for genotyped animals and G is the genomic relationship matrix created as described by [34]:

$$ G={ZDZ}^{\hbox{'}}q $$

where Z is the gene matrix containing allele frequency adjustment; D is the matrix that have the SNP weight (initially D = I); and, q is a weighting/standardization factor. According to [35], such factors can be obtained by ensuring that the G average diagonal is next to A 22 . The pedigree file has a total of 1196 animals and the last three generations of animals with records were considered. A linear model was used to analyze HCT, WBC, RBC, HGB, PLT, and EPGlog. The model can be represented by the following matrix equation:

$$ y=\mathrm{X}\upbeta + Za+e $$

where y is the observations vector; β is the vector of fixed effects; a is the additive direct vector; X is the incidence matrix; Z is the incidence genetic random effects additive direct matrix (the β vector associated with the y vector); e is the residual effect vector. The priori distributions of vectors y, a and e were given by:

$$ {\displaystyle \begin{array}{c}\hfill y\sim MVN\left( X\beta + Za\right)\hfill \\ {}\hfill a\mid G\sim MVN\left( 0,H\otimes G\right)\hfill \\ {}\hfill e\mid R\sim \mathrm{MVN}\left( 0,I\otimes R\right)\hfill \end{array}} $$

where H is the relationship coefficients matrix among animals obtained from the single-step analyses (single-step); R is the residual variance matrix; I is the identity matrix; G is genetic additive variance matrix and is the Kronecker product. An inverted chi-square distribution was used for the prior values of the direct and residual genetic variances. A uniform distribution was used a priori for the fixed effects.

A threshold model was applied to analyze the FAM trait. The scores of FAM for each individual i, were defined by Ui in the underlying scale yi = (1) t 0  < Ui < t 1 ; (2) t 1  < Ui < t 2 ; (3) t 2  < Ui < t 3 ; (4) t 3  < Ui < t 4 ; i = 1, … n, where n is the number of observations; t 1 to t 4 were the threshold values; and U is the unobservable continuous variable, in underlying scale, limited between two unobservable thresholds. After specifying the thresholds t 0 to t 4 it is necessary to adjust one of the thresholds (from t 1 to t 4 ) into an arbitrary constant. In the present study, it was assumed that t 1  = 0, in such a way that the vector of estimable thresholds was defined as t = t 2 , t 3 , t 4 . After defining the model parameters, the link between categorical and continuous scales could be established based on the contribution of the probability of an observation that fitted the first category, which is proportional to:

$$ P\left( yv=0|t,0\operatorname{}\right)=P\left( Uv<t|t,0\operatorname{}\right)=\Phi \left[\left(t-{w}^{'} vo\right)/{\sigma}_e\right] $$

where yv is the response variable for the vth observation; t is the threshold value arbitrarily assigned as the true value is unobservable; Uv is the value of the underlying variable for the vth observation; (ϕ) is the cumulative distribution function of a standard normal variable; and w’v is the scale of the incidence matrix that linked θ to the vth observation. As the observations are conditionally independent, given θ, the likelihood function is defined by the product of contributions of each record.

The analyses were performed using the GIBBS2F90 and THRGIBBSF90 programs [33, 36]. The a posteriori variance component estimates were obtained using the POSTGIBBSF90 program [36]. The analysis consisted of a single chain of 500,000 cycles discarding the first 20,000 cycles, taking a sample at every 100 iterations. Thus, 48,000 samples were used to obtain the parameters. The data convergence was verified with the graphical evaluation of sampled values versus interactions according to the criteria proposed by several authors [37,38,39], using the Bayesian Output Analysis (BOA) of the R 2.9.0 software [40].

Linkage disequilibrium estimation

The linkage disequilibrium (LD) between two SNPs was evaluated using r 2 as follows:

$$ {r}^2={\frac{\left( freq.{AB}^{\ast } freq. ab- freq.{Ab}^{\ast } freq. aB\right)}{\left( freq.{A}^{\ast } freq.{a}^{\ast } freq.{B}^{\ast } freq.b\right)}}^2=\frac{(D)^2}{\left( freq.{A}^{\ast } freq.{a}^{\ast } freq.{B}^{\ast } freq.b\right)} $$

where,

$$ D= freq. AB- freq.{A}^{\ast } freq.B $$

And

$$ {D}^{\hbox{'}}=\left\{\begin{array}{c}\hfill \frac{D}{\min \left( freq.{A}^{\ast } freq.b, freq.{a}^{\ast } freq.B\right)}\kern0.24em if\kern1em D<0\hfill \\ {}\hfill \frac{D}{\min \left( freq.{A}^{\ast } freq.B, freq.{a}^{\ast } freq.b\right)}\kern0.24em if\kern1em D<0\hfill \end{array}\right\} $$

where freq.A, freq.a, freq.B and freq.b are the frequencies of alleles A, a, B and b, respectively, and freq.AB, freq.ab, freq.aB and freq.Ab are the frequencies of haplotypes Ab, ab, aB and Ab in the population, respectively. The expected frequency of haplotype AB (freq.AB) is calculated as the product between freq.A and freq.B. The r 2 takes values close to zero when the alleles A and B segregate independently. The freq.AB higher or lower than the expected value indicates that these two loci in particular tend to segregate together and are in LD, with a maximum value for r 2 of one.

Principal component analysis

The principal component analyses (PCAs) were obtained from the genomic relationship matrix through the preGS90 program. The preGS90 is an interface program to process genomic information for the BLUPF90 family [36]. Efficient methods for the creation of the genomic relationship matrix, relationship based on genomic data, and their inverses are described by [41]. The PCAs applied to genotype data can be used to calculate principal components (PCs) that explain differences among individual samples in the genetic data. The top PCs are viewed as continuous axes of variation that reflect genetic variation due to ancestry in the sample. Individuals with similar values for a particular principal component will have a similar ancestry for that axes.

Genome-wide association analysis

The genome-wide association analysis for each trait was performed using the single-step GWAS (ssGWAS) methodology [42]. The same linear animal model for HCT, WBC, RBC, HGB, PLT and EPGlog, and the threshold model for FAM used to estimate the variance components were applied. The effects were decomposed in genotyped (ag) and non-genotyped (an) animals as describe by [42], considering the effect of genotyped animals as:

$$ {\mathrm{a}}_g=Z\mathrm{u}, $$

where Z is a matrix that relates genotypes of each locus and u is a vector of marker effects, and the variance of animal effects was assumed as:

$$ \operatorname{var}\left({\mathrm{a}}_g\right)=\operatorname{var}\left(\mathrm{Zu}\right)={\mathrm{ZDZ}}^{\hbox{'}}{\sigma}_u^2={G}^{\ast }{\sigma}_{\mathrm{a}}^2 $$

where D is a diagonal matrix of weights for variances of markers (D = I for GBLUP), \( {\sigma}_{\mathrm{u}}^2 \) is the genetic additive variance captured by each SNP marker when no weights are present, and G * is the weighted genomic relationship matrix. The ratio of covariance of genetic effects (a g ) and SNPs (u) is:

$$ \operatorname{var}\left[\begin{array}{c}\hfill {a}_t\hfill \\ {}\hfill \mathrm{u}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill {\mathrm{ZD}\mathrm{Z}}^{\hbox{'}}\hfill & \hfill {\mathrm{ZD}}^{\hbox{'}}\hfill \\ {}\hfill {\mathrm{DZ}}^{\hbox{'}}\hfill & \hfill D\hfill \end{array}\right]{\sigma}_{\mathrm{u}}^2, $$

Sequentially:

$$ {\mathbf{G}}^{\ast}\kern0.5em =\frac{\operatorname{var}\left({\mathrm{a}}_{\mathrm{g}}\right)}{\sigma_{\mathrm{a}}^2}=\frac{\operatorname{var}\left(\mathrm{Zu}\right)}{\sigma_{\mathrm{a}}^2}={\mathrm{ZDZ}}^{\hbox{'}}\lambda $$

where λ is a variance ratio or a normalizing constant. According to [34],

$$ \lambda =\frac{\sigma_{\mathrm{a}}^2}{\sigma_{\mathrm{a}}^2}=\frac{1}{\sum_{\mathrm{i}=1}^{\mathrm{M}}2{\mathrm{p}}_{\mathrm{i}}\left(1-{\mathrm{p}}_{\mathrm{i}}\right),} $$

where M is the number of SNP and p i is the allele frequency of the second allele of the i th SNP. According to [43], the markers effect can be described by:

$$ \widehat{\mathrm{u}}=\frac{\upsigma_u^2}{\upsigma_u^2}\mathrm{DZ}\hbox{'}\mathbf{G}{\ast}^{-1}{\widehat{\mathrm{a}}}_{\mathrm{g}}={\mathrm{DZ}}^{\hbox{'}}{\left[{\mathrm{ZDZ}}^{\hbox{'}}\right]}^{-1}{\widehat{\mathrm{a}}}_{{\mathrm{g}}^{\hbox{'}}} $$

The estimated SNP effects can be used to estimate the variance of each individual SNP effect [44] and apply a different weighting for each marker, such as:

$$ {\hat{\upsigma}}_{\hat{\mathrm{u}},\mathrm{i}}^2={\hat{\mathrm{u}}}_i^22{\mathrm{p}}_{\mathrm{i}}\left(1-{\mathrm{p}}_{\mathrm{i}}\right) $$

The following iterative process described by [42] was used considering D to estimate the SNP effects:

  1. 1.

    D = I;

  2. 2.

    To calculate the matrix G = ZDZ'qG = ZDZ q

3. To calculate GEBVs for all animals in data set using ssGBLUP;

4. To calculate the SNP effect: \( \widehat{\mathrm{u}}={\uplambda \mathrm{DZ}}^{'}{\mathrm{G}}^{\ast -1}{\widehat{\mathrm{a}}}_g \)

5. To calculate the variance of each SNP:\( {\mathrm{d}}_{\mathrm{i}}={\widehat{\mathrm{u}}}_{\mathrm{i}}^22{\mathrm{p}}_{\mathrm{i}}\left(1-{\mathrm{p}}_{\mathrm{i}}\right) \), where I is the i-th marker;

6. To normalize the values of SNPs to keep constant the additive genetic variance;

7. Exit, or loop to step 2.

The markers effects were obtained by two iterations from step 2 to 7. The percentage of genetic variance explained by i-th region was calculated as described by [45].

$$ \frac{\mathrm{Var}\left({\mathrm{a}}_{\mathrm{i}}\right)}{\upsigma_{\mathrm{a}}^2}\times 100=\frac{\mathrm{Var}\left({\sum}_{\mathrm{j}=1}^{10}{\mathrm{Z}}_{\mathrm{j}}{\widehat{\mathrm{u}}}_{\mathrm{j}}\right)}{\upsigma_{\mathrm{a}}^2}\times 100 $$

where a i is genetic value of the i-th region that consists of 10 consecutive SNPs, \( {\upsigma}_a^2 \) is the total genetic variance, Z j Z j is vector of gene content of the j-th SNP for all individual, and \( {\widehat{\mathrm{u}}}_j \) \( {\widehat{u}}_j \) is marker effect of the j-th within the i-th region. Considering a length of the windows size not exceeding 3 to 4 Mb, the results were presented as the proportion of genetic additive variance explained by windows of 10 consecutive SNPs.

Search for genes

The chromosome regions that explained more than 1.0% of the additive genetic variance were selected to explore and determine possible quantitative trait loci (QTL). The Map Viewer tool of ovine (Ovis aries) genome was used for identification of genes, available at “National Center for Biotechnology Information” (NCBI - http://www.ncbi.nlm.nih.gov) database, using the bank references of HuRef assembly, CHM1 1.0, CRA TCAGchr7v2 and Ensembl Genome Browser (http://www.ensembl.org/index.html). The classification of the genes according to its biological function, identification of metabolic pathways and gene enrichment was performed by DAVID tool v6.8 [46] and GeneCards (http://www.genecards.org/). Gene ontology (GO; P < 0.01) and Kyoto Encyclopedia of Genes and Genomes (KEGG; P < 0.05) pathways were identified by DAVID tool v6.8 [46]. All annotated genes in the Ovis aries genome were used as background.

Results and discussion

Linkage disequilibrium and population structure

The Fig. 1 represents the first two principals (PC1 and PC2) component analysis (PC) based on the genomic relationship matrix. Despite the animals being originated from flocks located in different regions of Brazil, it is important to note that there was no subpopulation structure in this population.

Fig. 1
figure 1

First two principal components (pc1: 3.16%;pc2:15%) of the genomic relationship matrix of genotyped animals

The predicted values of LD vs. linkage distance between genetic markers were presented in Fig. 2. On the basis of this figure, it is possible to state that when considering a distance between markers lower than 1 Mb, the level of LD indicates that the Ovine SNP12k BeadChip may be a suitable tool for identifying genomic regions associated with traits related to gastrointestinal parasites resistance. Most tightly linked SNP pairs have the highest r 2 and average r 2 rapidly decreases as linkage distance increases (Fig. 2). The overall LD mean between markers considering windows up to 10 Mb was 0.07. The LD mean between adjacent SNPs across autosomes ranged from 0.02 to 0.10 (Table 2). Several authors have study the pattern of LD between markers in the genome of various sheep breeds [47,48,49,50,51,52]. It is difficult to compare the level of LD obtained in different studies since different sample sizes, LD measures, marker types, marker densities and historical population demographics [51] were used. However, the level of LD obtained in the present study was similar to those reported in previous studies for wool sheep breeds using the Ovine SNP50 BeadChip [50,51,52].

Fig. 2
figure 2

Linkage disquilibrium (r 2) between markers considering windows up yo 10 Mb

Table 2 Summary of SNP markers analyzed and average linkage disequilibrium (LD; r 2) between synthetic adjacent markers obtained for each autosome (OAR)

Genetic parameters estimates

The descriptive statistics for FAM indicates an incidence of animals with some degree of anemia (mainly levels 4 and 5) and animals that are not affected (levels 1 and 2) (Table 1). In the present study, animals were free of other sources that lead to anemia, i.e. fluke or minerals deficiencies such as copper, so, it is possible to asseverate that the presence of anemia in these animals was probably due to the prevalence of a Haemonchus population and the susceptibility or not of these ewes to infection. Similarly, the descriptive statistics for EPGlog indicates an infection by gastrointestinal parasites (Table 1). Most of the gastrointestinal parasites belonging to the Trichostrongylidae family have similar shape and eggs size, and unless the EPGlog is accompanied by a larvae cultivar, the eggs can belong to any specie. Moreover, the larva 4 and 5 of Haemonchus develops in the wall of the stomach with two particularities: (i) they feed with blood and (ii) do not lay eggs, which make FAM more precise than EPGlog to predict infection by Hemonchus contortus in all damaging stages.

The variance components and heritability estimates for EPGlog, RBC, FAM, WBC, PLT, HCT, and HGB were described in Table 3. The parameter estimates convergence was verified by inspection of trace-plots [37, 38] in which the convergence diagnosis indicates the convergence of the chain. Thus, the burn-in period used was considered enough to achieve the convergence of the estimates of all parameters. The marginal posterior distributions of heritability estimates for the traits showed accurate values, tending to a normal distribution, since the mean and the median were similar (Table 3). Symmetric distributions of central tendency statistics were an indicative that the analysis is reliable.

Table 3 Estimates of additive genetic variance (Va), residual variance (Vr), and heritability (h 2) for log10 of eggs per gram of feces (EPGlog), red blood cells (RBC), famacha (FAM), platelets (PLT), white blood cells (WBC), hematocrit (HCT), and hemoglobin (HGB)

Heritability estimates were low for EPGlog, moderate for RBC, PLT, HCT, HGB and WBC, and high for FAM. Studies involving Santa Inês breed are scarce in the literature and most of them refers to genetic parameters for traits related to parasites resistance performed in wool sheep breeds. Bisset et al. [53] reported heritability estimate for resistance and resilience of sheep to Haemonchus contortus in a South African Merino flock and observed high heritability for EPG, FAM and HCT, with values ranging from 0.47 to 0.55. Contrasting with our results, moderate heritability estimates for EPG among different studies were reported. Bishop et al. [54] studying Texel lambs observed a weighted average heritability of 0.26 and 0.38 for Strongyle and Nematodirus nematode resistance, respectively. Pickering et al. [25] observed an estimate of 0.25 in New Zealand sheep, and more recently, Benavides et al. [55] reported a heritability estimate of 0.36 in Australian Merino flock. Riley et al. [56] found a low heritability estimate for FAM in a Merino flock in South Africa, with values ranging from 0.06 to 0.24. The heritability estimate obtained for FAM pointed out that genetic progress for this trait is feasible, and this trait should be included in Santa Inês breeding programs. Moreover, the FAM has low cost and it is easily measured. In general, all traits showed genetic variability, therefore, it is important to investigate the presence of genomic regions or genes affecting these traits so as to elucidate and better understand their genetic architecture, especially for Haemonchus contortus, since it is one of the most predominant, highly pathogenic and economically important gastrointestinal parasite in sheep 5.

GWAS, genomic regions and enrichment analysis

The SNPs windows regions which accounted for more than 1% of the genetic additive variance were used to search for candidate genes (CG), which are described in Additional file 1. A total of 22, 21, 23, 20, 26, 25 and 23 windows for EPGlog, FAM, WBC, RBC, PLT, HCT and HGB traits were identified, respectively. Among the associated windows, 10 are common for HCT and HGB traits, located on OAR1, OAR2, OAR3, OAR5, OAR8 and OAR15 (Additional file 1: Table S9 and S10 – Figs. 3 and 4).

Fig. 3
figure 3

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for hematocrit (HCT) trait

Fig. 4
figure 4

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for hemoglobin (HGB) trait

The genomic regions associated with EPGlog are detailed in Additional file 1: Table S4 and illustrated in Fig. 5. Genes such as CYP11A1 and CYP1A1 located on chromosome 18, and CYP19 and SFXN1 located on chromosome 7 have functions associated with transportation or construction of iron molecules, and its absence in the blood can be an indicator of anemia in animals. The B2M, SFXN1, IL25, BMP4, TSHR, CCL28, PIK3R1, FGF10, IL15, IL2, TP-1, BPMG, BCL10, HSPD1, and MALT1 genes are described to have functions related to the body’s immune and defense response. Indeed, these genes participate in metabolic pathways related to the development of the immune system and in the regulation of the immune process effects. Atlija et al. [26] also observed that the IL25 gene is involved in EPG trait in adult sheep, with functions also linked to immune response. The enrichment analysis revealed the PI3K-Ark signaling pathway, which controls several cellular responses and has important functions in the immune system by regulating many key events in the inflammatory response related to damage and infection [57]. The genomic regions associated with FAM are described in Additional file 1: Table S5 and illustrated in Fig. 6. The identified CGs have functions related to the development of the immune system and the body’s defense response. Genes, such as ADAM10, IL6ST, TNFRF13B, SIVA1, JUN, PAX1, PIK3R1, SIT1, and AKT1 are related to the differentiation of T cells, a group of white blood cells (leukocytes) responsible for defending the body against antigens. The GUCY1A2 gene, which interacts selectively and non-covalently with iron Fe and the SIVA1 gene, both are expressed in different subpopulations of T and B cells and provide co-stimulatory signals for the proliferation of T and B cells and immunoglobulin production by B cells [58]. The GUCY1A2 gene operates in the process of regulating the proliferation and elimination of T cells, maintaining its number stable in the absence of external stimulus.

Fig. 5
figure 5

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for egg count per gram of feces (EPGlog) trait

Fig. 6
figure 6

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for famacha (FAM) trait

The results of the enriched genes and functional grouping analyses showed that genes associated with FAM are involved in functions related to body’s immune response. The biological processes (Additional file 2 – Table S11 and S12) showed that these genes are significantly (P < 0.05) involved in lymphocytes and leukocytes homeostasis.

The CCL28 gene located on OAR16 identified both for EPG and FAM acts as a chemotactic for CD4 and CD8 inactive T cells [59]. The metabolic pathways associated with the CCL28 gene and with others associated with EPGlog and FAM traits, whose functions were found related to immunoglobulin synthesis in the intestine, are presented in Fig. 7.

Fig. 7
figure 7

Metabolic pathway involving genes present for egg count per gram of feces (EPG) and famacha (FAM) traits

The genomic regions associated with WBC are described in Additional file 1: Table S6 and illustrated in Fig. 8. The genes found in the SNP windows encode proteins that are involved in many immune system processes, responding to any potential internal or invasive threat. The vertebrate’s immune system is composed by the innate immune system (defense responses mediated by germline encoded components that directly recognize components of potential pathogens) and by the adaptive immune system which consists of T- and B-lymphocytes [60]. Many different innate immune mechanisms are deployed for host defense, a unifying theme of innate immunity is the use of germline-encoded pattern recognition receptors for pathogens or damaged self-components, such as the toll-like receptors, nucleotide-binding domain leucine-rich repeat (LRR)-containing receptors, retinoic acid-inducible gene I-like RNA helicases and C-type lectin receptors [61], like the CRCP (CGRP receptor component (CRCP)).

Fig. 8
figure 8

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for white blood cells (WBC) trait

Some of these genes (Additional file 1: Table S6) have already been documented in others species and been associated with immune response. Beside the functions described above, they also present some particularities, i.e. the EPBH1 (EPH Receptor B1), which is involved in the GPCR pathway and developmental biology. The first process involves the G protein–coupled receptors (GPCRs) that constitute a large protein family of receptors that recognize molecules outside the cell and activate inside signal transduction pathways and, finally, cellular responses [62].

The GATA3 (GATA binding protein 3) gene also acts on the defense response, reacting to the presence of a foreign body or to the occurrence of an injury, restricting the injury/damage extension or preventing/recovering it from the infection caused by the foreign body. The GATA3 gene also positively regulates T cell differentiation and any process that activates or increases the frequency, rate or extent of T cell differentiation. An additional gene having a peculiar function is the SAMHD1 (SAM domain and HD domain 1), which is also associated with the regulation of the innate immune responses. The SLA2 (Src-like-adaptor 2) gene acts in the T cell regulation and negative regulation of B cell activation, and in any process that stops, prevents, or reduces the frequency, rate or extent of B cell activation.

The genes in the enrichment pathways (Additional file 2), such as the MYPN (myopalladin), PDGFRL (platelet derived growth factor receptor like), ROR1 (receptor tyrosine kinase-like orphan receptor 1), CNTN1 (Contactin 1) and FANK1 (fibronectin type 3, and ankyrin repeat domains 1) are associated with the Immunoglobulin-like domains that are related in both sequence and structure, and can be found in several diverse protein families. Ig-like domains are involved in a variety of functions, including cell-cell recognition, cell-surface receptors, muscle structure and the immune system [63].

Additional file 1: Table S7 and S8 show the genes with functions associated with immune response, innate, and adaptative immune response or those participating in the regulation of innate immune response for the genomic regions associated with RBC and PLT traits, respectively

The genomic regions and genes associated with RBC are presented in Additional file 1: Table S7 and illustrated in Fig. 9. The KEGG pathays analysis revealed PI3K-Ark and toll-like receptor signaling pathways significantly enriched for WBC trait (Additional File 2). The CD109 gene (CD109 antigen) is a GPI-linked cell surface antigen expressed by CD34+ acute myeloid leukemia cell lines, T-cell lines, activated T lymphoblasts, endothelial cells, and activated platelets [64]. This gene was also reported by [26], indicating that CD109 gene potentially contributes to resistance to gastrointestinal nematodes in sheep.

Fig. 9
figure 9

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for red blood cells (RBC) trait

The IL2RB gene (interleukin 2 receptor), which is involved in T cell-mediated immune responses, is present in three forms with respect to ability to bind interleukin 2. The low affinity form is a monomer of the alpha subunit and it is not involved in signal transduction. The intermediate affinity form consists of an alpha/beta subunit heterodimer, while the high affinity form consists of an alpha/beta/gamma subunit heterotrimer. Both the intermediate and high affinity forms of the receptor are involved in receptor-mediated endocytosis and transduction of mitogenic signals from interleukin 2 [provided by RefSeq, Jul 2008]. Kondo et al. [65] showed that a clonogenic common lymphoid progenitor, a bone marrow-resident cell that gives rise exclusively to lymphocytes (T, B, and natural killer (NK) cells), can be redirected to the myeloid lineage by stimulation through exogenously expressed interleukin (IL)-2 and GM-CSF (granulocyte/macrophage colony-stimulating factor) receptors. As the IL2RB is one of the receptors that act in IL-15 expression, it has been proposed to be a critical cytokine for NK cell development. This gene can affect this redirections in cytokine signaling and can regulate cell-fate decisions. Additionally, a critical step in lymphoid commitment is downregulation of cytokine receptors that drive myeloid cell development.

Poliovirus receptor-like proteins (PVRLs), such as PVRL4, are adhesion receptors of the immunoglobulin superfamily and function in cell-cell adhesion [66]. The encoded protein contains two immunoglobulin-like (Ig-like) C2-type domains and one Ig-like V-type domain. It is involved in cell adhesion through trans-homophilic and -heterophilic interactions. It is a single-pass type I membrane protein [provided by RefSeq, Jan 2011].

The CXCR4 gene may influence the immune system under physiologic and pathologic conditions through negative regulation of MHC class II expression, possibly through PKA and SRC kinase [67]. In a study, Contento et al. [68] observed that while human T cell activation by antigen-presenting cells is taking place, the CCR5 and CXCR4 chemokine receptors are recruited into the immunological synapse, where they deliver costimulatory signals. This gene also participates in the intestinal immune network for IgA production pathway.

The ASCC3, CRLF2RB, FNDC3A, FNDC4, IFNAR1, IFNAR2, LEPR, MYLK, NCAM2, PVRL4 SLAMF1, TXNIP, CD109, ACAN, NTRK3, CCDC141, LEPR, MYLK, PVLR4, IL12RB2, IGSF10, and JAM2 genes (Additional file 1: Table S7) are domains with an Ig-like fold that can be found in several proteins in addition to immunoglobulin molecules. For example, Ig-like domains occur in several different types of receptors (such as various T-cell antigen receptors), several cell adhesion molecules, MHC class I and II antigens, as well as the hemolymph protein hemolin, and the muscle proteins titin, telokin and twitchin (IPR013783).

The KEGG pathway analysis revealed platelet activation, toll-like receptor and PI3K-Ark signaling pathways for PLT (Additional file 2). The CXCL1, CXCL8 and CXCL10 genes identified for PLT (Additional file 1: Table S8; Fig. 10) belong to a subfamily of chemokines that basically are structurally related to molecules that regulate cell trafficking of various types of leukocytes through interactions with a subset of 7-transmembrane and G protein-coupled receptors. They also play a fundamental role in the development, homeostasis, and functionality of the immune system [69].

Fig. 10
figure 10

Manhattan plot of the additive genetic variance explained by windows of 10 adjacent SNPs for platelet (PLT) trait

Different genes were also identified for PLT, in these regard, the Fc-gamma receptors (FCGRs), such as FCGR1A, FCGR3A, and FCGR2B, that are integral membrane glycoproteins which exhibit complex activation or inhibitory effects on cell functions after aggregation by complexed immunoglobulin G (IgG) [70]. These genes encode proteins that play an important role in the immune response, a receptor for the Fc portion of immunoglobulin G, and participate in the removal of antigen-antibody complexes from the circulation, as well as other antibody-dependent responses. It is also involved in the phagocytosis of immune complexes and in the regulation of antibody production by B-cells, respectively [provided by RefSeq, Jul 2008; RefSeq, Jun 2010].

Atlija et al. [26] observed some significant chromosome-wise QTL detected by linkage analysis and from the combined LD and linkage analysis. These authors reported some CG involved in immune response. Theirs findings support our results, since genes such as MBN, BTC, CXCL1, CXCL10, EREG, RASSF6, SCARB2, TMPRSS11D, AMBN, and AREG were observed in both studies, indicating that these genes can be used and deply studied as CG for resistance to gastrointestinal nematodes in sheep. Some of the genes cited in Additional file 1: Table S9 and S10 have functions related to immune response or to immunoglobulin development, and were associated with HCT and HBG traits, respectively.

For HCT trait (Table S9), genes like IGF1R, IFNAR1, IFNAR2, IL12RB2, IL2RB, C2F2RB, PVRL4, LEPR, MYLK, SLAMF1, ADGRF3, and TXNIP were identified. These genes encode domains with an Ig-like fold that can be found in several proteins in addition to immunoglobulin molecules. For example, Ig-like domains occur in several different types of receptors (such as T-cell antigen receptors), many cell adhesion molecules, MHC class I and II antigens, as well as the hemolymph protein hemolin, and the muscle proteins titin, telokin and twitchin.

The IL1A and IL1B genes have the protein encoded as a member of the interleukin 1 cytokine family. This cytokine is a pleiotropic cytokine involved in several immune responses, inflammatory processes, and hematopoiesis [provided by RefSeq, Jul 2008]. Also the IL1B gene, initially discovered as the major endogenous pyrogen, induces prostaglandin synthesis, neutrophil influx and activation, T- and B-cells activation, cytokine production, antibody production, fibroblast proliferation and collagen production. Also it was observed a function related to ligand-binding domain that displays similarity to C2-set immunoglobulin domains (antibody constant domain 2) in the LEPR gene. It has a specific effect on T lymphocyte responses, differentially regulating the proliferation of naive and memory T-cells.

For HCT and HGB traits, genes related to immune response were identified, i.e. the RAC2 gene, that seems to increase the resistance to parasites. In order to understand the function of RAC2 GTPase in regulating the cellular immune response, Williams et al. [71] assayed the effects of hemocytes in parasitized RAC2 mutant larvae, as well as characterized the effect of over-expressing RAC2 gene in hemocytes. These authors reported that this gene has an important role in the cellular immune response, being necessary for hemocyte spreading and cell-cell contact formation during immune surveillance against the parasitoid L. boulardi. When an invading parasitoid is recognized as a foreign body, circulating hemocytes should recognize them and attach to the egg chorion. After the attachment, the hemocytes should then spread out to form a multilayered capsule surrounding the invader. In RAC2 mutants this process was disrupted. Besides of that, this gene also acts on the B cell receptor signaling pathway, like the PPP3CA gene.

The HSPD1 gene encodes a mitochondrial protein that plays a role as a signaling molecule in the innate immune system. Zanin-Zhorov et al. [72] reported that this gene, as well as a synthetic peptide derived from it, acted as a costimulatory of human regulatory CD4-positive/CD25 (IL2RA) positive T cells (Tregs), which inhibit lympho proliferation and IFNG and TNF secretion by CD4-positive and CD8-positive T cells. The authors concluded that the self-molecule HSPD1 can downregulate adaptive immune responses by upregulating Tregs through TLR2 signaling.

The CD80 and CD86 genes identified in almost all traits (RBC, PLT and HCT) probably participate effectively in the activation of T cells, which requires engagement of two separate T-cell receptors. The antigen-specific T-cell receptor (TCR) binds foreign peptide antigen-MHC complexes, and the CD28 receptor binds to the B7 (CD80/CD86) costimulatory molecules expressed on the surface of antigen-presenting cells (APC). The simultaneous triggering of these T-cell surface receptors with their specific ligands results in the activation of this cell. Many in vitro and in vivo studies demonstrated that both CD80 and CD86 ligands have an identical role in the activation of T cells. Recently, functions of B7 costimulatory molecules in vivo have been investigated in B7-1 and/or B7-2 knockout mice, and the authors concluded that CD86 could be more important for initiating T-cell responses, while CD80 could be more significant for maintaining these immune responses [73]. Recently, [24] observed some important genes matching to important pathways involved in host immune response against parasites. The PDGFRA gene (OAR6: 67,950,121–69,892,816) was also observed by Benavides et al. [24], whose findings have pointed out as a key gene in cytokine signaling.

Conclusions

The traits indicating gastrointestinal parasites resistance shown adequate genetic variability to respond to selection in Santa Inês breed, and it is expected a higher genetic gain for famacha when compared to the others traits. The level of LD estimated for markers separated by less than 1 Mb indicates that the Ovine SNP12k BeadChip will likely be a suitable tool for identifying genomic regions associated with those traits related to gastrointestinal parasite resistance.

Several candidate regions related to immune system development and activation, inflammatory response, regulation of lymphocytes, and leukocytes proliferation were found in this study. These candidate regions and CG may help in the selection of animals with higher resistance to parasites, and consequently, reduce the anthelmintic drugs costs, also the production losses linked to the use of it. Furthermore, when reducing the use of anthelmintic drugs, a potential reduction in waste problems regarding meat and milk discharge due to drugs residues should be minimized, reducing the impact upon the environment.

Abbreviations

BCS:

Body condition score

EPGlog :

Egg counts per gram of feces

FAM:

Degree of anemia assessed by the famacha card

FC:

Feces consistency

FS:

Fur score

HCT:

Hematocrit

HGB:

Hemoglobin

LD:

Linkage disequilibrium

PLT:

Platelets

RBC:

Red blood cell count

WBC:

White blood cell count

References

  1. Zygoyiannis D. Sheep production in the world and in Greece. Small Rumin Res. 2006;62:143–7.

    Article  Google Scholar 

  2. Perry BD, Randolph TF, Mcdermott JJ, Sones KR, Thornton PK. Investing in Animal Health Research to Alleviate Poverty. ILRI (International Livestock Research Institute); 2002.

  3. Liu J, Zhang L, Xu L, Ren H, Lu J, Zhang X, et al. Analysis of copy number variations in the sheep genome using 50K SNP BeadChip array. BMC Genomics. 2013;14:229.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  4. Da Silva M V, Lopes PS, Guimarães SE, De R, Torres A. Utilização de marcadores genéticos em suínos. I. Características reprodutivas e de resistência a doenças The use of genetic markers in swine. I. Reproductive and disease resistance traits. 2002;

  5. Amarante AFT, Bricarello PA, Rocha RA, Gennari SM. Resistance of Santa Ines, Suffolk and Ile de France sheep to naturally acquired gastrointestinal nematode infections. Vet Parasitol. 2004;120:91–106.

    Article  CAS  PubMed  Google Scholar 

  6. O’Connor LJ, Walkden-Brown SW, Kahn LP. Ecology of the free-living stages of major trichostrongylid parasites of sheep. Vet Parasitol. 2006;142:1–15.

    Article  PubMed  Google Scholar 

  7. Rocha RA, Amarante AFT, Bricarello PA. Comparison of the susceptibility of Santa Inês and Ile de France ewes to nematode parasitism around parturition and during lactation. Small Rumin Res. 2004;55:65–75.

    Article  Google Scholar 

  8. Bricarello PA, Amarante AFT, Rocha RA, Cabral Filho SL, Huntley JF, Houdijk JGM, et al. Influence of dietary protein supply on resistance to experimental infections with Haemonchus contortus in Ile de France and Santa Ines lambs. Vet Parasitol. 2005;134:99–109.

    Article  CAS  PubMed  Google Scholar 

  9. Costa RLD, Bueno MS, Veríssimo CJ, Cunha EA, Santos LE, Oliveira SM, et al. Performance and nematode infection of ewe lambs on intensive rotational grazing with two different cultivars of Panicum maximum. Trop Anim Health Prod. 2007;39:255–63.

    Article  CAS  PubMed  Google Scholar 

  10. Woolaston RR, Piper LR. Selection of Merino sheep for resistance to Haemonchus contortus: genetic variation. Anim Semit. 2016;62:451–60.

    Article  Google Scholar 

  11. Woolaston RR, Windon RG. Selection of sheep for response to Trichostrongylus colubriformis larvae: genetic parameters. Anim Sci. 2016;73:41–8.

    Article  Google Scholar 

  12. 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. New Zeal J Agric Res. 2005;48:10.

    Google Scholar 

  13. 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 

  14. Mcewan JC. Developing genomic resources for whole genome selection. Proceedings of the New Zealand Society of Animal Production. 2007. p. 148–53.

  15. Kemper KE. The implications for the host-parasite relationship when sheep are bred for enhanced resistance to gastrointestinal nematodes. PhD thesis. Melbourne: The University of Melbourne; 2010.

    Google Scholar 

  16. Pickering NK, Blair HT, Hickson RE, Dodds KG, Johnson PL, McEwan JC. Genetic relationships between dagginess, breech bareness, and wool traits in New Zealand dual-purpose sheep. J Anim Sci. 2013;91:4578–88.

    Article  CAS  PubMed  Google Scholar 

  17. Mcrae KM, Stear MJ, Good B, Keane OM. The host immune response to gastrointestinal nematode infection in sheep. Parasite Immunol. 2015;37:605–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Miller JE, Horohov DW. Immunological aspects of nematode parasite control in sheep. J Anim Sci. 2006;84(Suppl):124–32.

    Article  Google Scholar 

  19. Bolormaa S, van der Werf JHJ, Walkden-Brown SW, Marshall K, Ruvinsky A. A quantitative trait locus for faecal worm egg and blood eosinophil counts on chromosome 23 in Australian goats. J Anim Breed Genet. 2010;127:207–14.

    Article  CAS  PubMed  Google Scholar 

  20. Alba-Hurtado F, Muñoz-Guzmán MA. Immune Responses Associated with Resistance to Haemonchosis in Sheep. Biomed Res. Int. 2013. doi: 10.1155/2013/162158.

  21. 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.

    Article  CAS  PubMed  Google Scholar 

  22. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  23. 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.

    Article  PubMed  PubMed Central  Google Scholar 

  24. 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 x Dorper backcross population. PLoS One. 2015;10:1–20.

    Google Scholar 

  25. Pickering NK, Auvray B, Dodds KG, Mcewan JC. Genomic prediction and genome-wide association study for dagginess and host internal parasite resistance in New Zealand sheep. BMC Genomics. 2011;16:1–11.

    Google Scholar 

  26. Atlija M, Arranz J-J, Martinez-Valladares M, Gutiérrez-Gil B. Detection and replication of QTL underlying resistance to gastrointestinal nematodes in adult sheep using the ovine 50K SNP array. Genet Sel Evol. 2016;48:4.

    Article  PubMed  PubMed Central  Google Scholar 

  27. Ueno H, Gonçalves PC. Manual para diagnóstico das helmintoses de ruminantes. 1998. p. 149.

  28. Van Wyk JA, Cabaret J, Michael LM. Morphological identification of nematode larvae of small ruminants and cattle simplified. Vet Parasitol. 2004;119:277–306.

    Article  PubMed  Google Scholar 

  29. Gordon H McL. The diagnosis of helminthosis in sheep. Med. Vet. Rev.1967; 67(140–168).

  30. Russel AJF, Doney JM, Gunn RG. Subjective assessment of body fat in live sheep. J Agric Sci, Camb. DigiTop - USDA’s Digital Desktop Library. 2016;72:451–4.

  31. Bath GF, Hansen JW, Krecek RC, Van Wyk JA, Vatta A. Sustainable approaches for managing haemonchosis in sheep and goats. Final Report of Food and Agricultural Organization (FAO) Technical Cooperation Project in South Africa. Food and Agriculture Organization of the United Nations; 2001.

  32. Gianola D, Fernando RL. Bayesian methods in animal breeding theory. J Anim Sci. 1986;63:217–44.

    Article  Google Scholar 

  33. Aguilar I, Misztal I, Johnson DL, Legarra A, Tsuruta S, Lawlor TJ. Hot topic: a unified approach to utilize phenotypic, full pedigree, and genomic information for genetic evaluation of Holstein final score. J Dairy Sci. 2010;93:743–52.

    Article  CAS  PubMed  Google Scholar 

  34. VanRaden PM, Van Tassell CP, Wiggans GR, Sonstegard TS, Schnabel RD, Taylor JF, et al. Invited review: reliability of genomic predictions for north American Holstein bulls. J Dairy Sci. 2009;92:16–24.

    Article  CAS  PubMed  Google Scholar 

  35. Vitezica ZG, Aguilar I, Misztal I, Legarra A. Bias in genomic predictions for populations under selection. Genet. Res. (Camb). 2011;93:357–66.

    Article  CAS  PubMed  Google Scholar 

  36. Misztal I, Tsuruta S, Strabel T, Auvray B, Druet T, Lee DH. BLUPF90 and related programs (BGF90). Proc. 7th World Congr. Genet. Appl. to Livest. Prod. 2002;28:21–2

  37. Geweke J. Evaluating the accuracy of sampling-based approaches to the calculation of posterior moments. In: Bernardo JM, Berger JO, David AP, Smith AFM, editors. Bayesian statistics. New York: Oxford University; 1992. p. 625-31. Cap. 4

  38. Heidelberg P, Welch P. Simulation run length control in the presence of an initial transient. Oper Res. 1983;31:1109–14.

    Article  Google Scholar 

  39. Raftery AE, Lewis S. How many iterations in the Gibbs sampler? Bayesian Stat. 1992:763–73.

  40. R Development Core Team. R: A language and environment for statistical computing. Vienna: R Foundation for Statistical Computing; 2008. ISBN 3–900051–07-0. URL http://www.R-project.org. Accessed 25 Aug 2016.

  41. Aguilar I, Misztal I, Legarra A, Tsuruta S. Efficient computation of the genomic relationship matrix and other matrices used in single-step evaluation. J Anim Breed Genet. 2011;128:422–8.

    Article  CAS  PubMed  Google Scholar 

  42. Wang H, Misztal I, Aguilar I, Legarra A, Muir WM. Genome-wide association mapping including phenotypes from relatives without genotypes. Genet Res (Camb). 2012;94:73–83.

    Article  CAS  Google Scholar 

  43. Strandén I, Garrick DJ. Technical note: derivation of equivalent computing algorithms for genomic predictions and reliabilities of animal merit. J Dairy Sci. 2009;92:2971–5.

    Article  PubMed  Google Scholar 

  44. Zhang Z, Liu J, Ding X, Bijma P, de Koning DJ, Zhang Q. Best linear unbiased prediction of genomic breeding values using a trait-specific marker-derived relationship matrix. PLoS One. 2010;5:1–8.

    Google Scholar 

  45. Wang H, Misztal I, Legarra A. Differences between genomic-based and pedigree-based relationships in a chicken population, as a function of quality control and pedigree links among individuals. J Anim Breed Genet. 2014;131:445–51.

    Article  CAS  PubMed  Google Scholar 

  46. Huang DW, Sherman BT, Lempicki RA. Bioinformatics enrichment tools: paths toward the comprehensive functional analysis of large gene lists. Nucleic Acids Res. 2009;37:1–13.

    Article  Google Scholar 

  47. McRae AF, McEwan JC, Dodds KG, Wilson T, Crawford AM, Slate J. Linkage disequilibrium in domestic sheep. Genetics. 2002;160:1113–22.

    CAS  PubMed  PubMed Central  Google Scholar 

  48. McRae AF, Pemberton JM, Visscher PM. Modeling linkage disequilibrium in natural populations: the example of the Soay sheep population of St. Kilda, Scotland. Genetics. 2005;171:251–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Meadows JRS, Chan EKF, Kijas JW. Linkage disequilibrium compared between five populations of domestic sheep. BMC Genet. 2008;9:61.

    Article  PubMed  PubMed Central  Google Scholar 

  50. García-Gámez E, Sahana G, Gutiérrez-Gil B, Arranz J-J. Linkage disequilibrium and inbreeding estimation in Spanish Churra sheep. BMC Genet. 2012;13:43.

    Article  PubMed  PubMed Central  Google Scholar 

  51. Mastrangelo S, Di Gerlando R, Tolone M, Tortorici L, Sardina MT, Portolano B. Genome wide linkage disequilibrium and genetic structure in Sicilian dairy sheep breeds. BMC Genet. 2014;15:108.

    Article  PubMed  PubMed Central  Google Scholar 

  52. Zhao F, Wang G, Zeng T, Wei C, Zhang L, Wang H, et al. Estimations of genomic linkage disequilibrium and effective population sizes in three sheep populations. Livest Sci. 2014;170:22–9.

    Article  Google Scholar 

  53. Bisset SA, Van Wyk JA, Bath GF, Morris CA, Stenson MO, Malan FS. Phenotypic and genetic relationships amongst FAMACHA score, faecal egg count and performance data in Merino sheep exposed to Haemonchus contortus infection in South Africa. In: Proceedings of the 5th International Sheep Veterinary Congress, Stellenbosch, 2011.

  54. Bishop SC, Jackson F, Coop RL, Stear MJ. Genetic parameters for resistance to nematode infections in Texel lambs and their utility in breeding programmes. Anim Sci. 2004;78:185–94.

    Google Scholar 

  55. Benavides MV, de Souza CJH, Moraes JCF, Berne MEA. Is it feasible to select humid sub-tropical Merino sheep for faecal egg counts? Small Rumin Res. 2016;137:73–80.

    Article  Google Scholar 

  56. Riley DG, Van Wyk JA. Genetic parameters for FAMACHA© score and related traits for host resistance/resilience and production at differing severities of worm challenge in a Merino flock in South Africa. Vet Parasitol. 2009;164:44–52.

    Article  CAS  PubMed  Google Scholar 

  57. Hawkins PT, Stephens LR. PI3K signalling in inflammation. Biochim Biophys Acta - Mol Cell Biol Lipids. 1851;2015:882–97.

    Google Scholar 

  58. Prasad KVS, Ao Z, Yoon Y, Wu MX, Rizk M, Jacquot S, et al. CD27, a member of the tumor necrosis factor receptor family, induces apoptosis and binds to Siva, a proapoptotic protein. Immunology. 1997;94:6346–51.

    CAS  Google Scholar 

  59. Wang W, Soto H, Oldham ER, Buchanan ME, Homey B, Catron D, et al. Identification of a novel chemokine (CCL28), which binds CCR10 (GPR2). J Biol Chem. 2000;275:22313–23.

    Article  CAS  PubMed  Google Scholar 

  60. Junior Janeway CA, Travers P, Walport M, Shlomchik MJ. Immunobiology. 5th ed. New York: Garland Science; 2001.

    Google Scholar 

  61. Cooper MD, Herrin BR. How did our complex immune system evolve? Nat Rev Immunol. 2010;10:2–3.

    Article  CAS  PubMed  Google Scholar 

  62. Trzaskowski B, Latek D, Yuan S, Ghoshdastider U, Debinski A, Filipek S. Action of molecular switches in GPCRs--theoretical and experimental studies. Curr Med Chem. 2012;19:1090–109.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  63. Teichmann SA, Chothia C. Immunoglobulin superfamily proteins in Caenorhabditis elegans. J Mol Biol. 2000;296:1367–83.

    Article  CAS  PubMed  Google Scholar 

  64. Lin M, Sutherland DR, Horsfall W, Totty N, Yeo E, Nayar R, et al. Cell surface antigen CD109 is a novel member of the alpha(2) macroglobulin/C3, C4, C5 family of thioester-containing proteins. Blood. 2002;99:1683–91.

    Article  CAS  PubMed  Google Scholar 

  65. Kondo M, Scherer DC, Miyamoto T, King AG, Akashi K, Sugamura K, et al. Cell-fate conversion of lymphoid-committed progenitors by instructive actions of cytokines. Nature. 2000;407:383–6.

    Article  CAS  PubMed  Google Scholar 

  66. Reymond N, Fabre S, Lecocq E, Adelaide J, Dubreuil P, Lopez M. Nectin4/PRR4, a new Afadin-associated member of the Nectin family that trans-interacts with Nectin1/PRR1 through V domain interaction. J Biol Chem. 2001;276:43205–15.

    Article  CAS  PubMed  Google Scholar 

  67. Sheridan C, Sadaria M, Bhat-Nakshatri P, Goulet R, Edenberg HJ, McCarthy BP, et al. Negative regulation of MHC class II gene expression by CXCR4. Exp Hematol. 2006;34:1085–92.

    Article  CAS  PubMed  Google Scholar 

  68. Contento RL, Molon B, Boularan C, Pozzan T, Manes S, Marullo S, et al. CXCR4-CCR5: a couple modulating T cell functions. Proc Natl Acad Sci. 2008;105:10101–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  69. Zlotnik A, Yoshie O. Chemokines: a new classification system and their role in immunity. Immunity. 2000;12:121–7.

    Article  CAS  PubMed  Google Scholar 

  70. Shanaka WW, Rodrigo I, Jin X, Blackley SD, Rose RC, Schlesinger JJ. Differential enhancement of dengue virus immune complex infectivity mediated by signaling-competent and signaling-incompetent human Fc-gamma-RIA (CD64) or Fc-gamma-RIIA (CD32). J Virol. 2006;80:10128–38.

    Article  Google Scholar 

  71. Williams MJ, Ando I, Hultmark D. Drosophila melanogaster Rac2 is necessary for a proper cellular immune response. Genes Cells. 2005;10:813–23.

    Article  CAS  PubMed  Google Scholar 

  72. Zanin-zhorov A, Cahalon L, Tal G, Margalit R, Lider O, Cohen IR. Heat shock protein 60 enhances CD4 + CD25 + regulatory T cell function via innate TLR2 signaling. J Clin Invest. 2006;116:2022–32.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  73. Vasilevko V, Ghochikyan A, Holterman MJ, Agadjanyan MG. T-cell proliferation after activation with suboptimal doses of PHA. DNA Cell Biol. 2002;21:137–49.

    Article  CAS  PubMed  Google Scholar 

Download references

Acknowledgements

To Fapesp, (Sao Paulo Research Foundation, grants # 2010/05516-7, #2011/00396-6 and #2014/07566-2). MP Berton, E Peripolli received scholarship from Coordination for the Improvement of Higher Education Personnel (CAPES; Coordenação de Aperfeiçoamento de Pessoal de Nível Superior) in conjunction with the Postgraduate Program on Genetics and Animal Breeding, Universidade Estadual Paulista, Faculdade de Ciências Agrárias e Veterinárias (FCAV, Unesp). F. B held productivity research fellowships from The Brazilian National Council for Scientific and Technological Development (CNPQ)

Availability of data and materials

The phenotypic and genomic information used in this study belongs to a private sheep breeding program company, so we do not have authorization to share the data.

Funding

Sao Paulo Research Foundation – FAPESP grants # 2010/05516–7, #2011/00396–6 and #2014/07566–2.

Authors’ contributions

MPB, RMOS, EP, NBS did the manuscript writing and discussion of the results. JFM, MSA, BVG, in Genomic modeling and genomic analysis. GB in concepts, writing and modeling. PSO in phenotypic data collection and genotyping. JPE in data management. FB JBSF in conception, funding, modeling, genomic analysis and coordination. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Consent for publication

Not applicable.

Ethics approval

All experimental procedures involving animals were approved by the Ethical Committee of FZEA/USP CEUA n°7,718,021,216.

Author information

Authors and Affiliations

Authors

Corresponding author

Correspondence to Mariana Piatto Berton.

Additional files

Additional file 1: Table S4.

Genomic regions associated with egg counts per gram of feces (EPGlog) in Santa Inês sheep; Table S5. Genomic regions associated with famacha (FAM) in Santa Inês sheep; Table S6. Genomic regions associated with blood cell count (WBC) in Santa Inês sheep; Table S7. Genomic regions associated with red blood cells (RBC) in Santa Inês sheep; Table S8. Genomic regions associated with platelet (PLT) in Santa Inês sheep; Table S9. Genomic regions associated with hematocrit (HCT) in Santa Inês sheep; Table S10. Genomic regions associated with blood cell count (HGB) in Santa Inês sheep. (XLS 74 kb)

Additional file 2: Table S11.

KEGG pathways (P < 0.05) for counting eggs per gram of feces (EPGlog), famacha (FAM), white blood cells (WBC), red blood cells (RBC), platelets (PLT), hematocrit (HCT), and hemoglobin (HGB) traits obtained from the DAVID software; Table S12. Gene Ontology terms for biological processes significantly (P < 0.01) related with counting eggs per gram of feces (EPGlog), famacha (FAM), white blood cells (WBC), red blood cells (RBC), platelets (PLT), and hematocrit (HCT) traits. (XLS 68 kb)

Rights and permissions

Open Access 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.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Berton, M.P., de Oliveira Silva, R.M., Peripolli, E. et al. Genomic regions and pathways associated with gastrointestinal parasites resistance in Santa Inês breed adapted to tropical climate. J Animal Sci Biotechnol 8, 73 (2017). https://doi.org/10.1186/s40104-017-0190-4

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40104-017-0190-4

Keywords