Unraveling the genetic variability of host resilience to endo and ectoparasites under natural infestation, in Nellore cattle

Background: 25 Host resilience (HR) to parasites can affect growth in pastured raised cattle. This study is a detailed 26 investigation of the genetic mechanisms of HR to ticks (TICK), gastrointestinal nematodes (GIN), and 27 Eimeria spp. (EIM) under natural infestation. HR was defined as the slope coefficient of random 28 regression models of body weight (BW) when TICK, GIN, and EIM burdens were used as 29 environmental gradients. The BW was evaluated in five measurement events (ME): when animals were 30 331, 385, 443, 498, and 555 days old on average. 7307 BW records were available from 1712 animals 31 weighted at least in one ME. Out of those, 1075 animals had valid genotypic information after quality 32 control analysis that were used in genome-wide association studies (GWAS) and GWAS meta-analyses 33 to identify genomic regions associated with HR. Both the genetic correlations between intercept and HR to each parasite, and the genetic correlations between BW measured in animals submitted to different parasite burden indicated that there was genotype x parasite burden interaction for BW, and selection for BW under environment with 38 controlled parasite burden might be an efficient strategy to improve both, BW and HR. Furthermore, there was no impact of age of measurement on genetic variance estimates for HR to different parasites. However, genetic correlation between HR to the same parasite measured in different ages ranged from 41 low to moderate in magnitude, with a posteriori means (high posterior density interval with 90% of 42 samples) varying from 0.13 (-0.05; 0.35) to 0.40 (0.15; 0.63) for TICK, from 0.11 (-0.06; 0.29) to 0.52 43 (0.37; 0.67) for GIN and from 0.25 (0.07; 0.43) to 0.56 (0.34; 0.77) for EIM. These results indicate the 44 importance of age of measurement in studies on HR. and a 47 multiple trait selection method that combine BW and HR to parasites should be used in parasitic 48 endemic areas to avoid economic losses due parasitic diseases.


1
Background 53 Ecto and endoparasites such as ticks (TICK), gastrointestinal nematodes (GIN), and Eimeria spp. 54 (EIM) are endemic in tropical countries and they are also responsible for several economic and 55 productivity losses in cattle production systems [1]. Moreover, parasitic loads represent an important 56 challenge to the sustainability of cattle production, especially in tropical countries, such as Brazil. 57 The negative impact of ticks on cattle production is due to the direct effects of feeding, such as weight 58 loss, anaemia, and damage of leather, and indirect effects, such as the transmission of tick-borne 59 pathogens [2]. Gastrointestinal parasites like GIN and EIM also negatively affect the cattle 60 performance due to both direct and indirect effects: competition for nutrients, physical tissues damage, 61 and the host immune response to parasite invasion [3][4][5]. Reduced food intake, weight loss, diarrhea, 62 and dehydration are the main symptoms of these intestinal parasitosis [4,6,7]. For EIM infections, 63 anaemia is also an important symptom [3]. 64 The sustainability of cattle production in endemic areas depend on the animal's ability of respond to 65 stressor factor as parasite loads. The terminology around the different mechanisms of host's response 66 to disease is still confuse. In general, animals can respond to disease using two complementary 67 mechanisms: resistance and tolerance [8]. The host resistance can be characterized by the ability of a 68 host to limit parasite burdens while host tolerance can be defined as the ability to limit the damage 69 caused by a given parasite burden (Råberg et al., 2009). Therefore, the host resilience (HR) is the 70 phenotype that captures these two mechanisms against pathogens [8,9], and can be defined as the 71 animal capability of maintain a relatively undepressed production level when subjected to 72 environmental parasite burdens [10,11]. 73 The HR can be estimated as a continuous trait using reaction norm models of performance on 74 environmental parasitic load [8]. In this case, the additive variance of dependent variable, is divided 75 into two coefficients: the intercept that is the additive component of the variability in performance and 76 the slope that is the HR [12]. Moreover, when linear regressions are used, the genetic correlation 77 between the intercept and slope coefficients quantify the genetic association between host fitness and 78 HR [13]. Significative correlations between these two parameters indicate, thus, the presence of 79 genotype x environmental interaction [14]. 80 Reaction norm models have been used to estimate HR in milk production and fertility traits to Fasciola 81 hepatica in Irish cattle [15]. The authors used the prevalence of Fasciola hepatica in each herd and 82 year to define the environmental gradient. In our study, we estimated HR in body weight to TICK, 83 GIN, and EIM burdens to estimate genetic parameters of HR in beef cattle. The slope solutions were 84 considered as the estimated breeding values for HR and used as a phenotype to perform genome-wide 85 association studies (GWAS) for HR to each parasite. Therefore, the main objectives of the present 86 study are to estimate genetic parameters for HR of cattle to different parasites, evaluate the trends of 87 these traits with aging, and suggest possible mechanisms influencing HR to parasites in Nellore cattle. 88 In short, this study is an investigation of the genetic mechanisms of HR to endo and ectoparasites under 89 natural infestation. The genetic parameters of HR and its association with BW uncovered here are 90 evidence for discussing the inclusion of HR as selection criteria in cattle breeding programs. 91 The bulls were weighed at six measurement events (ME): at day 1 of the performance test, at the end 110 of the adaptation period (day 70) and in intervals of 56 days until the end of the test. BW information 111 registered at day 1 was not used in this study. The intervals of 56 days defined five ME. The animals 112 were on average 331, 385, 443, 498, or 555 days old in each of the five ME. 113 The parasite counts used in the present study were obtained at each ME through counts of engorged 114 female ticks (length size > 4.5 mm -TICK) on the right side of each animal. The length size of the 115 engorged female ticks were defined according to the technique proposed by Wharton and Utech [16]. 116 Furthermore, fecal samples were collected directly from the animals' rectum using properly identified Sheater's solution (500 g of sugar, 6.5 ml of phenol and 360 ml of water). Then, 0.15ml of the final 124 solution was used to fill the McMaster chamber used to perform the counts for eggs and oocysts. 125 There is a practice of multiplying tick counts by two to estimate the tick burden in the entire animal 126 [18]. For the modified McMaster technique, a multiplicative factor can also be used to infer on the 127 animal's parasite burden. This multiplicative factor depends on the dilution used to prepare the test. 128 For the dilution we used, the adequate multiplicative factor is 50 [19,20]. Ticks, eggs, and oocysts 129 counts used in the present study are the real counts observed on the right side of each animal or at the 130 McMaster chamber, without multiplication by any constant. The bulls included in the present study 131 were subjected to natural parasite infestation because they were raised in a herd with commercial 132 purposes. Approximately 65% of the bulls were dewormed at the beginning of the performance tests 133 (day 1 of adaptation period) with Ivermectin 4% (1ml of Ivermectin per 50Kg of live BW -Master LP, 134 Ouro Fino Saúde Animal, Cravinhos, SP). The dewormed bulls were randomly chosen by 135 contemporary group, in such a way that we had some entire groups dewormed or not. The 136 contemporary groups were defined as the group of animals that were raised together in the same 137 paddock. 138 5 Blood samples were collected with sterilized syringes into vacuum tubes of 3.5 ml containing 9NC 139 Coagulation Sodium Citrate 3.2%, to conserve the host's DNA . Blood samples were frozen and  140  transferred in chilled coolers to the Laboratory of Genetics at EV-UFMG, where they were stored in  141  freezers at -20ºC. 1230 blood samples were selected for genotyping with a low-density DNA array: the  142  Z-chip v2 (Neogen, Lincoln, Nebraska, EUA, which genotypes -27533 SNPs). Most of the genotyped  143  bulls were from the performance tests with more than 20 animals per group, as described above, and  144 each animal had information regarding the three parasites, in at least four ME. Some genotyped animals 145 did not have phenotypic measures, but they were representative sires of this herd (31 genotyped sires  146 without phenotypic records with an average of 25.26 offspring in the relationship matrix). 147

Phenotypic data editing 148
The data set of phenotypes had information on BW, TICK, GIN, and EIM. For each ME we have only 149 considered animals that had information regarding the four phenotypes. Cohorts were defined by the 150 combination of contemporary group and ME. Bulls belonging to cohorts with less than five animals 151 were not considered in the study. At the end of the data editing process, 1712 animals had information 152 in at least one ME (Table S1) The environmental pathogen load is a crucial component of resilience that must be considered [8]. 167 Thus, we used the median counts of TICK, GIN, and EIM in each cohort as the environmental gradient, 168 that is formed by the combination of animals raised in the same paddock (contemporary groups) on the 169 same period of the year (ME). Each combination of parasite x ME was used to generate a different 170 dataset, for which the genetic parameters were estimated. Cohorts with median parasite count equal to 171 zero do not indicate they are parasite free cohorts because at least one animal in each cohort had parasite 172 counts greater than zero. In short, every animal in our dataset were exposed to natural infestation of 173 the three different evaluated parasites. 174

Genetic parameters for body weight and host resilience 175
We used single trait linear random regression models (STM) where BW in each ME was considered 176 as a dependent variable (trait) and the median counts of parasites as an independent variable. It is 177 important to emphasize that in this study the random slope coefficient was considered as the HR to 178 parasite. The analysis was performed for each ME and each parasite separately, and no effect of co-179 infection was considered. The STM was performed using Bayesian inference methodology, and can be 180 described as: 181 where ijkl y , represents the weight of the animal i, evaluated at the cohort j, with the age k and submitted 183 to an environmental parasite burden l; j C , is the systematic effect of the cohort; 1 d , is the slope to fit 184 the effect of age in which each animal was evaluated; ( The breeding values (EBV) estimated for HR to each parasite in each ME using STM were considered 217 as the phenotypes of HR for GWAS analysis. No fixed effect was considered since all the known 218 environmental effects were considered in the previous analysis to estimate the breeding values for HR. 219 The model used for the GWAS analysis for HR can be described as: 220 Where HR is the vector of EBVs for host resilience to TICK, GIN, or EIM in each ME and the other 222 terms as previously described. 223 In the GWAS mixed models, for BW and HR, we used a genomic relationship matrix (GRM) [ and a full dosage compensation correction was applied to include the X chromosome markers in the 227 estimate for the GRM. This software uses the algorithm proposed by Taylor [23]. The effect of sex 228 over the solutions for SNPs located in the X chromosome was also considered in this GWAS analysis. 229 The heritability for BW and HR were estimated based on the variances obtained from the markers 230 (using GRM), so we refer to that estimate as "SNP-derived heritability". It is important to highlight 231 that, in the present study, the heritability for BW was estimated based in both GRM and the pedigree-232 based relationship matrix, so for BW we presented both the conventional heritability and the SNP-233 derived heritability estimates. 234 The correlation matrix among BW and HR to TICK, GIN, and EIM measured in each ME was 235 computed by the Pearson correlations between the SNP effects (SNP correlations) for each one of the 236 traits, as proposed by Fortes et al. [24]. To compute the SNP correlation matrix, we first standardized 237 the SNP effects estimates by its standard error. 238

Genetic variances for host resilience across ages 239
We used two-trait linear random regression models (TTM) to model BW in each ME in function of an 240 intercept and a slope defined by the median counts of parasites. The ME information (ME.331, 241 ME.385, ME.443, ME.498, and ME.555) were analyzed by two-trait analysis, to achieve convergence 242 of the parameters of the model. The TTM can be described by the same equation of STM. For this 243 model, the covariance matrix for the random effects (G1) was assumed as: 244 the P-values, effect direction, and sample size of each GWAS, a Z-score and an overall P-value for 262 each marker were calculated. The meta-analysis was performed for markers that had valid solutions in 263 at least 2 studies and no genomic control was performed during the meta-analyses. Afterwards, we 264 used a Bonferroni correction for multiple testing to define the threshold for SNP effect significance. 265 Only SNPs with P-values < 2.31x10 -6 were considered as significant. Another threshold of P-values < 266 10 -4 was used to infer suggestive SNPs, a common practice in GWAS [26,27]. 267 Based on meta-analysis results, quantitative trait locus (QTL) associated with each trait were described. 268 The QTL boundaries were defined as follow: First, for each bovine chromosomes (CHR), an initial 269 peak SNP was defined as the SNP with the lowest P-value; Second, around the peak SNP, regions of 270 0.5Mbp up and downstream were searched for other significant SNPs. If inside this interval we 271 identified other significative SNP, the boundaries of the QTL were expanded to include the SNP and 272 another 0.5Mbp (up and downstream) was investigated. The process was repeated until there was no 273 more significative SNPs in these 0.5Mbp windows. Finally, a new peak SNP was called if there was a 274 significative SNP in the same CHR but outside of the boundaries of the first QTL. The process was 275 repeated within each CHR until no more peak SNPs could be identified. 276 To provide additional evidence to QTLs associated to HT, we imposed one more criterion. Only 277 regions with at least four significative or suggestive SNPs were considered as a QTL (adapted from 278 van den Berg et al., 2016). Suggestive SNPs had a P-value <10 -4 and to be considered as supporting 279 evidence for the QTL they had to satisfy one more condition: to have above average LD with the peak 280 or other significant SNP in the QTL. The LD between SNPs was evaluated by the D prime (D') 281 estimated using the expectation-maximization method at pairwise analysis carried with the SNP & 282 Variation Suite v8.8.3 [21]. SNPs were considered in high LD when D' was greater than the mean + 2 283 standard deviations of the D' computed between all combinations of markers for each CHR. For the 284 traits in which no QTLs could be described, we identified the genes around isolated significant SNPs 285 (0.5Mbp downstream and upstream the SNP position). 286 We looked for genes located inside the QTL boundaries using the ARS-UCD1.2 bovine genome 287 assembly (available at www.ncbi.nlm.nih.gov/assembly/GCA_002263795.2). The search for genes 288 was made using the GALLO package [29] of R software [30]. Genes located inside QTLs or around 289 isolated significant SNPs were considered as candidate genes. Candidate genes formed target gene lists 290 that were confronted with a trained gene list in functional analysis for candidate gene prioritization, as 291 described below. 292 The trained list of genes was constructed using keywords (Table S2) that described each of evaluated  293 phenotypes (BW and HR to the three parasites). These lists were built on GUILDify v2.0 [31], which 294 9 is a web application for phenotypic characterization of genes. GUILDify searches for genes starting 295 from user-provided keywords in the Biologic Interaction and Network Analysis (BIANA) knowledge 296 database. These genes associated with the keywords are used as seeds to generate the protein interaction 297 networks, for the selected organism, analyzed with graph theory algorithms to prioritize new disease 298 genes [31]. In the present study, the selected model organism was Homo sapiens, since bovine was not 299 an option. The Netscore prioritization algorithm from the GUILD package was used (with repetition = 300 3 and interaction = 2; default values of GUILDify). The output of GUILDify is a trained list of genes, 301 ranked according to the interaction network. The first 100 genes were used as the trained gene list for 302 each studied trait. 303 Candidate gene prioritization analysis were conducted with ToppGene Suite [32]. These analyses were 304 performed in two-steps. First, for each trait, a functional enrichment analysis was performed to verify 305 if the trained gene list was enriched for any functional category or parameter. We used Gene Ontology 306 (Molecular function, Biological process, and Cellular component), Human phenotype, Mouse 307 phenotype, Pathway, PubMed, Transcription factor binding site, Co-expression, and Disease as 308 training databases to identify over-representative terms from the trained gene list. The P-value cut-off 309 for each training parameter was 0.05 with a False Discovery Rate correction. After this step a 310 representative profile of the trained gene list was obtained. 311 In the second step a similarity score was generated for each gene in our candidate gene lists. This score 312 is created by functional annotation of the candidate gene followed by a comparison of its function to 313 each enriched term, learned in the training step. The similarity score calculation and the P-values 314 associated to them are described in Chen et al. [32]. In summary, a fuzzy-based similarity measure is 315 applied for categorical terms [33], and Pearson correlation between the test gene and the enriched gene 316 lists is applied for quantitative functional parameters. In the case of a missing value (for instance, lack 317 of one or more annotations for a test gene), the score is set to −1. Otherwise, it is a real value in [0, 1] 318 [32]. 319 The overall scores and P-values are obtained with a meta-analyses that considers all the functional 320 categories annotated [32]. The prioritized genes were considered those with an overall P-value ≤ 0.05. 321 For the candidate gene prioritization analysis, we used the default setting in ToppGene Suite that is a 322 background gene set from the genome for computing the P-value with 5000 coding genes and two 323 features to be considered for prioritization. 324 3 Results 325

Environmental parasite burden 326
The environmental parasite loads observed here were low (Table 1). It might be a consequence of 327 preventive treatment applied to 65% of animals randomly chosen on the first day of each performance 328 test. A common strategy to estimate the total count of parasite per animal is to multiply by 2 the number 329 of engorged females observed in the right side of each animals and it is important to highlight that the 330 parasite counts reported here were not multiplied by any constant. For instance, ticks' count represents 331 the number of ticks observed on the right side of the animals, and the eggs and oocysts' counts represent 332 the number of parasites observed on the McMaster chamber (Table 1). 333

Genetic parameters for body weight and genotype x parasite burden interaction 334
In general, there is no significant difference between genetic parameters of intercept and slope 335 coefficients estimated by STM (Table 2) and TTM (Table S3). Except for genetic additive variances of intercept coefficient at ME.555 when EIM parasite burden was considered in the model.  (Table S3). Therefore, the results of STM are presented and discussed 340 in the main text (Table 2, Figure 2) and the results for TTM are presented on the supplementary Table  341 S3 and Supplementary Figure 2. 342 The SNP-derived heritability (average ± standard error) of BW (Table 3)   for GIN and 10.5 for EIM), the HPD90 related to those posterior means were large showing no 358 significant differences between them (Figure 2). 359 The SNP correlations between the BW measured in each ME are high, with averages (standard errors) 360 ranging from 0.705 (0.005) between BW at ME.331 and BW at ME.555 to 0.882 (0.003) between BW 361 at ME.498 and BW at ME.555 ( Figure 3). These high SNP correlations between BW measured from 362 331 to 555 days old are in accordance with GWAS results for BW in which markers located at 363 chromosome 6 and 14 were suggestively associated with BW measured at all evaluated ME 364 (Supplementary Figure 3). Also, BW measured at different ages could be considered repeated measures 365 of the same phenotype, instead of being perceived as independent traits since animals that are heavier 366 in the beginning of performance tests tend to be heavier in the end. 367 The genetic correlations of BW measured at each ME between different parasite burden varied. They 368 ranged from high and positive, between animals submitted to similar parasite burden, to moderate and 369 negative, between animals submitted to extreme different parasite burden (Figures 4, 5, 6 It is possible to verify that negative correlations between BW measured at zero and maximum counts 377 occurred at the ME in which negative and significant covariances were estimated between intercept 378 and slope ( Table 2). The genetic correlations between intercept and slope ( verified. For breeding purposes, the identification of this interaction might help to define the 384 environmental conditions and management practices to which candidate animals will be submitted. 385

Genetic parameters for host resilience and their association with body weight 386
In this study, HR was estimated as a genetic component of BW, the slope of random regression models. 387 Therefore, genetic variance estimates were obtained for HR, but heritability was not estimated. There 388 was no impact of age of measurement on HR genetic variances to TICK, GIN, or EIM regardless of 389 the model used (STM and TTM) since HPD90 of genetic variance for HR across ME overlapped (Table  390 2, Table S3). Also, when we compare STM and TTM results, for the same ME, the HPD90 of HR 391 genetic variance to TICK, GIN or EIM overlapped again. These overlaps are evidence for the fact that 392 the HR's genetic variances were similar across age groups and statistical models. 393 The SNP-derived heritabilities for HR to TICK, GIN, and EIM at each ME were computed through 394 GWAS analyses when the slopes solutions were considered as HR phenotype. As expected, these 395 estimates presented high magnitude ( The genetic association between BW and HR was measured throughout SNP correlations between the 398 traits ( Figure 3). In our study, both positive and zero correlations were considered as favourable. 399 Negative and unfavourable correlations were observed between BW and HR to TICK at ME.331 (-400 0.648±0.005), ME.443 (-0.307±0.006), and ME.498 (-0.148±0.007), between BW and HR to GIN at 401 ME.385 (-0.038±0.007), and between BW and HR to EIM at ME.555 (-0.081±0.007 - Figure 3). 402

Host resilience to parasitic burden at five ages are independent traits 403
The posterior means of genetic variances for HR to TICK, GIN, and EIM were similar in each ME 404 ( with GWAS results where suggestively associated markers for TICK, GIN, or EIM were age-specific 418 ( Therefore, the trend of HR to parasites changes at different ages, for instance HR should be considered 420 as an independent trait. 421 3.5 Candidate genes and pathways associated with host resilience to TICK, GIN, and EIM 422 The search for genes associated with HR was carried only at the regions defined as QTL or in the 423 vicinity of significant SNP, considering only the meta-analysis results, and not the results for each age 424 separately (Figure 8). The aim was to identify candidate genes that might explain the genetic 425 correlations observed between HR to each parasite across the different ages. In short, we focused on 426 the associations (and therefore candidate genes) that had stronger evidence by combining in one 427 analysis all the HR data available for each parasite. 428 There were no relevant QTLs identified for HR to TICK through the meta-analysis. However, we 429 identified candidate genes nearby significant SNPs. A total of 52 genes formed the candidate list for 430 HR to TICK: 11 on CHR 2, 6 on CHR16, and 35 on CHR19. (Table 4). Out of these candidates, 21 431 were prioritized in our candidate gene prioritization analyses. Information about genes prioritized for 432 HR to TICK is presented at Table S4. 433 There were significant SNPs associated to HR to GIN on CHR 9, CHR 14 and CHR 28 ( Brazil´s total milk yield production). They concluded that only 195 farms (22%) do not use any 468 prophylactic drug to control parasites load, including Eimeria spp. and gastrointestinal nematodes. 469 Furthermore, more than one-half of beef cattle farmers in USA usually dewormed different categories 470 of animals one or more times per year [38]. Therefore, environment with controlled parasite burden 471 through prophylactic treatment represents the reality of commercial farms. 472 Even though under low parasite load challenge we could observe the impact of parasite burden on the 473 body weight, and a raising trend for the heritability of body weight as the loads increased. It is expected 474 that more challenging environments, this means, higher parasite loads, can lead to more significant 475 effects parasite burden on both BW and genetic parameters estimates for HR (Falconer, 1990). 476 Considering this, the low parasite loads observed on the present study might be a limiting aspect of our 477 study for the study of the genetic architecture of HR to different parasites, and further genomic regions 478 associated with these traits could be found when applying the same methodology we used on a 479 population submitted to higher burdens. to other studies that might characterize a commercial farm environment and should be used to evaluate 495 the genetic parameters of health traits in beef cattle industry. 496

Genetic parameters for body weight and genotype x parasite burden interaction 497
The similarity between BW heritabilities estimated with genomics (SNP-derived heritability), STM, 498 and TTM serve as evidence for the adequacy of the low-density SNP Mundo Novo might explain the lower heritability estimates obtained here. It is important to highlight 507 that those two genetic management practices, selection and non-inclusion of external candidates, 508 promoted significant benefits to this population that is highly regarded as a genetic disseminator of the 509 Nellore breed in Brazil. 510 In the present study we did not find significant changes in BW heritability estimates across varying 511 gradients of parasite burden. The natural infestations and overall low parasite burden observed in this 512 data might partially explain the constant heritability, because challenging environments might lead to 513 increases in heritability estimates [45,46]. For instance, heritability for FAMACHA score in ram and 514 ewe lambs submitted to high worms burden, mostly Haemonchus contortus, Trichostrongylus spp., 515 and Teladorsagia spp., were significantly higher than the heritability obtained in low and moderate 516 parasite burden (Riley and Van Wyk, 2009). FAMACHA score is a common veterinarian approach to 517 evaluate the parasite burden, based on the colour of the animals' eye mucosa (indicative of anaemia). 518 Moreover, an uprise trend in heritability estimates for milk yield was observed with increased 519 temperature-humidity index, a direct indicator of heat stress (Lee et al., 2011). Therefore, we expect 520 that BW evaluation on infested environments, with higher range parasite counts than those observed 521 here, might lead to variation in heritability estimates of BW across varying gradients of parasite burden. 522 We were able to verify a rising trend on heritability means (mainly for HR to TICK and GIN at ME.555, 523 and HR to EIM at all ME) however these values were followed by large HPD. It is important to 524 highlight that Mundo Novo farm has an efficient animal husbandry program, including strategic 525 parasite control, that correspond to the outstanding practices of a nucleus farm and it explains the low 526 parasite burden observed in our study, including the more uniform weight gain. 527 Our results showed low genetic correlations between BW measured at the lowest and the highest 528 parasite loads specially at younger ages (331 days old compared to 555 days old), indicating potential 529 genotype by parasite burden interaction when animals are raised in very low parasite burden (cohort's 530 median count equal zero) and infested environments. Hollema et al.
[47] also emphasize the importance 531 of considering the genotype x worm burden interaction for growth rate in Australian Merino sheep to 532 increase the efficiency of selection for animals that are more parasite resistant and more resilient to 533 environmental worm challenge. However, these authors verified significant decrease on the heritability 534 for growth rate with the increase of worm burden [47]. It becomes even more important for growth 535 selection on pastured systems in tropical areas that are typically under different levels of natural 536 infestation, and that apply different strategies for parasite control. 537

Genetic parameters for host resilience and their association with body weight 538
The SNP-derived heritability estimated for HS in the present study are high, but they do not indicate 539 that HS is highly heritable. In fact, these values are a statistical artefact since the phenotype of HS is 540 itself an estimated breeding value. However, the SNP-derived heritability indicate that breeding values 541 estimated using pedigree-based genetic evaluations can be efficiently explained by the genomic 542 similarity between individuals. There is genetic variance for HR to parasites and genetic improvement 543 for this trait can be achieved through selection, even though the genetic gain across generation might 544 be slow. The SNP correlations between HR and BW were obtained through standardized values, where 545 the SNP effects were divided by its standard error. Standardized values reduced the impact of large 546 differences on genetic variance of BW and HR obtained here when computing SNP correlations 547 estimates. Furthermore, the SNP correlations were particularly interesting since the effect of parasite 548 infestation was not considered in the BW's GWAS analyses. In this sense, while genetic correlation 549 between intercept and HR indicates the presence of genotype x parasite burden interaction for BW, the 550 SNP correlations indicates some genetic association between BW and HR to parasites. 551 Unfavourable correlations were observed mainly between BW and HR to TICK at ME.331, ME.443, 552 and ME.498. These correlations were moderate and negative at ME.331 but were reducing in 553 magnitude as animals aged. Two aspects must be discussed here. First, the fact that inside each age 554 category (ME), heavier animals are expected to be larger in size, and have a wider skin surface with a 555 denser vasculature, which was already suggested as a possible explanation for the associations between 556 parasite burden and BW [40], and might also explain the association between BW and HR. Second, 557 the immune response mechanisms might vary with age [48-51] and so younger animals (for instance 558 those evaluated in ME.331) might differ from older ones (like those evaluated in ME.555) in terms of 559 HR. Altogether, these two aspects can justify the differences observed in the SNP correlations between 560 HR to TICK and BW in the different ME. Further discussion about the age effect on HR expression is 561 presented below. 562

Host resilience to parasitic burden at five ages are independent traits 563
The low to moderate magnitude of genetic correlation between slopes at each ME in TTM and SNP 564 correlations for HR to parasites in each ME indicate that HR to TICK, GIN or EIM measured at 565 different ages are independent traits, showing no benefit of using indirect selection to improve HR in 566 older animals through selection of younger ones and vice and versa. It is already known that the severity 567 of infectious diseases can vary dramatically across ages, possibly due to the immaturity of the immune 568 system of young animals [52]. Similarly, it is possible that host resilience also changes across ages, 569 depending on the exposure to co-infestation and different parasites in each growing phase [53]. 570

Candidate genes and pathways associated with host resilience to TICK 572
We found here a significative association between the genes WNT4 and CFH and HR to TICK. The 573 WNT4 gene was recently associated with the suppression of type 2 immunity (Hung et  infestations. 582 Another homeostasis-related mechanism that appears to be relevant to the expression of HR to TICK 583 is apoptosis, suggested by the association of YWHAE and SCARF1 genes. Apoptosis is the 584 mechanisms that ensures that damaged, aged, or excess cells are deleted in a regulated manner that is 585 not harmful to the host and plays an essential role in the development and maintenance of all 586 mammalian tissues [59]. The insertion of the tick's mouthparts during tick blood feeding leads to 587 damages of epidermis and dermis cells [60]. In response to damage, pathogen-associated and or 588 damage-associated molecular patterns are detected, and leukocytes aggregate near to the site of lesion 589 [61]. After eliminating the initial threat, leukocyte recruitment ceases, and the previously recruited 590 cells are disposed, mainly mediated by apoptotic cells that are subsequently phagocyted [61]. A rapid 591 and immunologically 'clean' removal of apoptotic cells by neighbouring phagocytic cells is essential 592 for the maintenance of homeostasis and avoidance of inflammation [62]. This mechanism may explain 593 why genes that are relevant to apoptosis, like YWHAE, and SCARF1, may also be important for HR 594 to TICK. 595 16 YWHAE is a protein isoform of 14-3-3 family of eukaryotic proteins which have anti-apoptotic 596 activity, and regulate members of the mitochondrial apoptotic machinery, as well as a staggering 597 number of signalling molecules that mediate the transmission of survival and death signals to the 598 mitochondrial death machinery [63]. SCARF1 is one of the main genes mediating the clearance of 599 apoptotic cells [64]. These mechanisms are part of a delicate balance in order to maintain the 600 homeostasis. To understand how they take part on the expression of HR further studies are necessary. 601 It would be useful, for example, to evaluate if these genes are up or downregulated in more resilient 602 hosts upon infestation. 603 Alpha-2 pigment epithelium derived factors secretion (that are transcripts of SERPINF1 and 604 SERPINF2 genes) were related to HR to TICK. SERPINF1 was associated with mice full-thickness 605 cutaneous wound healing by promoting epithelial basal cell and hair follicle stem cell proliferation. 606 SERPINF2 acts on other aspects of the wound healing, in the clearance of the anticoagulant plasma 607 protein C, inhibiting its actions and enhancing the coagulation process [65,66]. Altogether, the 608 mechanisms controlled by SERPINF1 and SERPINF2 might indicate the relevance of fast wound 609 healing on the expression of HR to TICK. 610 The direct effect of tick saliva might also activate specific mechanisms of HR. Salp15 is the first protein 611 associated with the immunosuppressive activity of Ixodes scapularis tick saliva (Anguita et al., 2002). 612 CDC42 activation can act as a stimulus to F-actin polymerization, and the amount of F-actin was 613 reduced upon pretreatment of T cells with Salp15 in mice [67]. The complexity of events at the tick 614 host interface is increased by the process in which injection of saliva occurs alternatingly with uptake 615 of blood as well as of digested tissues at an increasing rate over the course of blood feeding [58]. It 616 was demonstrated that Salp15 inhibited the activation of CD4+ T-cells in mammalian hosts, resulting 617 in decreased activation of the transcription factor NF-κB [68]. nematode parasites are endemic, immunity to infection in previously exposed individuals is associated 623 with expression of T-helper type-2 (TH2) cytokines and an inverse association between TH2 cytokines 624 and susceptibility to GIN in younger children (between 4 and 13 years old) was verified [70]. 625 AGO2 was associated in the present study with HR to GIN  The use of AGO2-miRNA pathway in the regulation of immune response of cattle infected by GIN 636 was not described yet, however previous results presented here as well as the significant association of 637 AGO2 with HR to GIN might indicate the relevance of this pathway. 638 17 The restitution of intestinal epithelial barrier damage, that might be caused by GIN infections, is 639 another mechanism that appears to be relevant for the expression of HR. CXCL12 activates the 640 chemokine receptor CXCR4 and enhances intestinal epithelial wound healing through reorganization 641 of the actin cytoskeleton [78]. CXCL12 is a constitutive and inflammatory chemokine in the intestinal 642 immune system, expressed by normal intestinal epithelial cells [79]. It was up-regulated in both catfish 643 intestine submitted to following experimental infection of Edwardsiella ictalurid, the causative 644 bacterium of enteric septicemia of catfish [80], and human intestine of patients with inflammatory 645 bowel disease [79]. No associations of CXCL12 expression in the intestine of animals submitted to 646 gastrointestinal nematodes burden was published yet. However, our results indicate that this might be 647 an important candidate gene for the study of HR to GIN. 648

Candidate genes and pathways associated with host resilience to EIM 649
The chemokines also play an important role on the development of HR to EIM, as the association of 650 genes CXCL9, CXCL10, and CXCL11 suggests. Tables 712 Table 1       (ME). "Age" represents the mean age that animals had at each ME. "nb" is the number of bulls 762 evaluated at each ME and "nc" is the number of cohorts evaluated at each ME. Red arrow indicates a 763 70-day interval between evaluations, while blue arrows indicate a 56-day interval. 764   (TICK), evaluated in five different evaluation periods (ME.331, ME.385, ME.443, ME.498, ME.555). 774 331, 385, 443, 498, and 555 represent the mean ages in days that the animals had in each evaluation. 775 The x and y-axis of each plot correspond to the ticks' burden observed in each period (only the 776 minimum, maximum and the three quantiles of parasite burden are presented). 777 between host resilience to ticks (TICK), gastrointestinal nematodes (GIN) and Eimeria spp. (EIM) at 790 different measurement events (ME) of Nellore bulls. ME.331, ME.385, ME.443, ME.498, ME.555 are 791 evaluation periods when animals' age were 331, 385, 443, 498, and 555 days old in average. 792