Genome Variation in Tick Infestation and Cryptic Divergence in Tunisian Indigenous Sheep

doi:10.21203/rs.3.rs-789829/v1

Download PDF

Research Article

Genome Variation in Tick Infestation and Cryptic Divergence in Tunisian Indigenous Sheep

https://doi.org/10.21203/rs.3.rs-789829/v1

This work is licensed under a CC BY 4.0 License

Journal Publication

published 28 Feb, 2022

Read the published version in BMC Genomics →

You are reading this latest preprint version

Background: Ticks are obligate haematophagous ectoparasites considered second to mosquitos as vectors and reservoirs of multiple pathogens of global concern. Individual variation in tick infestation has been reported in indigenous sheep, but the genes regulating the trait are poorly understood.

Results: Here, we report 397 genome-wide signatures of selection overlapping 991 genes from the analysis, using four methods (ROH, LR-GWAS, XP-EHH, F_ST), of 600K SNP genotype data from 170 Tunisian sheep exhibiting high and low resistance to ticks. We considered 45 signatures detected by consensus results of at least two methods as high-confidence selection sweep regions. These spanned 104 genes which included immune system function genes, solute carriers and chemokine receptor. One region spanned STX5, that has been associated with tick resistance in cattle, implicating it as a prime candidate in sheep. We also observed RAB6B and TF in a high confidence candidate region that has been associated with growth traits suggesting natural selection is enhancing growth and developmental stability under tick challenge. The analysis also reveals fine-scale genome structure suggesting the existence of cryptic divergence in the Tunisian sheep.

Conclusion: Our findings provide genomic reference that could enhance our understanding of the genetic architecture of tick resistance and cryptic divergence in indigenous sheep.

Epigenetics & Genomics

Genome variation

tick infestation

cryptic divergence

Tunisian indigenous sheep

Ticks are among the most ancient acari and possibly the earliest to evolve blood-feeding capacity through the development of immune-active proteins and lipids that induce vasodilatory, anti-haemostatic and immunomodulatory activities that guarantee successful engorgement and acquisition of a blood meal [1–5]. Compared to other arthropod vectors, ticks transmit a greater variety of pathogenic microorganisms including protozoa, bacteria and viruses, implicated in severe infections in humans, wild and domestic animals [6, 7]. Globally, ticks and tick-borne infections (T-TBI) are a major health impediment to efficient livestock production and reproductive performance. These translates to reduced farm incomes, increased consumer and producer costs, and disruption of local and global markets, and trade through the imposition of quarantines and market bans.

Prophylactic use of acaricides is the commonly employed strategy to control and eradicate ticks [8, 9]. The (over)use of acaricides has however imposed selection pressure that has resulted in the development of acaricide-resistant tick strains, persistent environmental contamination, infection of farmers and treated animals, retention of chemical residues in livestock products, and increased costs of existing, and of developing and manufacturing of new more potent, acaricides [9–13]. Anti-tick vaccines are a promising option, but they do not confer protection against multiple tick species [7, 14, 15] which is a common occurrence in most regions, especially the (sub)tropics [16–18]. These factors are driving the search for alternative sustainable tick control strategies such as the use of tick-resistant animal genotypes. The host’s natural resistance to ticks can be exploited as a long-term sustainable control strategy. It offers a low-cost, permanent solution that demands no extra resources and costs to generate a unit of product [19]. Moreover, the inherent resistance targets all tick species.

Genetic variation in resistance to parasites has been demonstrated in many livestock species with investigations on individual variation in resistance to T-TBI being the subject of intense research in cattle (see review by [20]) compared to other livestock species. Most of the studies have shown that an effective immunological response (resistance/tolerance) to T-TBI, is genetically determined with a heritability of between 0.05 and 0.42 (see review by [20]). For haemopathogens, the immunological response has been associated with the ability to resist the development of anemia following infection [21, 22], variation in the expression of extracellular matrix metalloproteinase [23], differences in chemokines, toll-like and chemokine receptors [24, 25], variation in genes that limit the supply of blood-meal to ticks and enhance both innate and adaptive immune responses [26]. Several genes have been proposed as candidates for tick resistance. Lumican (LUM) was identified as a potential biomarker for tick resistance in cattle [26]. [27] reported two SNPs in the ELTD1 gene that are associated with tick burdens in both dairy and beef cattle, and a haplotype of nine tag SNPs and two others that are associated with tick counts in dairy cattle. The authors reported also haplotypes spanning the ITGA11 gene that were associated with tick burdens. Class I antigens of the bovine major histocompatibility complex have been reported to be associated with tick loads [28, 29]. Alleles at the bovine lymphocyte antigen (BOLA-DRB3) have also been linked with resistance to Rhipicephalus (Boophilus) microplus infestation [30–32].

Domestic sheep (Ovis aries) are the second most abundant ruminant livestock after cattle [33] and are an important component of livestock enterprises in tick-endemic areas. Ticks transmit to sheep several pathogens, including viruses (tick-borne encephalitis virus, Thogoto virus, Louping-ill virus, Crimean-Congo haemorrhagic fever virus), bacteria (Mycoplasma ovis, Anaplasma ovis, Borrelia burgdorferi s.l., Francisella tularensis, Dermatophilus congolensis, Coxiella burnetti) and protozoa (Babesia spp. and Theileria spp.) [34–40]. Some of these pathogens cause important zoonotic diseases such as Crimean Congo Haemorrhagic fever, Q fever and Human granulocytic anaplasmosis [41] that result in dramatic human and animal health impacts including significant socio-cultural and economic losses. Despite reports on individual variations in tick burdens in sheep including both prevalence and infestation intensity [42–45], little is known regarding the genetic basis of this phenotype.

Taking advantage of observed individual natural variation in tick infestation in two indigenous sheep populations (Barbarine (B; fat-tailed) and Queue Fine de ľOuest (Q; thin-tailed)) from Tunisia [40], we investigated the likely genetic basis of the variation by analysing genome-wide 600K SNP genotype data that was generated using 170 individuals showing high- and low-infestation/resistance (HR/LR) to ticks. Our results suggest that the resistance could be the outcome of multigene associations with RAB6B, TF, SLCO2A1, and STXBP6 the likely potential candidates and STX5 as the possible prime candidate driving the genetic variation in tick infestation and thus resistance in sheep.

Population genetic structure

Genetic structure and relationship were investigated with principal component analysis (PCA) and ADMIXTURE tool. The first and second PCs of the PCA explained respectively, 7.63% and 6.04% of the total genetic variation (Figure 1). They separate the study individuals into four genetic groups, herein named G1, G2, G3 and G4 (Figure 1a). These four groups do not correspond to the resistance status to ticks (HR/LR) (Figure 1b), the sampling region (northeast, northwest, southeast) (Figure 1c) and breed (B/Q) (Figure 1d).

Following ADMIXTURE analysis, the lowest CV error was at K = 4 (Figure 1e) suggesting the presence of at least four genetic clusters in the dataset. These clusters correspond to the four genetic groups revealed by PCA and do not correspond to the resistance status to ticks, sampling region and breed. Therefore, the results from the two algorithms corroborate each other and for consistency, we also named the four genetic clusters as G1, G2, G3 and G4 (Figure 1f). The G1 cluster occurs in a few individuals of G2, G3 and G4, suggesting the introgression of G1 into these three clusters.

Genome variation

Descriptive statistics, that provide insights on the levels of genetic diversity and estimated for the overall dataset, the two breeds (B/Q), the two cohorts (HR/LR) and the four genetic groups (G1, G2, G3, G4), are shown in Table 1. The overall eestimates of genetic diversity represented by observed (H_O) and expected (H_E) heterozygosity averaged 0.3305 ± 0.0305 and 0.3470 ± 0.00006, respectively. The individual values for each statistic are in general higher than 0.3286 ± 0.0344. The mean length of ROH was 3.406 ± 0.4072, ranging between 3.032 ± 0.5213 Mb (G1) and 3.728 ± 1.004 Mb (G4). The mean value of F_ROH was 0.0549 ± 0.0842 and ranged between 0.0239 ± 0.0527 (G1) and 0.1612 ± 0.1029 (G3). The average value of F was 0.0476 ± 0.0881 with a range of -0.0866 ± 0.1339 (G4) and 0.0525 ± 0.0991 (HR). The LR cohort shows higher values of H_O and H_E compared to the HR cohort. At the opposite, the LR cohort shows lower values of ROH size, F_ROH and F. For the four genetic groups, G4 had the highest values of H_O, H_E and mean length of ROH while the lowest values were observed in G1, G2 and G1, respectively. The F_ROH was highest in G3 and lowest in G1. Except G1, the other three genetic groups show negative values of F.

We assessed the decay in linkage disequilibrium (LD; Figure 2a) and the trends in effective population size (N_E; Figure 2b) for the overall dataset, the two cohorts (HR/LR) and the four genetic groups (G1-G4). The overall pattern of decay in LD as a function of genomic distance is the same across all the classes of datasets analysed. It generally reveals higher LD at shorter distances which decays rapidly and plateaus off at around 0.2Mb. The G2, G3 and G4 groups show persistently higher r² values (r² > 0.15) and thus higher LD compared to the overall dataset, HR, LR and G1 (r² < 0.15). The trends in N_E over generation time show different profiles for the data classes analysed. The overall dataset and G1 group show the highest N_E which increases gradually reaching the maxima around 350 generations ago after which it declines rapidly to the present time. The HR and LR cohorts show the next lowest estimates of N_E which increase gradually up to around 350 generations ago followed by a rapid decline to the present. The G2, G3 and G4, have the lowest N_E across all generations which declines gradually to the present time.

Genome-wide selection signature analysis

We assessed selection signatures using runs of homozygosity (ROH), Logistic regression- genome-wide association analysis (LR-GWAS), cross-population extended haplotype homozygosity (XP-EHH) and genetic differentiation (F_ST) to gain insights on the regions of the genome correlating with resistance to ticks (HR verses LR). We also used the XP-EHH and F_ST algorithms to investigate selection signatures that may explain the divergence between the four genetic groups, G1-G4, revealed by PCA and ADMIXTURE algorithms.

The ROH approach identified 110 and 105 ROH regions that spanned 280 and 281 annotated genes in the HR and LR cohorts (Figure 3a, 3b; Supplementary Table S1, S2), respectively. The LR-GWAS identified 242 candidate regions spanning 561 genes (Figure 4a), and F_ST identified 79 candidate regions, showing genetic differentiation between the HR and LR cohorts (Figure 4b), that spanned 243 genes. The XP-EHH identified 76 candidate regions under positive selection spanning 187 genes (Figure 4c). These 397 candidate regions overlapped 991 genes (Supplementary Table S3, S4, S5) of which 104 annotated genes were found in candidate regions that were commonly identified by at least two approaches.

We used the ROH analysis as an approach for identifying cohort-specific selection signatures. We therefore compared the ROH results of HR and LR (Figure 4a; 4b). This identified 30 ROH regions, spanning 57 annotated and 20 unannotated (prefixed “LOC”) genes that were specific to the HR cohort (Supplementary Table S6). Considering these HR-specific ROH regions and the regions identified by LR-GWAS, XP-EHH and F_ST, there were 45 candidate regions that were simultaneously identified by at least two methods across 17 autosomes and spanning 160 genes of which, 56 remain unannotated (Table 2). Among the putative candidate genes identified in the 45 candidate regions were immune system function genes (CDC42EP1, CD164, CD180, CDH18, MYO), solute carriers (SLCO2A1, SLC26A3, SLC24A3, SLC22A8) and a chemokine receptor (CCIN). Of these genes, SLCO2A1 (OAR1), MYO (OAR8), SLC24A3 (OAR13), CD180 (OAR16) and CDH18 (OAR16) occurred closest to the most significant SNP identified by LR-GWAS in the respective candidate regions. Four genes (RAB6B (OAR1), TF (OAR1), SLCO2A1 (OAR1), STXBP6 (OAR18) and STX5 (OAR21)) that occurred in our candidate regions are of particular interest with regard to the objectives of this study. They have direct or indirect association with body weight traits in various livestock species and resistance to ticks in cattle. Worth noting is that STX5 occurred in an ROH region that was specific to the HR cohort (Table 3a) and in the most significant region that differentiated HR from LR (Figure 4b).

Functional enrichment analysis was tested in the pool of 104 putative genes, excluding 56 that remain uncharacterized, that were present in the 45 candidate regions that overlapped between at least one HR-specific ROH region and the candidate regions that were identified by LR-GWAS, F_ST and XP-EHH (Table 2). The analysis resulted in nine KEGG Pathways and four GO terms that were significantly enriched (Supplementary Table S7). The two top-most significant terms were Lysosome (oas04142) and microRNAs in cancer (oas05206). These functions are represented by clusters of genes (such as GGA1, CD164, and BCL2L11) having roles in innate immunity and disease-related inflammation [46, 47].

To determine genomic regions associated with divergence between G1, G2, G3, and G4, we contrasted these groups with each other using XP-EHH and F_ST. The candidate regions that were identified by the two methods following comparative analysis involving each group are shown in Supplementary Table S8. The selection regions detected by XP-EHH were generally broader (Figure 5) compared to those identified by F_ST (Figure 6) which are narrower and perhaps more accurate. This may be because XP-EHH is LD-based, and LD is expected to extend over longer distances in regions under recent selection [48]. To narrow down the regions that overlapped between the two methods, we considered the ones identified by F_ST as the reference candidate regions. In total, 10, 18, 18 and 26 regions that overlapped in a comparative analysis involving at least one group-pair were identified for G1, G2, G3 and G4, respectively (Supplementary Table S8). These regions spanned 16, 44, 25 and 86 genes, respectively. The most significant region for G1, G2, G3 and G4 were located on OAR13, 7, 16 and 27, respectively. These regions spanned respectively 4, 2, 1 and 5 genes.

We performed functional enrichment analysis with the pool of 11, 38, 23 and 76 genes, excluding the uncharacterized ones, that were identified in each of the four genetic groups (G1, G2, G3, G4). The functional term clusters that were significantly enriched for the cluster of genes for each genetic group are shown in Supplementary Table S9.

Ticks and TBI impact negatively ruminant livestock production worldwide. Therefore, the control and ultimate elimination of T-TBI should be prioritized to minimize their impacts not only on animal health and production but also on long-term human and environmental wellbeing. Large variations in susceptibility to T-TBI points to genome-wide variability underpinning inter-animal variation in both tick infestation and the pathogens they transmit. Here, we generated ovine 600K SNP genotype data, from 170 animals of two breeds of Tunisian indigenous sheep graze natural communal pastures with no history of specific anti-tick prevention measures. The sheep comprised of individuals showing high and low resistance (HR/LR) to tick infestation and piroplasm infection under natural conditions [40]. We analysed the data with four approaches, ROH, LR-GWAS, XP-EHH and F_ST, to investigate signatures of selection underpinning variation in tick infestation and therefore resistance. As noted recently by [49], we also acknowledge that using naturally infected animals runs the risk of resistant animals being a cocktail of truly highly resistant animals and those that were never exposed to infestations. These factors may dilute the certainty that animals showing high or low resistance to different pathogens, including viruses, bacteria and parasites, could in fact be functionally dissimilar [49].

Compared to the HR cohort, LR showed a comparatively higher level of genetic diversity and lower levels of inbreeding, although the differences were not significant (P > 0.05). The differences in the levels of inbreeding in the two cohorts, as determined from F_ROH and F, were also not significant. The variability in the four genetic groups revealed by PCA and ADMIXTURE was of the same magnitude as that of the two cohorts. Comparable levels of genetic variability have been reported in Egyptian [50], Algerian [51], Tunisian [52] and Russian [53] breeds of sheep.

The PCA and ADMIXTURE provided corroborating evidence suggesting absence of genetic stratification that is consistent with the a priori classification of the study individuals based on their susceptibility/resistance status to tick infestation, geographic sampling regions and breeds. The lack of genetic differentiation based on susceptibility/resistance to tick infestation is not surprising. A similar lack of genetic stratification that corresponds to a priori classification of sheep based on prolificacy [54] or resistance to gastro-intestinal nematodes [49, 55] has been reported. [56] also observed a lack of clear differentiation between sub-populations of dual-purpose black and white and German Holstein cattle. Our result suggests that the animals that were identified to be of high- and low-resistance to ticks are not highly divergent for this phenotype. This can be attributed to various reasons including, weak selection pressure driving differences in tick susceptibility/resistance and thus any favourable alleles are likely to be very rare, low level of tick burden that does not result in detectable genomic differentiation, and absence of farmer-driven preferential use of tick resistant animals for breeding. One or more of these reasons may explain the lack of genetic differentiation between the two cohorts. The lack of genetic differentiation corresponding to breeds and sampling regions was not unexpected. Past and recent human-mediated translocations and dispersal of the two breeds across Tunisia has brought the two breeds in close geographic proximity [57] and is resulting in cross-mating that is homogenizing their genomes. A similar finding was reported by [52]. The fact that the genomes of the B and Q breeds could not be differentiated by ADMIXTURE in each cluster provides evidence at the genome level that supports past and on-going uncontrolled mating of the fat-tailed Barbarine with thin-tailed breeds to reduce the fat content, especially of the tail, in the Barbarine carcass [58].

Linkage disequilibrium was estimated for all marker pairs using the metric r² plotted as a function of increasing genomic distance. The overall pattern of decay in LD is the same for all the classes of datasets analysed. It reveals higher LD at shorter distances which declines rapidly and plateaus off at around 200 Kb. This differs from what has been observed in commercial breeds whose LD plateaus off at around 150 Kb [59] and is most likely the result of differences in breeding history. Specifically, the study populations comprise a broad genetic base while that of commercial breeds has been narrowed down by bottlenecks associated with breed formation and stringent artificial selection for traits of economic importance. Both the HR and LR cohorts showed identical patterns of LD decay over increasing genomic physical distances, and in the trend in N_E over the past 1,000 generations. Although the two cohorts showed differences in the prevalence of tick infestation [40], this result suggests that they share aspects of their past and recent population demographic, genetic and breeding history. This may be the case as their genomes show no genetic differentiation that is commensurate with their a priori classification based on tick susceptibility/resistance. We also assessed the demographic profiles of the four genetic groups that underly the genome architecture of the study individuals as revealed by PCA and ADMIXTURE. The G2, G3 and G4 genetic groups had the highest r² values (r² > 0.15) and thus highest LD over short genomic distances and a decay curve with persistently elevated values compared to G1, HR and LR cohorts. A similar pattern was observed in the primitive Soay sheep [59] and was attributed to its small effective population size due to its genetic isolation in Scotland’s St. Kilda Island. Indeed, G2, G3 and G4 show the lowest N_E which declines gradually over all generations. The possible reason(s) for the observed decline in N_E in these three groups remains unknown and is worth further investigation. In contrast, LD was lower (r² < 0.15) in HR, LR and G1 across all genomic distance intervals and they also showed the highest N_E sizes over all generations. A similar pattern and magnitude of LD was reported in the Iranian Qezel sheep and was attributed to high genetic diversity in the breed [59]. This explanation may not suffice to explain our results as there was no significant differences in the level of genetic diversity across the different classes of datasets that we analysed. Despite being unexplained, we postulate that the low LD and high N_E in HR, LR and G1, and similar magnitude of genetic diversity is unlikely to reflect any biological phenomenon.

In line with what we observe in our study, the analysis of SNP arrays shows that LD decays faster and at much shorter genomic distances in sheep [59, 60]. While LD estimates are not comparable between studies due to differences in sample sizes, number and type of markers used, it is evident that sheep are characterized by lower LD than other domestic animals. For instance, the extent of LD in various cattle breeds [61] is greater than that found in our study. [62] reported higher r² values of between 0.19 and 0.26 in four pig breeds. This observation may be the result of less intensive selection in sheep for traits of economic importance and that the species is derived from a larger initial gene pool compared to other domestic animals [63]. Indeed, multiple male and female lineages and introgression from wild relatives have been implicated to constitute the modern-day diversity of domestic sheep [64–66].

Since millennia, ticks have parasitized animals for a blood-meal to the extent of developing, through co-evolution, a sophisticated armoury that guarantees their biological success. Here, to explore genomic signatures associated with individual variation in tick infestation and thus host resistance to ticks, we segregated the study individuals into two extreme groups comprising animals showing high- and low-resistance (HR/LR) to tick infestations. The HR group comprised individuals showing neither tick infestation nor piroplasm infection, and the LR group included animals that were infested by ticks and infected by piroplasms. The differences in tick infestation/resistance phenotype between the two groups was significant [40] and their comparison was used to maximize the likelihood of identifying biologically meaningful and statistically significant results. We therefore performed a comparative analysis with ROH, LR-GWAS, XP-EHH and F_ST approaches to detect selection signatures and SNPs that are likely associated with variation and thus resistance to tick infestation.

The four analytical approaches revealed 45 candidate regions, that overlapped between at least two approaches, which spanned 104 annotated genes. We observed a selection signature on OAR21 that was identified by LR-GWAS and F_ST which spanned nine genes including STX5 (Syntaxin-5). The ROH analysis showed that this region was specific to the HR cohort and it was the most significant in differentiating the HR from the LR cohort in the F_ST analysis. In a study of Belmont Red cattle, STX5 was amongst 11 out of 14 genes that showed a significant increase in its expression levels in the skin of animals that were highly resistant (HR) to ticks at time zero post infestation compared to animals of low resistance (LR) [67]. The expression level was more pronounced in 3-hour skin samples, suggesting a response to tick attachment and this could contribute to host innate immunity and to a higher resistance to ticks. STX5 regulates the endoplasmic reticulum channel-release properties of polycystin-2, a member of the transient receptor potential cation channel family that can function as an intracellular calcium (Ca²⁺) release channel [68]. Based on the increased expression levels of the 11 out of 14 genes they tested, most of which are Ca²⁺ dependent genes, [67] suggested that the high mRNA transcription levels of Ca²⁺ signaling genes in skin of HR animals accounts for their greater resistance to ticks. Thus, previous tick exposure may prime animals that exhibit the HR phenotype to resist further infestations via increased expression of Ca²⁺ signaling genes. Based on these results, we suggest that the STX5 could be the prime candidate gene driving tick resistance in sheep.

Environmental changes can exert either positive or negative effects on thermoregulation mechanisms that could influence tick burdens [69]. It has been observed that tick burdens in cattle might be correlated with traits that influence thermal comfort [70]. For instance, traits such as skin thickness, hair density and skin secretions (quantity and quality) could influence tick resistance in domestic livestock and might also influence thermal comfort since they affect the ability of the animal to dissipate heat [71]. Observations in Colombian cattle showed that high temperature humidity index (THI) values were associated with lower tick burdens and a higher tick infestation would be expected when animals experience higher thermal discomfort [70]. These findings are of relevance in our study because one of our candidate regions found on OAR18 that was revealed by LR-GWAS and XP-EHH, spanned STXBP6 (Syntaxin Binding Protein 6) and another candidate region on OAR1, that was revealed by LR-GWAS, XP-EHH and F_ST, spanned SLCO2A1 (OATP2A1). Both genes were closest to the most significant SNP as revealed by LR-GWAS. STXBP6 was one of the genes found to be upregulated in the testes of roosters exposed to acute heat stress [72]. Transcripts of STXBP6 were also found to be among eight other genes that were correlated with the modified Rodnan skin thickness score and forced vital capacity in humans suffering from scleroderma and systemic sclerosis [73], a condition that is characterized by thickening and hardening of the skin. One of the physical barriers that can confer resistance to tick infestation(s) include skin thickness, with animals with thin skins having significantly lower tick burdens [67, 74, 75]. Furthermore, except for birds and canids, thin skins allow animals to dissipate more heat through sweating and evaporative cooling when ambient temperatures are above the thermoneutral zones [69], and at the same time, it decreases tick attachment rate to the host’s skin. SLCO2A1, a prostaglandin transporter, maintains an increased interstitial concentration of PGE2, a major chemical mediator of febrile response, in the hypothalamus which plays a key role in thermoregulation [76]. The function of these genes is most likely complemented by the GNE, the gene closest to the top-most SNP as identified by LR-GWAS on a candidate region on OAR2, that likely plays a role in adaptation to climate-mediated selective pressures in sheep [77, 78]. Taken together, and the fact that Tunisia is a relatively hot country with annual average temperature of 29 ^oC (range 5.7 ^oC-37.1 ^oC; www.climate-data.org), this information leads us to hypothesize that the STXBP6 could have a potential pleiotropic effect on skin thickness and thermoregulation in sheep that enhances tick and heat stress tolerance. SLCO2A1 on the other hand, can enhance tick resistance by regulating fever during infestation with ticks and infections by TBI as well as in thermoregulation. However, these hypotheses need to be validated with more data on appropriate phenotypes such as skin characteristics.

Birth weight is the earliest available body weight trait with considerable impact on lamb survivability and growth performance traits [79]. Our analysis revealed a candidate region on OAR1 that was identified by LR-GWAS, XP-EHH and F_ST that spanned among others the RAB6B and TF genes. This region was observed, following genome-wide association analysis, to be significantly associated with birth weight in sheep and could be considered potential candidates for the trait [80, 81]. Growth is essentially associated with bone development and it was found that STXBP6 had potential pleiotropic effect on bone tissue and fecundity traits in chickens [82] and was found to be in a region under selection in broilers [83] and layers [84] suggesting that STXBP6 may influence body weight. In several cases, T-TBI can destabilize host growth and development. To counter against this destabilization, we suggest that natural selection may be acting on the regions spanning RAB6B, TF and STXBP6 to ensure growth and development stability early in life as an adaptive strategy that ensures survival in the context of high T-TBI challenge.

The differentiation of the study individuals into four genetic groups independent of tick infestation status, geographic sampling and breed was surprising. It was also surprising that the G1 group was characterized by a different pattern of LD decay and the trend in N_E compared to G2, G3 and G4, implicating a different genetic and demographic history and/or profile for G1. To gain further insights into the genetic differentiation between these four genetic groups, we performed selection signature analysis with XP-EHH and F_ST approaches. These two analyses revealed several candidate regions spanning a number of putative genes showing differentiation between the four genetic groups. However, based on the functions of the genes, we are unable to suggest the exact genetic mechanisms that are driving the differentiation of the four genetic groups given that the identified candidate genes are associated with multiple functions some of which overlap. The lack of phenotypic data, such as production, reproduction, adaptation and physical traits, that describe the individuals of the four genetic groups makes it difficult to interpret our findings. We therefore suggest that there could be some underlying but yet unknown factors that are driving the genome divergence of the study populations. To determine these factors, more data as well as sampling of other breeds found in Tunisia and neighbouring countries will be required.

In this study, we investigated the possible genetic basis of resistance to tick infestation in two populations of Tunisian indigenous sheep. We therefore examined selection signatures in an a priori defined two groups of animals presenting contrasting phenotypes of high- and low-resistance to ticks, based on tick counts and piroplasm infections. The two cohorts were characterized by similar levels of genetic variation and a fine scale genomic structure that could not be explained by tick resistance status, geographic sampling region and breeds, suggesting that other unknown factors could be driving the genomic divergence. Four methods of detecting selection signatures identified regions of the genome that were most likely associated with differences in tick infestation and thus resistance with our results suggesting that STX5 could be a prime candidate gene for tick resistance in sheep. We further hypothesize that STXBP6 and SLCO2A1 could be potential candidate genes for tick resistance in indigenous sheep which should be investigated further. The occurrence of RAB6B and TF in a candidate region that was significantly associated with body weight traits indicates that natural selection could be enhancing growth and development stability as an adaptive strategy to tick infestation. Although the sample size used in our study is relatively small, the high density of marker loci used increases the confidence of our findings. Nevertheless, [85] found that a rather limited value of n = 15 diploid individuals, equivalent to 30 gametes, is sufficient to reach reasonable power when used with high density markers. This research confirms the importance of studying indigenous livestock genotypes within their production environments to gather information on genomic footprints of adaptive significance. Further studies on understanding the genetics of tick resistance can benefit from the inclusion of variables such as skin and coat characteristics in order to analyze their influence on tick infestation.

Study animals, sampling and genotyping

The blood samples used in this study were collected by veterinary surgeons, during routine animal health surveillance in farmers herds/flocks, using standard veterinary procedures in full respect of animal welfare and minimizing stress as per the guidelines for the care and use of animals of the Tunisian National Council of Veterinary Surgeons (TNCVS). The blood sampling was also done with the verbal permission, and in the presence, of animal owners. Therefore, permission from the Ethics Committee of the TNCVS was not required.

Two breeds, Barbarine (B) and Queue Fine de l’Ouest (Q), of Tunisian indigenous sheep were sampled. The animals graze predominantly on natural communal pastures year-round except during summer when they are grazed on crop residues. Prophylactic treatment is rare, but vaccinations against Brucellosis, Bluetongue, Foot and Mouth disease and Sheep-pox are done annually by the National Veterinary Services. The animals sampled for this study were from the northeast, northwest and southeast regions of Tunisia and were sampled as part of a larger cross-sectional study on T-TBI that was carried out between 2018 and 2020 [40]. The northeast and southeast regions are the homelands of the B breed although it has been dispersed throughout the country. It forms the foundation of the “Tunis” and “Barbaresca” breeds found respectively, in the USA and Italy [86]. The B is also found in Libya and eastern Algeria [87]. The Q breed predominates in the northwest region although they have been dispersed to the central, eastern and western plateaus of Tunisia. The breed is also found in mixed sheep-cropping systems in eastern Morocco and is genetically close to the Algerian Ouled Jellal breed [88].

For the cross-sectional study, 439 mature ewes were tagged and monitored through eight rounds of sampling. From each animal, whole blood (5 mls) was collected from each animal in EDTA coated vacutainers via jugular venipuncture and ticks were also collected from both ears, identified, counted and preserved in 70% ethanol. From the 439 ewes, 170 animals (B = 108; Q = 62) were purposely selected for this study based on tick scores and piroplasm infections. Genomic DNA was extracted from the blood samples with the Rapid Blood Genomic DNA Extraction Kit (Bio Basic, Markham, Canada). Tick scores, which were used to determine tick infestations, were visually estimated based on a scale from 0 (no ticks) to 5 (> 150 ticks) [89] while piroplasms were detected with the polymerase chain reaction technique [40].

The 170 samples were assigned to two cohorts based on extreme values of tick load and piroplasm infections viz: i) high tick-resistant (HR) - animals with no tick infestation(s) and/or no piroplasm infection(s) and ii) low tick-resistant (LR) - animals that are highly infested with ticks and/or infected by piroplasms (Table 1). The DNA from the 170 samples were genotyped with the Illumina 600K SNP BeadChip (Illumina Inc., San Diego, CA, USA) at Neogen GeneSeek, Lincoln NE, USA. The BeadChip comprises 606,006 probes targeting genome-wide SNPs, of which 577,401 are autosomal, 27,314 are on the X chromosome and 1,291 remain unassigned [90].

The raw genotypes were assessed for quality with PLINK1.9 [91] and then pruned with the following criteria: (i) one individual from a pair of highly related animals was discarded if they had an identity-by-state (IBS) score greater than 0.99, (ii) SNPs with minor allele frequencies (MAF) of no less than 0.01 were retained, (iii) individuals and SNPs with call rates lower than 90% and 95%, respectively were discarded and (iv) all unmapped SNPs were discarded. This generated a dataset of 537,214 SNPs comprising 74 HR (B = 43; Q = 31) and 96 LR (B = 65; Q = 31) individuals. This dataset was subjected to LD pruning with the parameters 50 5 0.5 which represent window size, step size and r² threshold, respectively resulting in 338,180 SNPs that were used for genetic structure and phylogenetic analysis.

Assessment of population genetic structure

Although the study individuals were classified a priori into two extreme cohorts, HR and LR, we first assessed population genetic structure and divergence to determine whether the underlying genome architecture corresponds to the two cohorts. We therefore performed principal component analysis (PCA) with PLINK v1.9 and the first two PC’s were plotted to visualize genetic relationships. We also inferred fine-scale genome structure and shared genome ancestry with the unsupervised mode of ADMIXTURE tool v1.30 [92], independent of any background information on the number and frequency of alleles in the ancestral gene pool. The ADMIXTURE tool was run with a K range of 1–8 clusters and we performed five replicate runs for each K to generate cross-validation (CV) errors. The lowest CV error was used to infer the optimal number of distinct genetic clusters in the dataset.

Assessment of genetic diversity and demographic dynamics

Expected (H_E) and observed (H_O) heterozygosity, patterns of LD decay across the genome and effective population size (N_E) across generation time were investigated for each breed, the HR and LR cohorts and the genetic groups revealed by PCA and ADMIXTURE algorithms. The H_E and H_O were calculated with PLINK v1.9. The extent of LD between pairwise SNPs was evaluated with the r² statistic calculated with PLINK v1.9. The N_E over generation time was estimated with the equation N_Et = (1/4c) (1/r² − 1) (Sved 1971), where N_Et is the N_E t generations ago (t = 1/2c); r² is the LD between pairwise SNPs; and c is the genetic distance in Morgans between pairs of SNPs.

Two coefficients of inbreeding were calculated with PLINK v1.9; the SNP based “F” and the runs of homozygosity (ROH) based “F_ROH”. The latter was computed as the ratio of the total length of ROH to the length of autosomes (2.45 Gb) [93]. The ROHs were identified for each animal with the following parameters: (i) the minimum number of SNPs in a sliding window was set to 50; (ii) the minimum ROH length was set to 1 Mb to eliminate the impact of strong LD; (iii) each ROH spanned a minimum of 80 consecutive SNPs; (iv) one heterozygous and five missing calls per window were allowed to avoid false negative results that can arise due to genotyping errors or missing genotypes; (v) the minimum SNP density was set to one SNP/100 kb, and (vi) the maximum gap between consecutive SNPs was set to 1 Mb.

Detection of selection signatures and association analysis

To detect selection signatures spanning regions associated with genomic targets for resistance to ticks, we first investigated the distribution of ROH across the genomes of the HR and LR cohorts. For each cohort, the frequency of a SNP within an ROH region was determined and a Manhattan plot visualising all the tested SNPs against their autosomal positions was generated. The most frequent SNPs occurring above the 50% cut-off threshold of the empirical distribution were taken as the most significant loci underlying ROH streatches that were possibly under selection. To identify the ROH streatches that are likely associated with tick resistance, we contrasted the HR and LR ROH regions and identified the ones that were specific to the HR cohort.

We used the population differentiation statistic, F_ST, to investigate regions of the genome showing divergence between the HR and LR cohorts. The unbiased pairwise F_ST [94] was computed for each SNP with the R package “HIERFSTAT” [95] using a window size of 100 Kb and a window-step size of 10 Kb. Windows with less than five SNPs were excluded from the analysis. The F_ST values were standardized by Z transformation following [96]. The windows that occurred within the top 0.001% of the F_ST values in each chromosome were considered the putative candidates under divergent selection.

We used the software developed by [97] to estimate the unstandardized XP-EHH statistics for all SNPs, after quality control, following the comparative analysis between the HR and LR cohorts. The unstandardized XP-EHH statistics were standardized using their means and variances. We estimated the p-values of the SNPs using the standard normal distribution following findings of previous studies [98–100]. We determined the candidate regions under positive selection by clustering the significant core SNPs (P-value < 0.05) within a distance of less than 200 kb from the top-most SNP.

We performed the logistic regression (LR) GWAS with PLINK v1.9 to explore further, the possible genomic regions and SNPs associated with variation in tick resistance. The HR and LR cohorts were used as the test and control groups, respectively. 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 age, sampling region and breed as covariates. The generated p-values were Bonferroni corrected to minimize false positive results among the SNPs that were likely under selection. The corrected p-values were standardized and the -Log10 (p-value) of 4.25 (the top 0.001) was set as the cut-off threshold to identify the candidate regions and associated SNPs. The estimations were summarized in 100 Kb window sizes and the genes and the top-most SNP found within the candidate regions were identified. The Manhattan plots showing the genome-wide distribution of the SNPs were generated with the R package “qqman” v3.5.1.

To determine the genomic regions associated with the divergence of the genetic groups revealed by PCA and ADMIXTURE algorithms, we used the XP-EHH and F_ST approaches. The former detects signatures arising from recent selection or those that are approaching fixation, and the latter identifies regions that show allele frequency differentiation which may be the result of long-term selection. The analyses were performed by contrasting one group with the others for example, G1 vs G2, G3, G4 etc. The approaches used to perform XP-EHH and F_ST while contrasting the HR and LR cohorts were also used in this analysis.

Functional annotation of candidate regions

The candidate regions revealed by ROH were analysed and the ones that were specific to the HR cohort identified. These HR-specific ROH regions and the regions identified by LR-GWAS, XP-EHH and F_ST were analysed, and those that overlapped between at least two methods were identified and merged using Bedtools v.2.28.0 [101]. The same analysis was performed to determine the regions that overlapped between the ones identified by XP-EHH and F_ST from the comparative analysis involving at least two genetic groups that were identified by the clustering algorithms. The genes spanned by the overlapping candidate regions were retrieved using the Biomart/Ensembl (http://www.ensembl.org/biomart) tool based on the Ovine v3.1 reference genome assembly. These set of genes were assessed for biological enrichment gene ontology (GO) and KEGG Pathway (www.kegg.jp/kegg/kegg1.html) terms compared to the full list of O. aries autosomal protein-coding genes with the functional annotation tool in DAVID v6.8 [102]. We further determined gene functions using the NCBI database (http://www.ncbi.nlm.nih.gov/gene/) and review of literature.

B = Barbarine; Q = Queue Fine de l’Ouest; HR = High resistant; LR = Low resistant; PCA = Principal component analysis; T-TBI = Ticks and tick-borne infections; ROH = runs of homozygosity; LR-GWAS = Logistic regression genome-wide association analysis; XP-EHH = cross-population extended haplotype homozygosity; F_ST = genetic differentiation statistic; OAR = Ovine chromosome

Ethics approval and consent to participate

The animals used in this study belonged to farmers who provided their verbal informed consent to their sampling which was performed as a routine veterinary monitoring activity (blood sampling and tick collection) by licensed veterinarians following the guidelines for the care and use of animals of the Tunisian National Council of Veterinary Surgeons; no ethical approvals were required.

Consent for publication

Not applicable

Availability of data and materials

The data used in this study is available from the communicating author on reasonable request upon signing a material transfer agreement.

Competing interests

The authors declare no competing interests in all aspects of the study

Funding

We thank all donors and organizations that globally support the work of the CGIAR Research Program on Livestock through their contributions to the CGIAR Trust Fund. This research was partly funded by the Arab Fund for Social and Economic Development received through the International Center for Agricultural Research in the Dr Areas (Grant number 131,001) and the Tunisian Ministry of Higher Education and Scientific Research through the “Laboratorie ďépidémiologie ďinfections enzootiques des herbivores, application à la lute” (Grant number LR16AGR01).

Author’s contributions

GM, RM, HA, RB and MJM conceived the study. KKM, RH and SL carried out the field work including sampling, DNA extraction and PCR under supervision of GM and RM. AAM analysed the data under supervision from MJM and both wrote the manuscript. All authors read and approved the final version of the manuscript.

Acknowledgements

We acknowledge the valuable contribution of veterinarians that facilitated sheep access and sheep owners for volunteering their flocks for the study.

Wikel SK. Host immunity to ticks. Annu Rev Entomol. 1996;41:1–22.
Wikel SK. Modulation of the host immune system by ectoparasitic arthropods: Blood-feeding and tissue-dwelling arthropods manipulate host defenses to their advantage. Bioscience. 1999;49:311–20.
Mans BJ, Neitz AWH. Adaptation of ticks to a blood-feeding environment: Evolution from a functional perspective. Insect Biochem Mol Biol. 2004;34:1–17.
Francischetti IMB, Sa-Nunes A, Mans BJ, Santos IM, Ribeiro JMC. The role of saliva in tick feeding. Front Biosci. 2009;14:2051–88. doi:10.2741/3363.
Oliveira CJF, Carvalho WA, Garcia GR, Gutierrez FRS, de Miranda Santos IKF, Silva JS, et al. Tick saliva induces regulatory dendritic cells: MAP-kinases and Toll-like receptor-2 expression as potential targets. Vet Parasitol. 2010;167:288–97.
Jongejan F, Uilenberg G. The global importance of ticks. Parasitology. 2004;129 SUPPL.
de la Fuente J, Kopáček P, Lew-Tabor A, Maritz-Olivier C. Strategies for new and improved vaccines against ticks and tick-borne diseases. Parasite Immunol. 2016;38:754–69.
Mendes EC, Mendes MC, Sato ME. Diagnosis of amitraz resistance in Brazilian populations of Rhipicephalus (Boophilus) microplus (Acari: Ixodidae) with larval immersion test. Experimental and Applied Acarology. 2013;61.
Abbas RZ, Zaman MA, Colwell DD, Gilleard J, Iqbal Z. Acaricide resistance in cattle ticks and approaches to its management: The state of play. Veterinary Parasitology. 2014;203.
Willadsen P. Tick control: Thoughts on a research agenda. Vet Parasitol. 2006;138:161–8.
Regitano LCA, Ibelli AMG, Gasparin G, Miyata M, Azevedo ALS, Coutinho LL, et al. On the Search for markers of tick resistance in bovines. In: Dev. Biol. 2008. p. 225–30.
Machado MA, S Azevedo AL, Teodoro RL, Pires MA, CD Peixoto MG, de Freitas C, et al. Genome wide scan for quantitative trait loci affecting tick resistance in cattle (Bos taurus × Bos indicus). BMC Genomics. 2010;11.
Domingos A, Antunes S, Borges L, do Rosário VE. Approaches towards tick and tick-borne diseases control. Revista da Sociedade Brasileira de Medicina Tropical. 2013;46:265–9.
de La Fuente J, Contreras M. Tick vaccines: Current status and future directions. Expert. Rev. Vaccines. 2015;14:1367–76.
Schetters T, Bishop R, Crampton M, Kopáček P, Lew-Tabor A, Maritz-Olivier C, et al. Cattle tick vaccine researchers join forces in CATVAC. In: Parasit. Vectors. 2016.
Kaiser MN, Sutherst RW, Bourne AS. Relationship between ticks and zebu cattle in southern Uganda. Trop Anim Health Prod. 1982;14.
de Castro JJ, Young AS, Dransfield RD, Cunningham MP, Dolan TT. Effects of tick infestation on Boran (Bos indicus) cattle immunised against theileriosis in an endemic area of Kenya. Research in veterinary science. 1985;39.
Tatchell RJ, Chimwani D, Chirchir SJ, Ong’are JO, Mwangi E, Rinkanya F, et al. A study of the justification for intensive tick control in Kenyan rangelands. The Veterinary record. 1986;119.
Frisch JE. Towards a permanent solution for controlling cattle ticks. Int J Parasitol. 1999;29:57–71.
Mapholi NO, Marufu MC, Maiwashe A, Banga CB, Muchenje V, MacNeil MD, et al. Towards a genomics approach to tick (Acari: Ixodidae) control in cattle: A review. Ticks. Tick. Borne. Dis. 2014;5:475–83.
Murray M, Trail JC, D’Ieteren GD. Trypanotolerance in cattle and prospects for the control of trypanosomiasis by selective breeding. Revue scientifique et technique (International Office of Epizootics). 1990;9:369–86.
PALING RW, MOLOO SK, SCOTT JR, GETTINBY G, MCODIMBA FA, MURRAY M. Susceptibility of N’Dama and Boran cattle to sequential challenges with tsetse-transmitted clones of Trypanosoma congolense. Parasite Immunol. 1991;13:427–45.
Wang YH, Reverter A, Kemp D, McWilliam SM, Ingham A, Davis CA, et al. Gene expression profiling of Hereford Shorthorn cattle following challenge with Boophilus microplus tick larvae. Aust J Agric Res. 2007;47:1397–407.
Piper EK, Jonsson NN, Gondro C, Lew-Tabor AE, Moolhuijzen P, Vance ME, et al. Immunological profiles of Bos taurus and Bos indicus cattle infested with the cattle tick, Rhipicephalus (Boophilus) microplus. Clin Vaccine Immunol. 2009;16:1074–86.
Piper EK, Jackson LA, Bagnall NH, Kongsuwan KK, Lew AE, Jonsson NN. Gene expression in the skin of Bos taurus and Bos indicus cattle infested with the cattle tick, Rhipicephalus (Boophilus) microplus. Vet Immunol Immunopathol. 2008;126:110–9.
Marima JK, Nel CL, Marufu MC, Jonsson NN, Dube B, Dzama K. A genetic and immunological comparison of tick-resistance in beef cattle following artificial infestation with Rhipicephalus ticks. Exp Appl Acarol. 2020;80:569–90.
Porto Neto LR, Bunch RJ, Harrison BE, Prayaga KC, Barendse W. Haplotypes that include the integrin alpha 11 gene are associated with tick burden in cattle. BMC Genetics. 2010;11.
Stear MJ, Nicholas FW, Brown SC, Holroyd RG. Class I antigens of the bovine major histocompatibility system and resistance to the cattle tick (Boophilus microplus) assessed in three different seasons. Vet Parasitol. 1989;31:303–15.
Stear MJ, Hetzel DJS, Brown SC, Gershwin LJ, Mackinnon MJ, Nicholas FW. The relationships among ecto- and endoparasite levels, class I antigens of the bovine major histocompatibility system, immunoglobulin E levels and weight gain. Vet Parasitol. 1990;34:303–21.
Martinez ML, Machado MA, Nascimento CS, Silva MVGB, Teodoro RL, Furlong J, et al. Association of BoLA-DRB3.2 alleles with tick (Boophilus microplus) resistance in cattle. Genet Mol Res. 2006;5:513–24.
Acosta-Rodríguez R, Alonso-Morales R, Balladares S, Flores-Aguilar H, García-Vazquez Z, Gorodezky C. Analysis of BoLA class II microsatellites in cattle infested with Boophilus microplus ticks: Class II is probably associated with susceptibility. Vet Parasitol. 2005;127:313–21.
Untalan PM, Pruett JH, Steelman CD. Association of the bovine leukocyte antigen major histocompatibility complex class II DRB3*4401 allele with host resistance to the Lone Star tick, Amblyomma americanum. Vet Parasitol. 2007;145:190–5.
Gilbert JA, Blaser MJ, Caporaso JG, Jansson JK, Lynch S v., Knight R. Current understanding of the human microbiome. Nat Med. 2018;24:392–400.
ben Said M, Belkahia H, Karaoud M, Bousrih M, Yahiaoui M, Daaloul-Jedidi M, et al. First molecular survey of Anaplasma bovis in small ruminants from Tunisia. Vet Microbiol. 2015;179:322–6.
Said M ben, Belkahia H, Alberti A, Abdi K, Zhioua M, Daaloul-Jedidi M, et al. First molecular evidence of Borrelia burgdorferi sensu lato in goats, sheep, cattle and camels in Tunisia. Ann Agric Environ Med. 2016;23:442–7.
Gharbi M. Epidemiological Study of Sheep Anaplasmosis (Anaplasma ovis Infection) in Kairouan, Central Tunisia. Adv Parasitol. 2015;2:30–4.
Rjeibi MR, Gharbi M, Mhadhbi M, Mabrouk W, Ayari B, Nasfi I, et al. Prevalence of piroplasms in small ruminants in North-West Tunisia and the first genetic characterisation of Babesia ovis in Africa. Parasite. 2014;21.
Rjeibi MR, Darghouth MA, Rekik M, Amor B, Sassi L, Gharbi M. First Molecular Identification and Genetic Characterization of Theileria lestoquardi in Sheep of the Maghreb Region. Transbound Emerg Dis. 2016;63:278–84.
Rjeibi MR, Darghouth MA, Gharbi M. Prevalence of Theileria and Babesia species in Tunisian sheep. Onderstepoort J Vet Res. 2016;83.
Khbou MK, Rouatbi M, Romdhane R, Sassi L, Jdidi M, Haile A, et al. Tick infestation and piroplasm infection in barbarine and queue fine de l’ouest autochthonous sheep breeds in Tunisia, North Africa. Animals. 2021;11:1–17.
Alessandra T, Santo C. Tick-borne diseases in sheep and goats: Clinical and diagnostic aspects. Small Rumin Res. 2012;106 SUPPL.
Ogore PB, Baker RL, Kenyanjui M, Thorpe W. Assessment of natural ixodid tick infestations in sheep. Small Rumin Res. 1999;33:103–7.
Granquist EG, Stuen S, Crosby L, Lundgren AM, Alleman AR, Barbet AF. Variant-specific and diminishing immune responses towards the highly variable MSP2(P44) outer membrane protein of Anaplasma phagocytophilum during persistent infection in lambs. Vet Immunol Immunopathol. 2010;133:117–24.
Stuen S, Grøva L, Granquist EG, Sandstedt K, Olesen I, Steinshamn H. A comparative study of clinical manifestations, haematological and serological responses after experimental infection with Anaplasma phagocytophilum in two Norwegian sheep breeds. Acta Vet Scand. 2011;53.
Grøva L. Tick-borne fever in sheep – production loss and preventive measures. Dept of Aquaculture and Animal Science. 2011;Ph.D.
Moumène A, Meyer DF. Ehrlichia’s molecular tricks to manipulate their host cells. Microbes and Infection. 2016;18:172–9.
Chen J, Ao L, Yang J. Long non-coding RNAs in diseases related to inflammation and immunity. Ann Transl Med. 2019;7:494–494.
Fu W, Dekkers JCM, Lee WR, Abasht B. Linkage disequilibrium in crossbred and pure line chickens. Genet Sel Evol. 2015;47.
Ahbara AM, Rouatbi M, Gharbi M, Rekik M, Haile A, Rischkowsky B, et al. Genome-wide insights on gastrointestinal nematode resistance in autochthonous Tunisian sheep. Sci Rep. 2021;11.
Kim ES, Elbeltagy AR, Aboul-Naga AM, Rischkowsky B, Sayre B, Mwacharo JM, et al. Multiple genomic signatures of selection in goats and sheep indigenous to a hot arid environment. Heredity. 2016;116:255–64.
Gaouar SBS, Lafri M, Djaout A, El-Bouyahiaoui R, Bouri A, Bouchatal A, et al. Genome-wide analysis highlights genetic dilution in Algerian sheep. Heredity. 2017;118:293–301.
Bedhiaf-Romdhani S, Baazaoui I, Ciani E, Mastrangelo S, Sassi M ben. Genetic structure of Tunisian sheep breeds as inferred from genome-wide SNP markers. Small Ruminant Research. 2020;191.
Deniskova TE, Dotsev A v., Selionova MI, Kunz E, Medugorac I, Reyer H, et al. Population structure and genetic diversity of 25 Russian sheep breeds based on whole-genome genotyping. Genet Sel Evol. 2018;50.
Dolebo AT, Khayatzadeh N, Melesse A, Wragg D, Rekik M, Haile A, et al. Genome-wide scans identify known and novel regions associated with prolificacy and reproduction traits in a sub-Saharan African indigenous sheep (Ovis aries). Mamm Genome. 2019;30:339–52.
Berton MP, Oliveira Silva RM, Peripolli E, Stafuzza NB, Martin JF, Álvarez MS, et al. Genomic regions and pathways associated with gastrointestinal parasites resistance in Santa Inês breed adapted to tropical climate. J Anim Sci Biotechnol. 2017;8.
Naderi S, Moradi MH, Farhadian M, Yin T, Jaeger M, Scheper C, et al. Assessing selection signatures within and between selected lines of dual-purpose black and white and German Holstein cattle. Animal Genetics. 2020;51.
Khaldi Z, Souid S, Rekik B, Haddad B. Genetic variability and structure of three tunisian ovine breeds based on microsatellite polymorphism. Indian J Anim Res. 2020;54:1459–64.
Bedhiaf-Romdhani S, Djemali M, Zaklouta M, Iniguez L. Monitoring crossbreeding trends in native Tunisian sheep breeds. Small Rumin Res. 2008;74:274–8.
Kijas JW, Porto-Neto L, Dominik S, Reverter A, Bunch R, McCulloch R, et al. Linkage disequilibrium over short physical distances measured in sheep using a high-density SNP chip. Anim Genet. 2014;45:754–7.
Al-Mamun HA, A Clark S, Kwan P, Gondro C. Genome-wide linkage disequilibrium and genetic diversity in five populations of Australian domestic sheep. Genet Sel Evol. 2015;47.
Porto-Neto LR, Kijas JW, Reverter A. The extent of linkage disequilibrium in beef cattle breeds using high-density SNP genotypes. Genet Sel Evol. 2014;46.
Badke YM, Bates RO, Ernst CW, Schwab C, Steibel JP. Estimation of linkage disequilibrium in four US pig breeds. BMC Genomics. 2012;13.
Kijas JW, Lenstra JA, Hayes B, Boitard S, Neto LR, Cristobal MS, et al. Genome-wide analysis of the world’s sheep breeds reveals high levels of historic mixture and strong recent selection. PLoS Biol. 2012;10.
Pedrosa S, Uzun M, Arranz JJ, Gutiérrez-Gil B, San Primitivo F, Bayón Y. Evidence of three maternal lineages in near eastern sheep supporting multiple domestication events. Proc R Soc Biol. 2005;272:2211–7.
Cao YH, Xu SS, Shen M, Chen ZH, Gao L, Lv FH, et al. Historical Introgression from Wild Relatives Enhanced Climatic Adaptation and Resistance to Pneumonia in Sheep. Mol Biol Evol. 2021;38:838–55.
Deng J, Xie XL, Wang DF, Zhao C, Lv FH, Li X, et al. Paternal Origins and Migratory Episodes of Domestic Sheep. Curr Biol. 2020;30:4085–4095.e6.
Tabor AE, Ali A, Rehman G, Garcia GR, Zangirolamo AF, Malardo T, et al. Cattle Tick Rhipicephalus microplus-host interface: A review of resistant and susceptible host responses. Front. Cell. Infect. Microbiol. 2017;7 DEC.
Geng L, Boehmerle W, Maeda Y, Okuhara DY, Tian X, Yu Z, et al. Syntaxin 5 regulates the endoplasmic reticulum channel-release properties of polycystin-2. Proc Natl Acad Sci USA. 2008;105:15920–5.
Katiyatiya CLF, Muchenje V, Mushunje A. Seasonal variation in coat characteristics, tick loads, cortisol levels, some physiological parameters and temperature humidity index on Nguni cows raised in low- and high-input farms. Int J Biometeorol. 2015;59:733–43.
Rocha JF, Martínez R, López-Villalobos N, Morris ST. Tick burden in Bos taurus cattle and its relationship with heat stress in three agroecological zones in the tropics of Colombia. Parasit Vectors. 2019;12.
Shyma KP, Gupta JP, Singh V. Breeding strategies for tick resistance in tropical cattle: a sustainable approach for tick control. J Parasit Dis. 2015;39:1–6.
Wang SH, Cheng CY, Tang PC, Chen CF, Chen HH, Lee YP, et al. Differential gene expressions in testes of L2 strain Taiwan country chicken in response to acute heat stress. Theriogenology. 2013;79.
Assassi S, Wu M, Tan FK, Chang J, Graham TA, Furst DE, et al. Skin gene expression correlates of severity of interstitial lung disease in systemic sclerosis. Arthritis Rheum. 2013;65:2917–27.
de Castro JJ, Capstick PB, Nokoe S, Kiara H, Rinkanya F, Slade R, et al. Towards the selection of cattle for tick resistance in Africa. Exp Appl Acarol. 1991;12:219–27.
Jonsson NN, Piper EK, Constantinoiu CC. Host resistance in cattle to infestation with the cattle tick Rhipicephalus microplus. Parasite Immunol. 2014;36:553–9.
Nakamura Y, Nakanishi T, Shimada H, Shimizu J, Aotani R, Maruyama S, et al. Prostaglandin transporter OATP2A1/SLCO2A1 is essential for body temperature regulation during fever. J Neurosci. 2018;38:5584–95.
Lv FH, Agha S, Kantanen J, Colli L, Stucki S, Kijas JW, et al. Adaptations to climate-mediated selective pressures in sheep. Mol Biol Evol. 2014;31:3324–43.
Mastrangelo S, Tolone M, Sardina MT, Sottile G, Sutera AM, di Gerlando R, et al. Genome-wide scan for runs of homozygosity identifies potential candidate genes associated with local adaptation in Valle del Belice sheep. Genet Sel Evol. 2017;49.
Ptáček M, Ducháček J, Hakl J, Fantová M. of multivariate relations among birth weight, survivability traits, growth performance, and some importAnalysisant factors in Suffolk lambs. Arch Anim Breed. 2017;60:43–50.
Zhang L, Liu J, Zhao F, Ren H, Xu L, Lu J, et al. Genome-Wide Association Studies for Growth and Meat Production Traits in Sheep. PLoS One. 2013;8.
Ghasemi M, Zamani P, Vatankhah M, Abdoli R. Genome-wide association study of birth weight in sheep. Animal. 2019;13:1797–803.
Johnsson M, Rubin CJ, Höglund A, Sahlqvist AS, Jonsson KB, Kerje S, et al. The role of pleiotropy and linkage in genes affecting a sexual ornament and bone allocation in the chicken. Mol Ecol. 2014;23:2275–86.
Fu W, Lee WR, Abasht B. Detection of genomic signatures of recent selection in commercial broiler chickens. BMC Genetics. 2016;17.
Qanbari S, Strom TM, Haberer G, Weigend S, Gheyas AA, Turner F, et al. A High Resolution Genome-Wide Scan for Significant Selective Sweeps: An Application to Pooled Sequence Data in Laying Chickens. PLoS One. 2012;7.
Ma Y, Ding X, Qanbari S, Weigend S, Zhang Q, Simianer H. Properties of different selection signature statistics and a new strategy for combining them. Heredity. 2015;115:426–36.
Djemali M, Aloulou R, Sassi M ben. Adjustment factors and genetic and phenotypic parameters for growth traits of Barbarine lambs in Tunisia. Small Rumin Res. 1994;13:41–7.
ben Salem H, Lassoued N, Rekik M. Merits of the fat-tailed Barbarine sheep raised in different production systems in Tunisia: Digestive, productive and reproductive characteristics. Trop Anim Health Prod. 2011;43:1357–70.
Rekik M, Aloulou R, ben Hammouda M. Small Ruminant Breeds of Tunisia. In: Iniguez L, editor. Characterisation of Small Ruminant Breeds in West Asia and North Africa. Beirut: International Center for Agricultural Research in the Dry Areas (ICARDA); 2005. p. 91–140.
Prayaga KC, Corbet NJ, Johnston DJ, Wolcott ML, Fordyce G, Burrow HM. Genetics of adaptive traits in heifers and their relationship to growth, pubertal and carcass traits in two tropical beef cattle genotypes. Animal Production Science. 2009;49.
Anderson R, McEwan J, Brauning R, Kijas J, Dalrymple J, Worley K, et al. Development of a high density (600K) Illumina Ovine SNP chip and its use to fine map the yellow fat locus.. In: January 2014 _ Plant and Animal Genome XXII Conference. San Diego, California USA: ISGC; 2014.
Purcell S, Neale B, Todd-Brown K, Thomas L, Ferreira MAR, Bender D, et al. PLINK: A tool set for whole-genome association and population-based linkage analyses. Am J Hum Genet. 2007;81:559–75.
Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655–64.
McQuillan R, Leutenegger AL, Abdel-Rahman R, Franklin CS, Pericic M, Barac-Lauc L, et al. Runs of Homozygosity in European Populations. Am J Hum Genet. 2008;83:359–72.
Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.
Goudet J. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Mol Ecol Notes. 2005;5:184–6.
Ahbara A, Bahbahani H, Almathen F, Abri M al, Agoub MO, Abeba A, et al. Genome-wide variation, candidate regions and genes associated with fat deposition and tail morphology in Ethiopian indigenous sheep. Front Genet. 2019;10 JAN.
Pickrell JK, Coop G, Novembre J, Kudaravalli S, Li JZ, Absher D, et al. Signals of recent positive selection in a worldwide sample of human populations. Genome Res. 2009;19:826–37.
Sabeti PC, Varilly P, Fry B, Lohmueller J, Hostetter E, Cotsapas C, et al. Genome-wide detection and characterization of positive selection in human populations. Nature. 2007;449:913–8.
Ma. Y, Zhang. H, Zhang Q, Ding X. Identification of selection footprints on the X Chromosome in pig. PLoS One. 2014;9.
Zhao F ping, Wei C hong, Zhang L, Liu J sen, Wang G kai, Zeng T, et al. A genome scan of recent positive selection signatures in three sheep populations. J Integr Agric. 2016;15:162–74.
Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.
Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

Table 1

Genetic diversity indices generated for the breeds, cohorts and genetic groups of Tunisian Sheep used in this study

	N	H_O (Mean ± Sd)	H_E (Mean ± Sd)	ROH (Mb) (Mean ± Sd)	F_ROH (Mean ± Sd)	F (Mean ± Sd)
Breed
Barbarine (B)	104	0.3313 ± 0.0299	0.3470 ± 0.00003	3.332 ± 0.4434	0.0515 ± 0.0826	0.0452 ± 0.0862
Queue Fine de l’Ouest (Q)	60	0.3291 ± 0.0322	0.3468 ± 0.00010	3.520 ± 0.5242	0.0625 ± 0.0881	0.0511 ± 0.0929
Cohort
LR	77	0.3321 ± 0.0289	0.3477 ± 0.00007	3.377 ± 0.4750	0.0507 ± 0.0803	0.0447 ± 0.0832
HR	77	0.3286 ± 0.0344	0.3468 ± 0.00006	3.439 ± 0.4117	0.0618 ± 0.0940	0.0525 ± 0.0991
Genetic group
G1	115	0.3416 ± 0.0192	0.3487 ± 0.00007	3.032 ± 0.5213	0.0239 ± 0.0527	0.0204 ± 0.0551
G2	18	0.3424 ± 0.0326	0.3331 ± 0.00003	3.667 ± 0.5966	0.1335 ± 0.0798	-0.0278 ± 0.0979
G3	20	0.3606 ± 0.0464	0.3397 ± 0.00007	3.435 ± 0.4571	0.1612 ± 0.1029	-0.0616 ± 0.1367
G4	12	0.3798 ± 0.0468	0.3495 ± 0.00002	3.728 ± 1.004	0.0850 ± 0.1105	-0.0866 ± 0.1339
Overall	170	0.3305 ± 0.0305	0.3470 ± 0.00006	3.406 ± 0.4072	0.0549 ± 0.0842	0.0476 ± 0.0881
Note: N: Samples size; H_O: observed heterozygosity; H_E: expected heterozygosity; ROH: Runs of homozygosity; F_ROH: ROH based inbreeding; F: inbreeding coefficient; LR: Low-resistant; HR: High-resistant

Table 2

Candidate regions, SNPs with topmost score in LR-GWAS, and associated genes found in candidate regions that were identified by at least two methods of selection signature analysis between the high- (HR) and low-resistance (LR) sheep cohorts.

Reg	OAR	Start	Stop	Size (Mb)	Method				Top SNP (LR-GWAS)	No. of genes	Genes (Top gene)**	P value*
Reg	OAR	Start	Stop	Size (Mb)	ROH	LR-GWAS	F_ST	XP-EHH	Top SNP (LR-GWAS)	No. of genes	Genes (Top gene)**	P value*
1	1	122490001	122680000	0.190	×	√	×	√	oar3_OAR1_122581225	4	*TIAM1*, LOC105616413, LOC105616409, LOC105616407	3.10E-06
2	1	226950001	227160000	0.210	×	√	√	×	oar3_OAR1_227042621	5	GFM1, LXN, MLF1, RSRC1, LOC105603699	2.55E-06
3	1	253380001	253660000	0.280	×	√	√	√	oar3_OAR1_253478525	4	LOC105604371, SLCO2A1, RAB6B, TF	2.70E-06
4	2	51910001	52140000	0.230	×	√	√	×	oar3_OAR2_52043939	7	RNF38, GNE, CLTA, CCIN, GLIPR2, LOC101111253, RECK	2.72E-06
5	2	53110001	53250000	0.140	×	×	√	√	oar3_OAR2_53171448	4	C2H9orf131, DNAJB5, KIAA1045, LOC101110212	1.79E-06
6	2	62970001	63160000	0.190	×	√	×	√	oar3_OAR2_63061644	3	ANXA1, LOC105607837, LOC105607840	2.68E-06
7	2	84550001	84740000	0.190	×	√	×	√	oar3_OAR2_84646297	4	*BNC2*, LOC105608204, LOC105608202, LOC105608200	3.47E-06
8	2	85810001	86140000	0.330	×	√	×	√	oar3_OAR2_86046574	2	*ADAMTSL1*, TRNAC-GCA	3.18E-06
9	2	201150001	201380000	0.230	×	√	√	√	oar3_OAR2_201281356	5	LOC101102297, C2H2orf69, TYW5, C2H2orf47, LOC105610167	2.72E-06
10	3	183020001	183260000	0.240	√	√	√	×	oar3_OAR3_183165082	2	TRNAC-GCA, DENND5B	2.76E-06
11	3	193980001	194170000	0.190	×	√	√	×	oar3_OAR3_194077322	1	PDE3A	2.65E-06
12	3	213290001	213550000	0.260	×	√	√	×	oar3_OAR3_213385724	11	PEX26, TUBA8, CDC42EP1, LGALS2, GGA1, SH3BP1, PDXP, LGALS1, NOL12, TRIOBP, TRNAG-CCC	2.70E-06
13	3	224020001	224277781	0.258	×	√	√	×	oar3_OAR3_224118862	11	NCAPH2, ODF3B, KLHDC7B, SYCE3, CPT1B, CHKB, MAPK8IP2, ARSA, SHANK3, ACR, RABL2B	3.35E-06
14	4	48760001	49000000	0.240	×	√	√	×	oar3_OAR4_48859309	5	CBLL1, SLC26A3, LOC105611149, DLD, LAMB1	2.82E-06
15	5	44110001	44300000	0.190	×	√	×	√	oar3_OAR5_44206370	3	LOC105607063, LOC101104574, LOC105607061	2.68E-06
16	5	105810001	106000000	0.190	×	√	√	×	oar3_OAR5_105908550	1	*MAN2A1*	2.56E-06
17	6	5580001	5770000	0.190	×	√	√	×	oar3_OAR6_5672674	2	LOC101115459, MAD2L1	2.60E-06
18	8	2520001	2710000	0.190	×	√	×	√	oar3_OAR8_2618920	2	*MYO*, IMPG1	2.55E-06
19	8	27780001	27970000	0.190	×	√	√	√	oar3_OAR8_27875909	7	AK9, ZBTB24, MICAL1, SMPD2, PPIL6, POLE4, CD164	2.72E-06
20	8	63260001	63450000	0.190	×	√	×	√	OAR8_68144751.1	1	*NHSL1*	2.59E-06
21	9	27550001	27740000	0.190	×	√	×	√	oar3_OAR9_27644689	3	*TRIB*, LOC105611791, NSMCE2	2.63E-06
22	9	33820001	34010000	0.190	×	√	√	×	oar3_OAR9_33912632	2	*PCMTD1*, ST18	2.60E-06
23	9	59370001	59560000	0.190	×	√	×	√	oar3_OAR9_59469003	1	*LOC105611701*	2.70E-06
24	11	6640001	6830000	0.190	×	√	×	√	oar3_OAR11_6733973	1	*ANKFN1*	2.82E-06
25	12	29800001	30030000	0.230	×	√	√	×	oar3_OAR12_29939034	2	*KIF26B*, LOC105616515	2.71E-06
26	12	54300001	54490000	0.190	×	√	√	×	oar3_OAR12_54393290	4	*TNN*, LOC105612227, KIAA0040, LOC105612225	2.58E-06
27	13	37980001	38190000	0.210	×	√	√	√	oar3_OAR13_38099599	2	*SLC24A3*, LOC105616699	2.30E-06
28	15	958579	1427669	0.469	√	×	×	√	s55467.1	5	MRE11A, ANKRD49, AASDHPPT, KBTBD3, MSANTD4	1.46E-06
29	15	5510001	5720000	0.210	×	√	√	×	oar3_OAR15_5627753	3	MMP27, MMP20, MMP7	2.73E-06
30	15	17290001	17690000	0.400	×	√	√	√	oar3_OAR15_17387492	3	NPAT, ATM, C15H11orf65	3.01E-06
31	15	20560001	20730000	0.170	×	√	√	√	oar3_OAR15_20686078	1	*LOC105602151*	1.71E-06
32	16	1420001	1610000	0.190	×	√	√	×	oar3_OAR16_1510505	4	*LOC105602434*, SPDL1, LOC105602438, DOCK2	2.68E-06
33	16	2280001	2470000	0.190	×	√	√	×	oar3_OAR16_2373107	3	KCNIP1, KCNMB1, LOC105602444	2.93E-06
34	16	12310001	12500000	0.190	×	√	×	√	oar3_OAR16_12408278	4	LOC105602490, CD180, LOC105602489, MAST4	2.55E-06
35	16	16060001	16250000	0.190	×	√	√	√	oar3_OAR16_16150781	1	*LOC105602500*, LOC105602501	3.40E-06
36	16	53590001	53780000	0.190	×	√	√	×	oar3_OAR16_53686329	2	*CDH18*, TRNAW-CCA	2.60E-06
37	18	23520001	23880000	0.360	√	√	√	×	oar3_OAR18_23760534	2	EFTUD1, MEX3B	2.65E-06
38	18	33990001	34180000	0.190	×	√	×	√	oar3_OAR18_34086919	1	*STXBP6*	2.59E-06
39	18	39710001	40000000	0.290	×	√	√	×	oar3_OAR18_39823086	1	*LOC105603171*	3.14E-06
40	19	880001	1070000	0.190	×	√	√	×	oar3_OAR19_977383	3	LOC105603323, LOC105603324, LOC105603325	2.60E-06
41	19	9910001	10100000	0.190	×	√	√	×	oar3_OAR19_10008059	1	*STAC*	2.64E-06
42	21	38480001	38690000	0.210	×	√	√	×	oar3_OAR21_38598176	7	LOC105601927, LOC101106045, LOC101104869, LOC101106295, LOC105604093, LOC101106556, LOC101106809	2.89E-06
43	21	38646061	38906932	0.261	√	√	√	×	oar3_OAR21_38573588	9	LOC101106556, LOC101106809, LOC101105115, LOC101107054, PAG6, LOC105604095, LOC443348, LOC101107317, LOC101107562	2.83E-06
44	21	40520001	40710000	0.190	×	√	√	×	oar3_OAR21_40617710	9	STX5, WDR74, LOC101116519, SLC3A2, CHRM1, LOC101112776, TRNAE-UUC, LOC101117286, SLC22A8	2.89E-06
45	25	19550001	19740000	0.190	×	√	√	×	oar3_OAR25_19649912	2	REEP3, LOC101122571	2.60E-06
										160
Significant markers following LR-GWAS (Bonferroni correction P < 0.001); * Top genes closest to the topmost SNP are indicated in bold.

No competing interests reported.

SupplementaryInformationTables.docx

Download PDF

Journal Publication

published 28 Feb, 2022

Read the published version in BMC Genomics →

Editorial decision: Major revision
13 Oct, 2021
Reviews received at journal
28 Sep, 2021
Reviewers agreed at journal
23 Sep, 2021
Reviewers invited by journal
23 Aug, 2021
Editor assigned by journal
23 Aug, 2021
Editor invited by journal
20 Aug, 2021
Submission checks completed at journal
20 Aug, 2021
First submitted to journal
07 Aug, 2021

You are reading this latest preprint version

Genome Variation in Tick Infestation and Cryptic Divergence in Tunisian Indigenous Sheep

Status:

Journal Publication

Version 1

Abstract

Figures

Background

Results

Population genetic structure

Genome variation

Genome-wide selection signature analysis

Discussion

Conclusions

Methods

Study animals, sampling and genotyping

Assessment of population genetic structure

Assessment of genetic diversity and demographic dynamics

Detection of selection signatures and association analysis

Functional annotation of candidate regions

Abbreviations

Declarations

Ethics approval and consent to participate

Consent for publication

Availability of data and materials

Competing interests

Funding

Author’s contributions

Acknowledgements

References

Tables

Additional Declarations

Supplementary Files

Status:

Journal Publication

Version 1