Genome-wide insights on gastrointestinal nematode resistance in autochthonous Tunisian sheep

Gastrointestinal nematode (GIN) infections have negative impacts on animal health, welfare and production. Information from molecular studies can highlight the underlying genetic mechanisms that enhance host resistance to GIN. However, such information often lacks for traditionally managed indigenous livestock. Here, we analysed 600 K single nucleotide polymorphism genotypes of GIN infected and non-infected traditionally managed autochthonous Tunisian sheep grazing communal natural pastures. Population structure analysis did not find genetic differentiation that is consistent with infection status. However, by contrasting the infected versus non-infected cohorts using ROH, LR-GWAS, FST and XP-EHH, we identified 35 candidate regions that overlapped between at least two methods. Nineteen regions harboured QTLs for parasite resistance, immune capacity and disease susceptibility and, ten regions harboured QTLs for production (growth) and meat and carcass (fatness and anatomy) traits. The analysis also revealed candidate regions spanning genes enhancing innate immune defence (SLC22A4, SLC22A5, IL-4, IL-13), intestinal wound healing/repair (IL-4, VIL1, CXCR1, CXCR2) and GIN expulsion (IL-4, IL-13). Our results suggest that traditionally managed indigenous sheep have evolved multiple strategies that evoke and enhance GIN resistance and developmental stability. They confirm the importance of obtaining information from indigenous sheep to investigate genomic regions of functional significance in understanding the architecture of GIN resistance.

Estimation of genetic diversity. Expected (H E ) and observed (H O ) heterozygosity, effective population size (N E ), and patterns of LD decay were investigated for the two cohorts (infected and non-infected). H E and H O were calculated using PLINK v1.9. Pair-wise LD was evaluated using the correlation coefficient (r 2 ) between alleles at two separate SNP loci using PLINK v1.9 with default settings. Following Sved 25 , N E over generation time was estimated with the equation N Et = (1/4c) (1/r 2 − 1), where N Et is the effective population size t generations ago (t = 1/2c); r 2 is the LD between pairwise SNPs; and c is the genetic distance in Morgan between pairs of SNPs.
For each cohort, two measures of inbreeding were calculated; (1) SNP based inbreeding coefficient (F) calculated with PLINK v1.9 and (2) runs of homozygosity (ROH) based inbreeding coefficient (F ROH ). For the latter, ROH streatches were first computed using the "detectRUNS" package in R 26 . The F ROH , was computed as the ratio of the total length of ROH to the length of autosomes (2.45 Gb) 27 .
Three estimates of ROH were calculated taking into account three genomic distance categories, ROH < 5 Mb, 5 Mb < ROH < 10 Mb, and ROH > 10 Mb, indicative of ancestral (more than ten generations), middle (5-10 generations), and recent (within five generations) inbreeding, respectively 28 . For the calculation, ROHs were identified for each individual using PLINK v 1.9 with the following parameters: (1) the minimum number of SNPs in a sliding window was 50; (2) the minimum ROH length was set to 1 Mb to eliminate the impact of strong www.nature.com/scientificreports/ LD; (3) each ROH spanned a minimum of 80 consecutive SNPs; (4) one heterozygous and five missing calls per window were allowed to avoid false negatives that may arise due to occasional genotyping errors or missing genotypes; (5) the minimum SNP density was set at one SNP every 100 kb, and the maximum gap between consecutive SNPs was set to 1 Mb.
Population structure analysis. The 335,070 SNPs were used to infer genetic structure and divergence using three methods. Principal component analysis (PCA) was performed using PLINK v1.9 and the first two PC's were plotted to visualize individual relationships. The proportion of shared ancestry between individuals was inferred with the unsupervised mode of ADMIXTURE tool v1.30 29 . This mode does not assume any background information on the number and frequency of alleles in ancestral populations. The ADMIXTURE tool was run with values of K increasing gradually from 1 to 6, to derive cross-validation (CV) errors. The lowest value of the CV error indicates the most likely number of ancestral populations. Five runs were performed for each K. To test for consistency in the ADMIXTURE results, we repeated the ADMIXTURE runs using three datasets of 6,000 SNPs, each drawn at random without replacement from the 335,070 SNPs. The pairwise allele sharing distance (ASD) 30 was computed using the program "asd". A distance matrix of all pairwise ASD dissimilarities, calculated as 1-ASD, was generated and subjected to hierarchical clustering with the Neighbour-Joining algorithm to yield an individual clustering dendrogram. ASD calculation does not require estimates of allele/ genotype frequencies making it valid when sample sizes are small. It is also suitable for detecting outliers and is robust to high LD among SNPs.

Detection of selection signatures and association analysis.
To investigate the molecular genetic basis underlying natural variation in GIN infection, we investigated genome-wide distributions of ROH in each individual in the infected and non-infected cohorts using PLINK v1. 9. For this analysis, we used the 540,528 SNPs that passed the quality thresholds. Using the "detectRUNS", we counted the number of times a given SNP occurred in the identified ROH in each cohort and presented a Manhattan plot of all the tested SNPs against their autosomal positions. The most frequently observed SNPs in ROHs occurring at the top 25% of the empirical distribution were taken as the most significant loci underlying an ROH under selection. To identify regions of ROH that are likely associated with variations in GIN infection, we compared the ROH regions between the infected and non-infected cohorts and identified the ones that were specific to the non-infected cohort. We performed the logistic regression (LR) GWAS, F ST and XP-EHH analyses to investigate further the genome regions associated with variation in GIN infection. LR-GWAS was performed with PLINK v1.9 using the noninfected cohort as the test sample and the infected cohort as the control. To obtain the 99% confidence intervals for the estimated parameters, the "-ci 0.99" and "-covar" options were invoked and Fisher's exact test was used to generate the p-values considering location, breed, sex, and age as covariates. The raw p-values were subjected to Bonferroni correction to control the likelihood of any false positive results among the markers identified to be under selection. The corrected p-values were standardized and the − Log 10 (p-value) of 4.25 (the top 0.001) was set as the cut-off threshold to identify candidate regions. The estimations were summarized in 200 Kb window sizes and the genes falling within the candidate regions were identified. The package ' qqman' in R version 3.5.1 was used to generate the Manhattan plot.
The population differentiation statistic, F ST , was used to investigate regions of the genome that have diverged between the two cohorts. The unbiased pairwise F ST 31 was computed using the HIERFSTAT Package of R 32 using a window size of 200 Kb and a window-step size of 100 Kb. Windows with less than five SNPs were excluded from the analysis. The F ST values were then standardized by Z transformation following Ahbara et al. 33 . Windows falling within the top 0.001% of the F ST values in each chromosome were considered the putative candidates under divergent selection.
The cross-population extended haplotype homozygosity (XP-EHH) 34 was used to compare expected haplotype homozygosity (EHH) and integrated haplotype score (iHS) between the two cohorts to detect selection and its direction. XP-EHH scans SNPs that are homozygous in one population but polymorphic in the other through pairwise comparison of EHH scores. Positive XP-EHH scores indicate selection in the test sample, while negative scores indicate selection in the control. The XP-EHH scores were estimated as: where I A is the integrated EHH value of the test population and I B is the integrated EHH value of the reference population. Haplotype phasing was inferred for each cohort simultaneously on all SNPs using BEAGLE v3.3.1 35 . The XP-EHH test was performed with the "rehh" package of R 36 and the raw XP-EHH scores were standardized to a distribution with zero mean and unit variance. Selection candidates were considered as the regions contained in any of the 200 Kb windows with a significance threshold of p < 0.001; this equates to an XP-EHH value of 4 at the default settings of "rehh" estimation.
Functional annotation of candidate regions. The candidate regions identified by ROH were analysed and the ones that were specific to the non-infected cohort identified. We also analysed the ROH regions of the non-infected cohort, LR-GWAS, F ST and XP-EHH candidate regions and the ones that overlapped between at least two approaches were identified and merged using Bedtools v.2.28.0 37 . Genes that were spanned by the ROH (non-infected cohort) and overlapping candidate regions were retrieved using the Biomart/Ensembl (http:// www. ensem bl. org/ bioma rt) tool based on the Ovine v3.1 reference genome assembly. The set of genes identified in the candidate regions were assessed for biological enrichment gene ontology (GO) and KEGG Pathway (www.  Table 2. The most and least frequent ROH length was observed for the 0-5 Mb and > 10 Mb length categories, respectively. The average length of ROH per animal was highest, in the 5-10 Mb length category for non-infected cohort and > 10 Mb length category for the infected cohort. The two cohorts show similar patterns of LD decay over genomic distance although the non-infected cohort showed overall lower LD (r 2 ) (Fig. 1a). In general, the pattern of LD decay shows higher LD at shorter distances which declines rapidly and plateaus off from 0.4 Mb for both cohorts (Fig. 1b). The trend in N E over generation time was the same for both cohorts (Fig. 1c). They are characterised by an increase in N E from 1000 generations ago, which attains maximum value at approximately 330 generations ago, and then declines to the present time. Generally, the non-infected cohort had higher N E across all generations.
Population structure analysis. Population structure and relationship was investigated using PCA ( Fig. 2a), ADMIXTURE tool (Fig. 2b,c) and ASD phylogeny (Fig. 2d). The first and second PC's of the PCA explained, respectively 1.74% and 1.45% of the total genetic variation. The study animals did not differentiate into distinct genetic groups/clusters that correspond to their infection status. Following ADMIXTURE analysis, the lowest CV error was at K = 1 (Fig. 2b) which suggests no genetic differentiation. The ADMIXTURE plot reveals that a similar pattern of genomic composition characterizes the two cohorts at 2 ≤ K ≤ 6 ( Fig. 2c). This pattern was replicated when ADMIXTURE was performed using three datasets comprising 6000 SNPs each, drawn at random without replacement from the 335,070 SNPs. The ASD phylogeny also showed no genetic stratification (Fig. 2d).
Genome-wide selection signature analysis. The ROH analysis identified 60 ROH regions in the two cohorts ( Fig. 3a,b) that spanned 311 genes. The LR-GWAS, F ST and XP-EHH identified 346, 32, and 68 regions (Fig. 4a-c), respectively which spanned 673, 152, and 295 genes. These 446 candidate regions overlapped with 645 genes (Supplementary Table S1) of which 71 were found in candidate regions that were identified by at least two methods (Fig. 5).
We considered the ROH analysis as a method that identifies selection signatures within a cohort. We therefore compared the ROH results of infected and non-infected cohorts (Fig. 3a,b). This identified 23 ROH regions www.nature.com/scientificreports/ that were specific to the non-infected cohort and which spanned 80 putative genes (Table 3) of which 30 remain uncharacterised/unannotated (prefixed with "ENSOARG"). Six of the 23 candidate regions occurred on OAR 1, 2, 3, 6, 10 and 14 and overlapped with FECGEN, LATRICH_2 and NFEC QTLs ( Table 3) that are associated with health traits and in particular, parasite resistance. The LR-GWAS was used to identify candidate regions and possible SNPs associated with GIN infection, while F ST , and XP-EHH were also used to investigate selection signatures by contrasting the non-infected and infected cohorts. These three approaches identified 30 candidate regions that overlapped between at least two methods across 17 autosomes (Table 4). When the ROH regions that are specific to the non-infected cohort are considered, the four methods identify 35 candidate regions overlapping between at least two approaches across 19 autosomes and span 121 genes including 11 which remain uncharacterized (Table 4). Of the 35 candidate regions, five that were identified by ROH to be specific to the non-infected cohort overlap with at least one region that was identified by either LR-GWAS, F ST and/or XP-EHH and span 13 genes (Table 4) while one region found on OAR 5 (region number 14) was identified by all the four methods (Table 4). This region spans four genes (LOC101104745, PCDHGA1, PCDHGA2, ENSOARG00000000218) and one QTL trait (BIRTH_WT, Body weight (birth)). None of the four genes are associated, directly or indirectly, with endoparasite resistance. Nineteen regions found across OAR 1, 3, 6, 8, 11, 12, 13, 14, 17, 18 and 21 span FECGEN, TFEC_1, HFEC, NFEC, LATRICH_2, IGA, SAOS, WORMCT, PEPSL and CEOSIN QTLs that are linked to response to GIN infection (Table 4). Ten regions overlapped with production (growth) QTL traits (BW, BIRTH_WT, BONE_WT) across OAR 2, 4, 5 and 6, four regions overlapped with meat and carcass traits associated with fatness QTLs (HCWT, FAT_WT, LMYP) across OAR 2 and OAR 10, and anatomy QTLs (BONE_WT) on OAR 24 (Table 4).
Functional enrichment analysis was first tested in the pool of 51 genes, excluding uncharacterized genes, that are present in the 23 ROH candidate regions that were specific to the non-infected cohort (Table 3). A second functional enrichment analysis was performed with the set of 110 genes, excluding 11 uncharacterized ones that were present in the 35 candidate regions that overlapped between at least one ROH, LR-GWAS, F ST and XP-EHH regions (Table 4). We found two functional term clusters that were significantly (enrichment score > 1.5) enriched for the genes present in the ROH regions of the non-infected cohort (Table 5). These were associated with "immune system process" and "cytokine receptor interaction". The second analysis resulted in three significantly (enrichment score > 1.5) enriched clusters. The first cluster was associated with the GO biological terms "Carboxylesterase type B" and "Hydrolase activity". The "Apical plasma membrane" term, which has a role www.nature.com/scientificreports/ in intestinal innate immunity, and "Integrin signalling" that functions in immune cells were among the most significant GO biological terms in the second and third clusters, respectively.

Discussion
GIN infections and associated gastroenteritis impact negatively the production efficiency of ruminant livestock, and their management is essential to meet future demand for animal source foods. The observation of large variations in GIN infection suggests variability at the genome level that underpins inter-animal variability in      23 study. The aim was to investigate signatures of variability in GIN infection and resistance in autochthonous Tunisian sheep that are managed under natural grazing and no history of anthelmintic intervention. We applied four methods to detect genomic regions associated with GIN resistance. Although our study yielded some interesting findings, we acknowledge that studies that use naturally acquired field exposure to deliver a challenge always run the risk of uninfected animals being a combination of truly highly resistant animals, animals that never saw infection and animals that have had infection cleared due to chemical treatment. Additionally, producer provided information is at times not always accurate. These factors may serve to dilute the certainty that the two groups are in fact functionally dissimilar. Although we expected the two cohorts to show differences in genetic diversity and structure, this was not the case; they showed similar levels of genetic diversity and inbreeding and the three clustering algorithms provided corroborating evidence of lack of stratification that was consistent with infection status. The lack of genetic differentiation and structure was also reported between prolific and non-prolific cohorts of Bonga sheep from Ethiopia 41 and in the Brazilian Santa Inês breed which shows variability in resistance to GIN infection 19 . The absence of genetic stratification appears to be a characteristic of sheep from the Maghreb as it has also been observed in sheep populations from Algeria 42 and Morocco 43 implicating extensive intermixing and cross-mating. There are four possible explanations for the lack of genomic stratification corresponding to infection status.
(1) The selection pressure driving the differences in GIN infection is/has not been stringent or long enough to Table 3. Candidate regions that are specific to the non-infected cohort as identified by ROH analysis. www.nature.com/scientificreports/ result in differentiation at the genome level; (2) the level of parasite infection may be too low to result in meaningful genomic variation; (3) it may point to a lack of farmer-driven selection that is biased towards the use of uninfected animals for breeding; and/or (4) there could be a high natural selection pressure whereby all animals effectively have the resistance genotype and therefore rather than a distribution, spanning susceptible to resistant, there is a skewed distribution that spans "moderately resistant" to "highly resistant" individuals. One or all of these reasons may explain the lack of genetic differentiation between the two cohorts. The level of genetic diversity in the two cohorts is higher than that observed in commercial breeds, but it falls within the range of values reported in other sheep found in Africa and China [44][45][46] .
We assessed the average length and number of ROH per individual at three genomic distance categories and the trends in N E over generation time for each cohort. The two cohorts had a high frequency of ROH in the shorter (0-5 Mb) length category reflecting older evolutionary events such as past or ancestral inbreeding. The two cohorts showed similar patterns and trends in LD decay and N E . However, the non-infected cohort had higher values of N E across all generations. This difference is difficult to explain but we presume it may be the  www.nature.com/scientificreports/ result of higher mortalities in the infected cohort. Both cohorts show a gradual increase in N E from 1000 generations ago. This is followed by a drastic decline from around 330 generations ago to the present time suggesting a bottleneck event. A similar trend was observed in Ethiopian and Sudanese local sheep 47 . Assuming a generation time of three years for traditionally managed local sheep, it can be inferred that the increase in N E started around 3000 years ago. The start of the decline in N E 330 generations back translates to around 990 years ago. Between 240 and 1140 years ago, three favourable climatic periods interspersed with short extreme dry spells were experienced in the continent 48 . We suggest that the sheep populations most likely thrived when conditions were optimal but shrunk during subsequent droughts. The footprints of these demographic perturbations appear to have been retained in the genomes of the indigenous sheep. Our analysis revealed 35 candidate regions, spanning 110 genes, that overlapped between at least two out of the four methods we used to detect genomic regions associated with differences in endo-parasite infections. The identification of overlapping candidate regions that are under selection by different approaches suggests plausible evidence for selective influences at the genome 49 . The convergence of our results points to the reliability of our findings and suggest that they are unlikely to be the outcome of chance events or analytical artefacts.
Both the within (ROH) and between-population (LR-GWAS, F ST , XP-EHH) approaches identified genomic regions that could be driving GIN resistance in autochthonous Tunisian sheep. At least two or more methods identified the same candidate regions that colocalized with the FECGEN, TFEC_1, HFEC, NFEC, LATRICH_2, IGA, SAOS, WORMCT, PEPSL and CEOSIN QTLs that have been associated with health traits and in particular parasite resistance, immune capacity and disease susceptibility in sheep showing resistance to GIN [50][51][52] . This result indicates that some common QTLs underly parasite resistance traits in sheep and support a role for convergent evolution in driving host GIN resistance. For instance, the FECGEN QTL which results in reduced FEC in unmanaged, naturally-parasitized domestic sheep 53 and in the primitive Soay sheep 20 , has been associated with a microsatellite allele (o(IFN)-γ 126 ) found in the first intron of the interferon gamma (IFN)-γ gene and with increased titre of Teladorsagia circumcincta-specific IgA. It has been shown that the effects of the o(IFN)-γ 126 allele and IgA on FEC are not independent, and that IgA may mediate the (IFN)-γ effect on FEC 20 . Furthermore, Table 5. Enriched functional term clusters following DAVID analysis of genes identified by ROH analysis in candidate regions that were specific to the non-infected cohort and those that overlapped between at least two methods of detecting selection signatures. www.nature.com/scientificreports/ the IGA QTL could be related to early response to incoming larvae, whereas the FECGEN QTL may be associated with the ability to avoid the development of adult parasites 52,54 23 . QTL and GWAS studies suggest that GIN infections can evoke several host responses that enhance innate and acquired immune responses, gastric mucosal protection, haemostasis pathways, delay in parasite development and reduction in number of eggs produced by GIN 16 . These manifestations depend on the nematode species, parasite exposure, and host-specific factors including age, sex, genetic make-up, hormonal and nutritional status 55 . The animals analysed in our study are grazed under natural open pastures and are exposed to a wide range of GIN. Our expectation, therefore, was that they would exhibit a large repertoire of host defence mechanisms which may be reflected in the functions of the putative genes present in the candidate regions. We observed two organic cation transporters, SLC22A4 (OCTN1) and SLC22A5 (OCTN2) which are carriers of hydroxyurea, in an overlapping candidate region on OAR 5 (Table 4). Hydroxyurea can inhibit ribonucleotide reductase (RR) and is retained at high concentrations in tissues of various mammalian species including sheep 56 . Compared to viral or bacterial ribonucleotide reductase, the mammalian RR is less susceptible to inhibition by hydroxyurea 56 . This action can be critical in inhibiting viral and bacterial replication without affecting mammalian cellular growth. Hydroxyurea can therefore offer a level of immunity against GIN infection as part of the innate immune defence system which may inhibit endoparasites in the gastrointestinal tract of non-infected animals. Some variants present in SLC22A4 and SLC22A5 have also been associated with Inflammatory Bowel and Crohn's disease's in humans 57,58 .
We observed BMPR2 (Bone Morphogenetic Protein Receptor Type 2) in an overlapping candidate region on OAR 2 ( Table 4). Ligands of this gene are members of the TGFβ superfamily whose homologues are key players in inducing immunological tolerance 59 . Elevated expression of TGFβ have been observed in mammalian hosts that mount an effective immune response against GIN [59][60][61][62] . BMPR2 was found to be highly expressed in mesenteric lymph nodes of cattle that were selected for resistance to intestinal nematodes 63 suggesting association with parasite resistance. BMPR2 also occurred in a cluster of genes that were upregulated following Eimeria acervulina infection in chicken and associated with the functional term "Inflammatory Response" in the "Disease and Disorders" category 64 . Several mechanisms may operate to increase the levels of TGFβ under parasite infection including, host homeostasis to minimize immunopathology under chronic infection; pathogens triggering TGFβ production or activation; or parasite mimicry of the host cytokine to drive the same pathway as the hosts TGFβ 59 . Examples of all three have been reported in diverse sheep breeds 59 and the one that operates in autochthonous Tunisian sheep remains unknown calling for further investigation.
The IL-4 (Interleukin-4) and IL-13 (Interleukin-13) occurred in a candidate region on OAR 5 that was specific to the non-infected cohort. IL-4 plays a crucial role in the differentiation of naive T helper (Th) cells into Th2 effector cells which promote humoral immunity and provide protection against intestinal helminths 65 . In sheep, Jacobs et al. 66 observed that impaired IL-4 signalling promoted the establishment of H. contortus and increased larval burden. An increase in IL-13 in intestinal lymph cells was observed in sheep selected for increased resistance to nematodes 67 . Two genes RUFY4 (RUN and FYVE Domain Containing Protein 4) and VIL1 (Villin 1) and two IL-8 receptors (CXCR1 and CXCR2) occurred in a strong candidate region on OAR 2 that was also identified by ROH to be specific to the non-infected cohort. RUFY4 is a positive regulator of autophagy and is expressed in a cell-specific manner or under specific immunological conditions associated with IL-4 expression 68 . VIL1, has been shown to protect gut epithelial cells by decreasing epithelial damage 69 . It was among the top 30 proteins that were differentially regulated between resistant and susceptible sheep and may play an important role in maintaining the epithelial integrity of abomasal mucosa in response to haemonchosis 70 . It can therefore be hypothesized that the interaction between RUFY4 and VIL1 with IL-4 may enhance long-term endoparasite resistance in sheep. The two IL-8 receptors were found to be expressed and distributed in normal and morphologically damaged large intestines, suggesting that IL-8 may play an important role in mediating immune response in gastrointestinal tract, beyond that of potentiating neutrophil recruitment and inflammation 71 following intestinal epithelial damage by endoparasites.
Protective immunity against GIN and subsequent parasite expulsion is known to be mediated by Th2 immune response that is orchestrated by the secretion of cytokines such as IL-4, IL-5, IL-10 and IL-13, which promote the recruitment of eosinophils, basophils, mast cells, goblet cell hyperplasia and concurrent mucus and antibody production [72][73][74] . Nematode expulsion is known to rely on smooth muscle contraction and increased mucus production 75 . The latter was found to be mediated by intelectin 1 (ITLN), a calcium dependent lectin, that is expressed by Paneth cells in the small intestine of mouse 76 and mucus neck cells in the abomasum of sheep 77 . The expression of ITLN in sheep goblet cells was found to be upregulated by IL-4 77 . High expression of ITLN was also found to be induced during nematode expulsion in Trichinella spiralis and Nocardia brasiliensis in rodents 78,79 and T. circumcincta infections in sheep 77,80 . Th2 response has also been described as a mediator for acute wound healing during helminth infection 81 . We therefore suggest that the upregulation of ITLN by IL-4 may play a role in enhancing GIN resistance in autochthonous Tunisian sheep through parasite expulsion and enhancing wound healing.
Other possible candidate genes that have been associated with GIN resistance and occur in our overlapping candidate regions include FGF2, FAM78B, SPATA5, SPPL2B and FAM96B. FGF2 also known as basic fibroblast growth factor (bFGF and/or FGF-β), is a growth factor and signaling protein that can synergize with IL-17 in the gut to activate the ERK pathway and induce genes for repairing damaged intestinal epithelium 82 . The gut epithelium is essentially the first line of defence against microbiota and pathogens and therefore, it plays a critical www.nature.com/scientificreports/ role in enhancing mucosal immunity. FAM78B (Family with Sequence Similarity 78-Member B) was in an overlapping region on OAR 1. In a genome-wide association study of endo-parasite phenotypes, FAM78B was found in a region on bovine chromosome 3 that spanned SNPs that were most strongly associated with antibody response to O. ostertagi 40 . SPATA5 (Spermatogenesis Associated 5) was among five genes that were included in the most significant functional term cluster linked with immunity-related and cell-proliferation processes in an investigation of genomic regions and genes for gastrointestinal parasite resistance in Djallonke sheep 83 . SPPL2B (Signal peptide peptidase-like 2B), a member of the signal peptide peptidase-like protease (SPPL) family, localizes to endosomes, lysosomes and plasma membrane and plays a role in enhancing innate and adaptive immunity by cleaving TNFα in activated dendritic cells that triggers IL-12 production 84 . In humans and chicken, there is strong evidence linking polymorphisms in FAM96B (family with sequence similarity 96 member B) with gastrointestinal and metabolic diseases, and with the development of the digestive system and disorder networks 85,86 . These findings suggest that several mechanisms may be employed to trigger and sustain GIN resistance in sheep under traditional grazing and exposure to multiple GIN species. Our analysis revealed several candidate genomic regions spanning a number of production (growth), and meat and carcass (fatness and anatomy) QTLs. This result was unexpected given our analytical strategy. Kipper et al. 87 reported a reduction of 5% and 31% in average daily feed intake and average daily weight gain, respectively in parasitized pigs. Endoparasites tend to limit host nutrient availability by reducing host food intake, digestion, absorption and nutrient assimilation, resulting in nutritional deprivation and destabilization of host growth and development 88,89 . To counter against these negative effects, we suggest that natural selection may be acting on the regions spanning growth QTLs to ensure growth and developmental stability of the study populations under GIN challenge. The parallel selection signatures overlapping the growth QTLs may therefore be adaptive strategy that ensures the survival of the study populations.
In conclusion, we assessed the diversity and population structure of two cohorts (infected and non-infected) of autochthonous sheep from Tunisia that were classified based on FEC levels. The two cohorts were characterised by similar levels of genetic diversity and the same genome background suggesting common history and genetic admixture. Four methods of detecting selection signatures identified regions of the genome that were most likely associated with GIN resistance suggesting that the animals have established a certain level of immunity under natural challenge. The functions of the putative genes and the overlapped QTLs suggest multiple strategies, including host immune response, intestinal epithelium damage repair, mucus production and parasite expulsion, play a role in GIN resistance in sheep. This may be due in part to the fact that GIN resistance is the net outcome of many physiological and immunological pathways, and thus resistance in animals could be owing to variation at a large number of loci. Furthermore, natural selection is also acting concurrently on regions spanning growth related QTLs to ensure developmental stability under GIN challenge. Although the data used in our study is relatively small, we believe this does not compromise the integrity of our findings as it is compensated for by the high density of the marker loci used. We have also limited our discussion to candidate regions that were specific to the non-infected cohort and those that overlap between at least two analytical approaches. This research confirms the importance of obtaining data from local sheep populations managed in their production environments to gather information on genomic regions of functional significance in GIN resistance.