Genome-Wide Association Study for Haemonchus contortus Resistance in Morada Nova Sheep

Among the gastrointestinal nematodes affecting sheep, Haemonchus contortus is the most prevalent and virulent, resulting in health problems and production losses. Therefore, selecting sheep resistant to H. contortus is a suitable and sustainable strategy for controlling endoparasites in flocks. Here, 287 lambs of the native Brazilian Morada Nova hair sheep breed were subjected to two consecutive artificial infections with H. contortus and assessed for fecal egg count (FEC), packed cell volume (PCV), and live weight (LW). Forty-four animals ranked as having extreme resistance phenotypes were genotyped using the Illumina OvineSNP50v3 chip. A case–control genome-wide association study (GWAS) detected 37 significant (p < 0.001) markers in 12 ovine chromosomes in regions harboring quantitative trait loci (QTL) for FEC, Trichostrongylus spp. adults and larvae, weight, and fat; and candidate genes for immune responses, mucins, hematological parameters, homeostasis, and growth. Four single-nucleotide polymorphisms (SNP; OAR1_rs427671974, OAR2_rs419988472, OAR5_rs424070217, and OAR17_rs401006318) genotyped by qPCR followed by high-resolution melting (HRM) were associated with FEC and LW. Therefore, molecular markers detected by GWAS for H. contortus resistance in Morada Nova sheep may support animal selection programs aimed at controlling gastrointestinal nematode infections in flocks. Furthermore, genotyping of candidate genes using HRM qPCR may provide a rapid and efficient tool for animal identification.


Introduction
Morada Nova is a hair sheep breed characterized by high heat tolerance and adaptation to tropical climatic conditions [1]. Traditionally bred in the Brazilian Northeast region in small farms under extensive production systems, this native breed is also raised in other conditions and regions of Brazil [2]. Despite their small body size, these animals present interesting traits such as prolificacy, early sexual maturity, maternal ability, lack of reproductive seasonality, and resistance to gastrointestinal nematodes [3].
The barber pole worm, Haemonchus contortus, is the most prevalent and virulent nematode species that affects small ruminants in tropical regions [4]. H. contortus is a blood-sucking parasite of the abomasum, which causes anemia, weight loss, and death [5]. Considering the negative impact of nematodes on animal health and production, and the consequent economic losses [4,6], parasite resistance is an important trait that should be explored for animal selection in sheep. Furthermore, selection for resistance is a sustainable measure for parasite control [7], as it reduces production losses, reduces the use of anthelmintics, and decreases the infectiousness of pastures and the subsequent larval challenge to animals in the flock [8,9]. Several gene polymorphisms associated with resistance to gastrointestinal nematodes have been detected in Morada Nova sheep by candidate gene genotyping [10,11]; however, no genome-wide association studies (GWAS) have been identified for parasite resistance in this breed. Parasite resistance is a quantitative complex trait, and several genes may contribute to the final phenotype [12]. Therefore, GWAS can detect associated genes and polymorphisms, which may help to elucidate the genetic mechanisms affecting the inheritance of parasite resistance.
Therefore, a GWAS of H. contortus-resistant and -susceptible Morada Nova sheep was performed using a case-control design to identify genomic regions involved in resistance to gastrointestinal nematodes. Furthermore, to our knowledge, this was the first GWAS related to H. contortus resistance in a Brazilian native parasite-resistant Morada Nova sheep breed resulting in 37 significant markers located on 12 ovine chromosomes harboring regions previously associated with resistance-related traits and functional candidate genes. In addition to providing guidance for studies in other breeds, knowledge of the molecular mechanisms involved in parasite resistance will support the use and maintenance of Morada Nova sheep as a genetic stock for gene introgression and direct use in production systems, aiming to control nematode infections in flocks.
sustainable measure for parasite control [7], as it reduces production losses, red use of anthelmintics, and decreases the infectiousness of pastures and the su larval challenge to animals in the flock [8,9].
Several gene polymorphisms associated with resistance to gastroi nematodes have been detected in Morada Nova sheep by candidate gene gen [10,11]; however, no genome-wide association studies (GWAS) have been iden parasite resistance in this breed. Parasite resistance is a quantitative complex t several genes may contribute to the final phenotype [12]. Therefore, GWAS ca associated genes and polymorphisms, which may help to elucidate the mechanisms affecting the inheritance of parasite resistance.
Therefore, a GWAS of H. contortus-resistant and -susceptible Morada No was performed using a case-control design to identify genomic regions inv resistance to gastrointestinal nematodes. Furthermore, to our knowledge, this wa GWAS related to H. contortus resistance in a Brazilian native parasite-resistant Nova sheep breed resulting in 37 significant markers located on 12 ovine chrom harboring regions previously associated with resistance-related traits and fu candidate genes. In addition to providing guidance for studies in other breeds, kn of the molecular mechanisms involved in parasite resistance will support the maintenance of Morada Nova sheep as a genetic stock for gene introgression an use in production systems, aiming to control nematode infections in flocks.
Using Haploview, linkage disequilibrium (LD) blocks were detected for markers located in OAR 2 and 8, and strong LD was observed for markers in OAR 7, 11, and 18, (Figure 2).

Figure 2.
Haploview-generated linkage disequilibrium (LD) plots of markers in GWAS for H. contortus resistance in Morada Nova sheep. Three haplotype blocks (bold) were identified for singlenucleotide polymorphism (SNP) markers in OAR 2 (two LD blocks) and OAR 8. Red diamonds without a number indicate complete LD for SNP markers in OAR 7,11,and 18. In the 2 Mbp upstream and downstream regions of each significant marker, 683 unique genes (Supplementary Table S2), including 431 protein-coding genes, were identified. The same genomic regions harbored reported quantitative trait loci (QTL; Supplementary Table S3) for Nematodirus, H. contortus, and Trichostrongylus colubriformis fecal egg counts (FEC), larvae and adults of Trichostrongylus spp. in the abomasum and small intestine, weight, and fat traits (Table 1).
Functional annotation analyses revealed enrichment of sulfotransferase activity and RNA polymerase II transcription factor activity, sequence-specific DNA binding terms for gene ontology molecular function, visual perception, response to stimulus, and phototransduction terms for biological processes, and homeobox for domain (Supplementary Table S4). Glycosaminoglycan biosynthesis-heparan sulfate/heparin, beta-alanine metabolism, histidine metabolism, and glycolysis/gluconeogenesis were enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (Supplementary  Tables S4 and S5).
The biological roles of genes (with localization on chromosomes) potentially involved in nematode resistance in Morada Nova sheep (Table 1)  In the 2 Mbp upstream and downstream regions of each significant marker, 683 unique genes (Supplementary Table S2), including 431 protein-coding genes, were identified. The same genomic regions harbored reported quantitative trait loci (QTL; Supplementary Table S3) for Nematodirus, H. contortus, and Trichostrongylus colubriformis fecal egg counts (FEC), larvae and adults of Trichostrongylus spp. in the abomasum and small intestine, weight, and fat traits ( Table 1).
Functional annotation analyses revealed enrichment of sulfotransferase activity and RNA polymerase II transcription factor activity, sequence-specific DNA binding terms for gene ontology molecular function, visual perception, response to stimulus, and phototransduction terms for biological processes, and homeobox for domain (Supplementary Table S4). Glycosaminoglycan biosynthesis-heparan sulfate/heparin, beta-alanine metabolism, histidine metabolism, and glycolysis/gluconeogenesis were enriched KEGG (Kyoto Encyclopedia of Genes and Genomes) pathways (Supplementary Tables S4 and S5).
The biological roles of genes (with localization on chromosomes) potentially involved in nematode resistance in Morada Nova sheep (Table 1)  Real-time polymerase chain reaction (qPCR) followed by high-resolution melting (HRM) analyses confirmed the genotypes attributed using the chip in all 44 animals for OAR1_rs427671974, OAR2_rs419988472, OAR5_rs424070217, and OAR17_rs401006318, demonstrating 100% agreement. In addition, these four SNPs were associated with FEC ( Figure 3A), and OAR1_rs427671974 and OAR2_rs419988472 were associated with live weight (LW) ( Figure 3B). No association with packed cell volume (PCV) was detected for any of the four SNP markers. Real-time polymerase chain reaction (qPCR) followed by high-resolution melting (HRM) analyses confirmed the genotypes attributed using the chip in all 44 animals for OAR1_rs427671974, OAR2_rs419988472, OAR5_rs424070217, and OAR17_rs401006318, demonstrating 100% agreement. In addition, these four SNPs were associated with FEC ( Figure 3A), and OAR1_rs427671974 and OAR2_rs419988472 were associated with live weight (LW) ( Figure 3B). No association with packed cell volume (PCV) was detected for any of the four SNP markers.

Discussion
The present study employed a case-control GWAS for H. contortus resistance using OvineSNP50v3 chip genotyping of 44 ranked phenotyped (considering FEC, PCV, and LW) animals from a population of 287 Morada Nova sheep, with data from two experimental parasite challenges. GWAS detected 37 significant molecular markers, some in strong LD, on 12 ovine chromosomes in genomic regions harboring several functional candidate genes and QTLs for related traits. GWAS was validated by HRM qPCR genotyping of four significant SNP markers that were associated with FEC and LW. Mixed gastrointestinal parasite species are present under natural or field infection [12], and uninfected animals can be confounded by highly resistant animals [13]. Then, the use of two experimental parasite challenges, comprising monospecific artificial infection with H. contortus and repetitive phenotype collection, may have contributed to the suitability of the employed GWAS approach.
Significant markers were located in genomic regions harboring QTLs for parasite resistance, such as H. contortus, T. colubriformis, and Nematodirus FEC and Trichostrongylus spp. adults and larvae in the abomasum and small intestine, in OAR 2, 3, 7, 8, 11, 17, and 18. In addition, superposition of production trait QTLs, such as weight and fat, was detected Pathogens 2022, 11, 939 6 of 12 in OAR 1, 2, 3, 5, 6, 8, 11, and 18. An association with production traits was expected in the present study, as live weight was used to rank animals with extreme phenotypes. However, several studies based solely on FEC and nematode infection rates have also detected regions spanning production QTLs, suggesting a correlation between live weight and growth with resistance [14], as well as a natural selection effect to ensure developmental stability under the challenge from gastrointestinal nematode infection [13]. Some genomic regions detected in GWAS for H. contortus resistance in Morada Nova sheep did not harbor reported QTLs for related traits, as observed for OAR 5, 15, 17, and 20, confirming that different molecular mechanisms may regulate resistance in different breeds. Consequently, due to the genetically fragmented nature of sheep and goat populations, GWAS information derived from one breed cannot be extrapolated to others without proper validation [7].
Regarding the enriched terms detected in functional analyses, homeobox has previously been associated with adaptive immune response in cattle [15], and the glycosaminoglycan biosynthesis heparan sulfate/heparin pathway has been associated with viral invasion (reviewed by [16]). In addition, changes in the metabolism of beta-alanine and other amino acids, mainly due to the effect of microbiota on the host, were observed following H. contortus infection [17].
Several mechanisms can disrupt parasite establishment and confer host resistance to gastrointestinal nematodes [18,19]. Briefly, while feeding in the abomasum, H. contortus secretes and excretes antigens that stimulate host inflammatory, humoral, and cellular immune responses. Consequently, the recruited T-helper cells release cytokines, mainly interleukins, which activate IgE synthesis, eosinophils, mast cells, and globular leukocytes in the mucosa. These events are followed by B-cell activation and antibody production (IgA and IgG1). Additionally, mast cell and eosinophil inflammatory products, such as histamines, proteases, leukotrienes, and prostaglandins, lead to mucus production and smooth muscle contraction, which induce parasitic paralysis, elimination, or death. Furthermore, events favoring host homeostasis and coping with parasitic loads, such as protein and energy metabolism, hematological parameters, and body weight, have also been associated with resistance and/or resilience to gastrointestinal nematodes [8,20]. Based on these physiological functions, candidate genes were investigated in the genomic windows 2 Mbp upstream and downstream of each significant SNP marker detected via GWAS.
Among candidate genes related to the immune response, the expression of BTLA, an immunoglobulin superfamily member, and LEF1, an enhancer of T-cell receptor-alpha, was increased in peripheral blood mononuclear cells in Suffolk sheep following exposure to H. contortus larvae antigens [21]. In addition, BTLA expression in host T CD4(+) cells and innate leukocytes affected intestinal immunity and Strongyloides ratti infection [22]. Increased expression of CFI, which regulates the complement cascade, was detected in the abomasum of H. contortus-resistant sheep breeds [23][24][25], and CFI was also associated with H. contortus FEC [26]. Increased expression of CD96, an immunoglobulin superfamily member, has been detected in nematode-susceptible goats [27] and cattle [28]. CD200 [29], a glycoprotein containing two immunoglobulin domains, and NFE2L2 (or Nrf2) [30], which is involved in the response to oxidative stress, affected macrophage regulation and Leishmania infection. TREML1 (or TLT-1), TREM2, and TREM1 genes, which are involved in inflammatory, innate, and adaptive immune responses, were associated with tick resistance in cattle [31], and the role of TRML1 in clot formation and inflammatory or immune-induced bleeding has been reported [32,33]. Furthermore, increased expression of TREM1 was detected in sheep peripheral blood mononuclear cells during chronic infection with Fasciola hepatica [34]. FER, a tyrosine kinase involved in leukocyte recruitment, regulated the intestinal epithelial lipopolysaccharide barrier in response to bacteria [35]. DPP4, which is involved in metabolism and immune regulation, was associated with the innate immune response to virus [36] and demonstrated a role in hypoxia response in sheep [37].
Regarding the gastric mucosa and mucins, the GALNT10 gene, which drives mucintype O-glycan synthesis, is a paralog of GALNTL6, which was found by GWAS to be associated with gastrointestinal parasite resistance in sheep [12]. In addition, decreased expression of PGC, a component of the gastric mucosa, has been detected in the abomasum of a sheep breed resistant to H. contortus [23]. In cattle, increased expression of GCNT3, which plays a role in mucin-type glycoproteins, has been detected in the small intestine in response to Cooperia oncophora [38] and in the abomasum following Ostertagia ostertagi infection [39].
Considering that H. contortus is a hematophagous parasite, genes affecting hematological parameters in hosts may have a potential role in resistance. DLX1, a homeobox transcription factor, regulated the TGF-β superfamily during blood production [40]; ED-NRA, which encodes the receptor for endothelin-1, affected vasoconstriction in yaks [41]; CLEC14A controlled angiogenesis in mice [42]; and SH2B3, a negative regulator of cytokine signaling, was found to be associated with erythrocyte traits in sheep by GWAS [43].
Regarding homeostasis, PLA2G1B, a phospholipase A2 that regulates energy metabolism and inflammation in the intestine, was considered an endogenous anthelmintic that induces Heligmosomoides polygyrus and Nippostrongylous brasiliensis death in mice [44]. HS3ST5, a cell-surface heparan sulfate, acted as a receptor facilitating Trypanosoma cruzi [45] and Toxoplasma gondii [46] invasion. In addition, the expression of CCDC80, which enables glycosaminoglycan binding activity, was affected by Trypanosoma cruzi infection [47].
Some of the identified genes in this study were associated with growth, muscle development, and fat traits. Homeobox genes HOXD1, HOXD3, HOXD10, HOXD12, and HOXD13 have been detected by GWAS for muscularity in cattle [48], and PITX2 was related to growth in sheep [49] and weight in cattle [50]. DPPA2 was involved in myogenesis [51], KLHL41 was associated with skeletal muscle differentiation [52], and MYF5, a myogenic factor, was associated with growth in sheep [53]. Furthermore, MAPK7 was involved in adipocyte differentiation [54], and SREBF1 affected fat metabolism and deposition in sheep [55].
The analytical strategy employed by GWAS to detect H. contortus resistance in Morada Nova sheep was validated by association analyses of four significant SNP markers genotyped by HRM qPCR with FEC and LW in animals. No association with PCV was detected. However, while FEC and LW were used as factors multiplied by 0.4 to rank animals in extreme phenotypes, PCV presented a lower weight (0.2) in the equation. This decision was based on the fact that PCV presented lower variability among animals compared with FEC and LW, suggesting that Morada Nova sheep are resilient, rather than fully resistant to H. contortus [56]. The BB alleles of OAR1_rs427671974 and OAR2_rs419988472 were associated with lower FEC and higher LW, and the AB allele of OAR5_rs424070217 and the BB allele of OAR17_rs401006318 were associated with lower FEC. These markers, in addition to validating the GWAS results, may be used for the selection of resistant Morada Nova sheep using HRM qPCR genotyping.

Parasitological Tests and Phenotypic Classification of Animals
For phenotypic evaluation, 2 g of feces was collected from the rectum, mixed with 28 mL of sodium chloride-saturated solution, and evaluated in a McMaster chamber [57]. The total number of eggs was multiplied by 50 to obtain the fecal egg count (FEC). To determine packed cell volume (PCV), blood was collected into a heparin microcapillary tube and centrifuged at 1200 rpm for 5 min to determine the percentage of erythrocytes.
The DNA samples and phenotypic data (FEC, PCV, and live weight) used to rank animals with extreme phenotypes were obtained in 2017 and 2018 [58]. Briefly, 287 Morada Nova lambs (146 males and 141 females), the progeny of 7 rams from the Embrapa Pecuária Sudeste flock, were treated with monepantel (2.5 mg/kg; Zolvix ® ) to remove natural infection with gastrointestinal nematodes. After two FEC at 7-day intervals, animals were experimentally infected with 4000 H. contortus third-stage larvae (L 3 ) from an anthelmintic susceptible isolate [58]. and, after 15 days, they were subjected to a second parasitic challenge with 4000 L 3 H. contortus of the same isolate, followed by the sampling protocol described previously. At 42 DPI of the second parasitic challenge, mean FEC, PCV, and LW were calculated from the collected data. Subsequently, animals were ranked as extreme phenotypes (extremely resistant and extremely susceptible to H. contortus) based on the following equation: (1)

Genome-Wide Association Study (GWAS)
Based on phenotype ranking, 44 lambs were selected considering homogeneous phenotype distribution (21 resistant and 23 susceptible), sex (23 females and 21 males), and number of progenies (7)(8)(9)(10)(11)(12) from four rams with a large number of descendants. This resulted in 12 resistant females, 11 susceptible females, 9 resistant males, and 12 susceptible males. In addition, differences in rank mean value varied from 8 to 36 times between extreme phenotype progenies from each ram.
DNA was extracted by saline precipitation [59] from venous blood samples obtained from lambs. DNA integrity was confirmed via 1% agarose gel electrophoresis, and the DNA concentration and purity (260/280 absorbance ratio between 1.8 and 2.0) were estimated using a NanoDrop 2000 spectrophotometer. The samples were genotyped using the Illumina OvineSNP50v3 chip (Supplementary Table S6) at the Centro de Genômica Funcional at ESALQ/USP in Piracicaba, SP, Brazil.

Bioinformatics and Functional Annotation
Genotyping chip data were subjected to quality control filtering, with 47,782 singlenucleotide polymorphisms (SNP) retained for analyses by case-control GWAS [60]. PLINK was used to filter markers, and SNPs were removed if they had a minor allele frequency (MAF) < 0.01, a call rate < 90%, a GC score < 0.6, and a deviation from Hardy-Weinberg equilibrium (HWE) < 10 −15 . An efficient mixed-model association expedited (EMMAX) was used to identify nominal p-values for phenotype-genotype interactions, and then PLINK (using -assoc, -perm, and -adjust functions) was used for Bonferroni testing. Linkage disequilibrium (LD) between SNP markers was calculated using Haploview [61].
ENSEMBL_Gene_ID was used as a gene list in the functional annotation tool of the DAVID Bioinformatics Resources 2021 update (https://david.ncifcrf.gov/home.jsp accessed on 29 April 2022) with Ovis aries genome as the background. DAVID was used to identify the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways [62]. Considering that an enrichment score (ES) of 1.3 is equivalent to p = 0.05 in Fisher's exact test [14], ES > 1.3 was used as the threshold to detect significantly enriched terms. Classification in GeneCards version 5.9 (https://www.genecards.org/ accessed on 9 May 2022) and a literature review were used to identify functional candidate genes.
The effects of SNPs on FEC, PCV, and LW were analyzed as described by [10]. Briefly, mean FEC values were normalized by orderNorm (ORD) using the bestNormalize package in R. Sex, age at weaning, group, ram, birth type, and age of dam were used as fixed effects, and each SNP was individually tested by ANOVA, using the aov() function, followed by Tukey's test adjusted for unbalanced group sizes (HSD.test()).

Conclusions
GWAS to determine H. contortus resistance profiles using 44 extreme phenotype Morada Nova sheep genotyped by Illumina OvineSNP50 chip resulted in 37 significant markers in OAR 1,2,3,5,6,7,8,11,15,17,18,and 20. Marker genomic regions harbored QTLs for FEC, adults and larvae of Trichostrongylus spp. in the abomasum and small intestine, weight, and fat traits, in addition to functional candidate genes related to the immune response, gastric mucin, hematological parameters, homeostasis, growth, and fat deposition. HRM qPCR genotyping of four molecular markers (OAR1_rs427671974, OAR2_rs419988472, OAR5_rs424070217, and OAR17_rs401006318) validated their association with FEC and live weight. Thus, the obtained results can be used for the selection of native Morada Nova sheep with the aim of controlling gastrointestinal nematode infection in flocks and increasing parasite resistance in production systems.
Supplementary Materials: The following supporting information can be downloaded at: https:// www.mdpi.com/article/10.3390/pathogens11080939/s1, Table S1: EMMAX p-values for phenotypegenotype interactions, Table S2: Genes located in the 2 Mbp upstream and downstream regions of each significant SNP marker on ovine chromosomes (OAR), Table S3: Quantitative trait loci (QTL) in the 2 Mbp upstream and downstream regions of each significant SNP marker on ovine chromosomes (OAR), Table S4: Enriched terms in functional annotation analyses using DAVID, Table S5: Enriched KEGG pathways, Table S6: OvineSNP50v3 Illumina, and Table S7: Fragment sizes and primers used in HRM qPCR for genotyping of SNP markers in four ovine chromosomes (OAR). Institutional Review Board Statement: The animal study protocol was approved by the Ethics and Animal Experimentation Committee of the EMBRAPA PECUÁRIA SUDESTE (protocol code 04/2017 approved on 08/01/2017 and 01/2020 approved on 09/15/2020).

Informed Consent Statement: Not applicable.
Data Availability Statement: Additional data, not presented in Supplementary Materials, are available on request from the corresponding author.

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