Skip to main content
Download PDF

Genome variation in tick infestation and cryptic divergence in Tunisian indigenous sheep

BMC Genomics volume 23, Article number: 167 (2022) Cite this article

Abstract

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 its genetic control remains unknown.

Results

Here, we report 397 genome-wide signatures of selection overlapping 991 genes from the analysis, using ROH, LR-GWAS, XP-EHH, and FST, of 600 K SNP genotype data from 165 Tunisian sheep showing high and low levels of tick infestations and piroplasm infections. We consider 45 signatures that are detected by consensus results of at least two methods as high-confidence selection 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 revealed fine-scale genome structure indicative of cryptic divergence in Tunisian sheep.

Conclusions

Our findings provide a genomic reference that can enhance the understanding of the genetic architecture of tick resistance and cryptic divergence in indigenous African sheep.

Peer Review reports

Background

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 [1,2,3,4,5]. These induce vasodilatory, anti-haemostatic and immunomodulatory activities that guarantees their successful engorgement and acquisition of a blood meal [1,2,3,4,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 and animals [6, 7]. Globally, ticks and tick-borne infections (T-TBI) are a major health impediment to livestock performance especially in the (sub-)tropics where the host, vector and pathogen overlap. So far, the impact of ticks on sheep production has not been investigated to the same extent as in cattle. The available estimates show that T-TBI affect about 80% of world’s cattle, with a global annual loss approximating US$7 billion [8, 9]. T-TBI negatively impact cattle production by reducing total milk yields by up to 23% [10] and the value of skins and hides by 20-30% [11, 12]. The negative effects of ticks on liveweights are also high. On average, an engorged female tick results in a body weight reduction of 1.37 ± 0.25 and 1.18 ± 0.21 g in Bos taurus and B. taurus x B. indicus crossbred cattle [13], respectively. In China estimates show that the total annual losses due to TBI in small ruminants approximates US$70 million [14].

Prophylactic use of acaricides is the most common strategy of controlling and eradicating ticks [15, 16]. The (over) use of acaricides has however imposed selection pressure resulting in the development of acaricide-resistant tick strains, environmental contamination, retention of chemical residues in livestock products and increased costs of existing, and of developing and manufacturing new and more potent, acaricides [16,17,18,19,20]. Anti-tick vaccines are a promising option, but they do not confer protection against multiple tick species [7, 21, 22] which is a common occurrence especially in the (sub-)tropics [23,24,25]. These factors are driving the search for alternative control strategies such as the use of tick-resistant animal genotypes as the host’s natural resistance to ticks can be exploited as a long-term sustainable control option targeting most 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 [26]). Most of the studies have shown that an effective immunological response (resistance/tolerance) to T-TBI, is genetically determined with a heritability range of 0.05-0.42 (see [26]). For haemopathogens, the immunological response has been associated with the ability to resist the development of anemia following infection [25, 26], variation in the expression of extracellular matrix metalloproteinase [27], differences in chemokines and their receptors and toll-like receptors [28, 29], variation in genes that limit the supply of blood-meal to ticks and genes that enhance innate and adaptive immune responses [30]. Several genes have been proposed as candidates for tick resistance. Lumican (LUM) was identified as a potential biomarker for tick resistance in cattle [30, 31] reported two SNPs in ELTD1 associated with tick burdens in both dairy and beef cattle, and a haplotype of nine tag SNPs and two others associated with tick counts in dairy cattle. The authors also reported haplotypes spanning ITGA11 associated with tick burdens. Class I antigens of the bovine major histocompatibility complex have also been associated with tick loads [32, 33] and alleles at the bovine lymphocyte antigen (BOLA-DRB3) have been linked with resistance to Rhipicephalus (Boophilus) microplus infestation [34,35,36].

Domestic sheep (Ovis aries) are the second most abundant ruminant livestock after cattle [37] 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.) [38,39,40,41,42,43,44]. Some of these pathogens cause important zoonotic diseases such as Crimean Congo Haemorrhagic fever, Q fever and Human granulocytic anaplasmosis [45] resulting in negative impacts on human and animal health, and significant socio-cultural and economic losses. Despite reports on individual variations in tick burdens in sheep including prevalence and infestation intensity [46,47,48,49], little is known regarding the genetic basis of this phenotype. In this study, taking advantage of observed individual natural variation in tick infestation in two indigenous sheep (Barbarine (B; fat-tailed) and Queue Fine de ľOuest (Q; thin-tailed)) from Tunisia [44], we investigated the possible genetic basis of the trait through the analysis of genome-wide 600 K SNP genotype data from 165 individuals showing high- and low-infestation (HR/LR) to ticks and piroplasm infections. Our results suggest tick resistance could be the outcome of multigene associations with STX5 being the possible prime, and RAB6B, TF, SLCO2A1 and STXBP6 being the likely potential, candidate genes driving genetic variation for tick infestation in sheep.

Results

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 (Fig. 1). They separated the study individuals into four genetic groups, herein named G1, G2, G3 and G4 (Fig. 1a). These four groups did not correspond to the resistance status to ticks (HR/LR) (Fig. 1b), the sampling region (northeast, northwest, southeast) (Fig. 1c) and breed (B/Q) (Fig. 1d).

Fig. 1
figure 1

Genetic structure of the two cohorts (HR (high resistant), LR (low resistant)) of Tunisian sheep. a, b, c and d PCA cluster analysis showing PC1 and PC2; e Cross-validation plot for admixture analysis;f Admixture analysis plot showing the genetic backgrounds present in the study cohorts for 2≤K≤5

Full size image

Following ADMIXTURE analysis, the lowest CV error was at K = 4 (Fig. 1e) suggesting at least four genetic clusters in the dataset. The clusters corresponded to the four genetic groups revealed by PCA but not to the resistance status to ticks, sampling region and breed. We therefore named the four genetic clusters as G1, G2, G3 and G4 (Fig. 1f). The G1 cluster occurs in a few individuals of G2, G3 and G4, pointing to its possible introgression into the latter.

Genome variation

Descriptive statistics, providing insights on genetic diversity were estimated for the overall dataset, the two breeds (B/Q), the two cohorts (HR/LR) and the four genetic groups (G1, G2, G3, G4) (Table 1). The overall estimates of genetic diversity represented by observed (HO) and expected (HE) heterozygosity averaged 0.3305 ± 0.0305 and 0.3470 ± 0.00006, respectively. The individual values for each statistic were in general higher than 0.3286 ± 0.0344. The mean length of ROH was 3.406 ± 0.4072 and ranged between 3.032 ± 0.5213 Mb (G1) and 3.728 ± 1.004 Mb (G4). The mean value of FROH 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 between -0.0866 ± 0.1339 (G4) and 0.0525 ± 0.0991 (HR). The LR cohort showed higher values of HO and HE compared to HR. At the opposite, the LR cohort showed lower values of ROH size, FROH and F. For the four genetic groups, G4 had the highest values of HO, HE and mean length of ROH while the lowest values were observed in G1, G2 and G1, respectively. The FROH was highest in G3 and lowest in G1. Except G1, the other three genetic groups showed negative values of F.

Table 1 Indices of genetic diversity generated for the two breeds, two cohorts and four genetic groups of Tunisian Sheep analysed in this study
Full size table

We assessed the decay in linkage disequilibrium (LD; Fig. 2a) and the trends in effective population size (Ne; Fig. 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 was the same across all the classes of datasets analysed. It generally revealed higher LD at shorter distances which decayed rapidly plateauing off around 0.2 Mb. The G2, G3 and G4 showed persistently higher r2 values (r2 > 0.15) and thus higher LD compared to the overall dataset, HR, LR and G1 (r2 < 0.15). The trends in Ne over generation time showed different profiles for the data classes analysed. The overall dataset and G1 group showed the highest Ne which increased gradually reaching maxima around 350 generations ago and then declined rapidly to the present time. The HR and LR cohorts showed the next lowest estimate of Ne which increased gradually up to around 350 generations ago and was followed by a rapid decline to the present. The G2, G3 and G4, had the lowest Ne across all generations which declined gradually to the present time.

Fig. 2
figure 2

Trends in LD decay (a) and Ne(b) across 1000 generations for the two cohorts (HR (high resistant), LR (low resistant)) and four genetic backgrounds (G1, G2, G3 and G4) of Tunisian sheep

Full size image

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 (FST) to gain insights on the regions of the genome correlating with resistance to ticks (HR verses LR) considering equal sample size of 74 animals in both cohorts.

The ROH approach identified 110 and 105 ROH regions spanning 280 and 281 annotated genes in the HR and LR cohorts (Fig. 3a and b; Supplementary Table S1, S2), respectively. The LR-GWAS identified 242 candidate regions spanning 561 genes (Fig. 4a), and FST identified 79 candidate regions, showing genetic differentiation between HR and LR (Fig. 4b), and spanning 243 genes. The XP-EHH identified 76 candidate regions spanning 187 genes (Fig. 4c). These 397 candidate regions overlapped 991 genes (Supplementary Table S3, S4, S5).

Fig. 3
figure 3

Manhattan plots showing the genome-wide distribution frequency of SNPs in stretches of ROH regions. The dashed lines indicate the 50% threshold for each cohort (HR (high resistant), LR (low resistant)) of Tunisian sheep investigated here

Full size image
Fig. 4
figure 4

Manhattan plots showing the genome-wide distribution of SNPs following a LR-GWAS, b FST and c XP-EHH analysis using the two cohorts (HR (high resistant), LR (low resistant)) of Tunisian sheep analysed in this study

Full size image

We used the ROH approach to identify cohort-specific selection signatures. We therefore compared the ROH results of HR and LR (Fig. 4a and b) and identified 30 ROH regions, spanning 57 annotated and 20 uncharacterized (prefixed “LOC”) genes that were specific to the HR cohort (Supplementary Table S6). Considering these HR-specific ROH regions and those identified by LR-GWAS, XP-EHH and FST, 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 uncharacterized (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 as identified by LR-GWAS in the respective candidate regions. Four genes (RAB6B (OAR1), TF (OAR1), SLCO2A1 (OAR1), STXBP6 (OAR18) and STX5 (OAR21)) found in our candidate regions are of particular interest with regard to the objectives of this study. They have been observed to be directly or indirectly associated with body weight traits in various livestock species and resistance to ticks in cattle.

Table 2 Candidate regions, SNPs with topmost score in LR-GWAS, and associated genes that were identified by at least two methods of selection signature analysis between the high- (HR) and low-resistant (LR) sheep cohorts
Full size table

Functional enrichment analysis was performed in the pool of 104 putative genes found in the 45 candidate regions that overlapped between at least one HR-specific ROH region and the ones identified by LR-GWAS, FST and XP-EHH, after excluding 56 genes that remain uncharacterized (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 terms were represented by clusters of genes (such as GGA1, CD164, and BCL2L11) having roles in innate immunity and disease-related inflammation [50, 51].

Discussion

T-TBI result in negative impacts on ruminant livestock production. Therefore, the control and ultimate elimination of T-TBI should be prioritized to minimize impacts not only on animal health and production but also on human and environmental wellbeing. Large variations in susceptibility to T-TBI points to genome-wide variability that underpins inter-animal variation in tick infestation and the pathogens they transmit. Here, we generated ovine 600 K SNP genotype data, from 165 animals of two breeds of Tunisian indigenous sheep that graze natural communal pastures and with no history of anti-tick prevention intervention. The sheep comprised of individuals showing high and low resistance (HR/LR) to tick infestation and piroplasm infection [44]. We analysed the data with ROH, LR-GWAS, XP-EHH and FST, to investigate signatures of selection associated with variation in tick infestation and thus possible resistance. As noted by [52], we also acknowledge that using naturally infected animals to study the genetic basis of resistance runs the risk of resistant animals being a cocktail of truly highly resistant individuals as well as those that were never exposed to infestations/infections. These factors may dilute the certainty that animals showing high or low resistance to different pathogens and parasites, could in fact be functionally dissimilar [52].

Compared to the HR, the LR cohort showed comparatively higher levels of genetic diversity but lower levels of inbreeding, though the differences were not significant (P > 0.05). 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 variation have also been reported in Egyptian [53], Algerian [54], Tunisian [55] and Russian [56] breeds of sheep.

The PCA and ADMIXTURE provided corroborating evidence suggesting the absence of genetic stratification that is consistent with the a priori classification of the study individuals based on their susceptibility/resistance to tick infestation/piroplasm infection, geographic sampling regions and breeds. The lack of genetic differentiation corresponding to susceptibility/resistance to tick infestation/piroplasm infection is not surprising. A similar lack of genetic stratification corresponding to a priori classification of sheep based on levels of prolificacy [57] or resistance to gastro-intestinal nematodes [52, 58] has been reported. [59] also observed no clear genetic differentiation between sub-populations of dual-purpose Black and White and German Holstein cattle. Our result suggests that the animals comprising the HR and LR cohorts are not highly divergent for the tick-resistance phenotype. This can be attributed to one or more of the following factors, (i) weak selection pressure driving differences in tick susceptibility/resistance and therefore any favourable alleles are likely to be rare, (ii) low level of tick burden that does not result in detectable genomic differentiation, and (iii) absence of farmer-driven preferential use of tick-resistant animals for breeding. 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 different sheep breeds in close geographic proximity and contact [60] resulting in cross-mating that is homogenizing their genomes. A similar finding was reported by [55]. The fact that B and Q could not be differentiated by ADMIXTURE provides evidence at the genome level that supports past and on-going mating of the fat-tailed Barbarine with thin-tailed sheep to reduce the tail-fat content in the Barbarine carcass [61].

LD was estimated for all marker pairs using the r2 metric plotted as a function of increasing genomic distance. The overall pattern of LD decay 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 100 Kb. This differs from what has been observed in commercial breeds whose LD plateaus off at around 150 Kb [62] and is most likely the result of differences in breeding history. Specifically, the study populations comprised a broad genetic base while that of commercial breeds has been narrowed down by bottlenecks arising from breed formation and stringent artificial selection for economic traits. Both the HR and LR cohorts showed identical patterns of LD decay, and trends in Ne. Although the two cohorts showed differences in tick infestation [44], this result suggests that they share aspects of their past and recent population demography and breeding history. This may be the case as their genomes showed no differentiation aligning with their classification based on tick susceptibility/resistance. We also assessed the demographic profiles of the four genetic groups underlying the genome architecture of the study individuals as revealed by PCA and ADMIXTURE. The G2, G3 and G4 genetic groups had the highest r2 values (r2 > 0.15) and thus highest LD over short genomic distances and a decay curve with persistently elevated values compared to G1. A similar pattern was observed in the primitive Soay sheep [62] 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 Ne which declines gradually over all generations. The possible reason(s) for the gradual decline in Ne in the three groups remains unknown and is worth investigating. In contrast, LD was lower (r2 < 0.15) in G1 across all genomic distance intervals and the group also showed the highest Ne 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 [62]. This explanation may however not explain our results as there was no significant differences in the level of genetic diversity across the different classes of datasets analysed herein. Despite being unexplained, we postulate that the low LD and high Ne in G1, and similar magnitude of genetic diversity is unlikely to reflect any biological phenomenon.

Since millennia, ticks have parasitized animals for blood-meal to the extent of developing, through co-evolution, a sophisticated armoury that guarantees their biological success. To explore genomic signatures encoding individual variation in tick infestation and thus host resistance to ticks, we segregated Tunisian sheep into two extreme groups comprising animals showing high- and low-resistance (HR/LR) to tick infestations and piroplasm infections following the results of a previous study on the same animals [42]. 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 [44] and their comparison was thus 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 FST approaches to detect selection signatures and SNPs that are likely associated with variation and thus resistance to tick infestation.

The four approaches revealed 45 candidate regions, that overlapped between at least two approaches, and which spanned 104 characterised genes. The LR-GWAS and FST identified a selection signature on OAR21 spanning nine genes, one of which was STX5 (Syntaxin-5). In a study of Belmont Red cattle, STX5 was amongst 11 of 14 genes that showed a significant increase in its level of expression in the skin of animals that were highly resistant to ticks at time zero post infestation compared to animals of low resistance [63]. The expression level was more pronounced in 3-hour skin samples, suggesting a response to tick attachment which could contribute to host innate immunity and 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 (Ca2+) release channel [64]. Based on the increased expression levels of the 11 (most of which are Ca2+ dependent genes) out of 14 genes that they tested, [63] suggested that the high mRNA transcription levels of Ca2+ signaling genes in skin of HR animals may explain their increased resistance to ticks. Previous tick exposure may prime animals that exhibit the HR phenotype to resist further infestations via increased expression of Ca2+ signaling proteins and based on these results, we suggest that STX5 could be a prime candidate gene driving tick resistance in sheep.

Environmental changes can exert positive or negative effects on mechanisms of thermoregulation that can influence tick burdens [65]. It has been observed that tick burdens in cattle might be correlated with traits that influence thermal comfort [66]. For instance, traits such as skin thickness, hair density and skin secretions can influence tick resistance and thermal comfort in domestic livestock as they affect the ability of the animal to dissipate heat [67]. 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 [66]. These findings are relevant to our study because one of our candidate regions that was revealed by LR-GWAS and XP-EHH on OAR18 spanned STXBP6 (Syntaxin Binding Protein 6) and another one that was revealed by LR-GWAS, XP-EHH and FST on OAR1 spanned SLCO2A1 (OATP2A1). Both genes were closest to the most significant SNP as revealed by LR-GWAS. STXBP6 was one among a cluster of genes found to be upregulated in the testes of roosters exposed to acute heat stress [68]. Transcripts of STXBP6 were also found to be among eight of other genes that were correlated with the modified Rodnan skin thickness score and forced vital capacity in humans suffering from scleroderma and systemic sclerosis [69], a condition that is characterized by thickening and hardening of the skin. Skin thickness is a physical barrier that can confer resistance to tick infestation(s); animals with thin skins having significantly lower tick burdens [63, 70, 71]. With the exception of birds and canids, thin skins allow animals to dissipate more heat through sweating and evaporative cooling when ambient temperatures are above thermoneutrality [65], 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 [72]. The function of these genes is most likely complemented by 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 [73, 74]. Taken together, this information led 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 T-TBI as well as thermoregulation. However, these hypotheses need to be validated with more data on appropriate phenotypes such as skin and coat characteristics.

Birth weight is the earliest available body weight trait with considerable impact on lamb survival and growth performance [75]. Our analysis revealed a candidate region on OAR1 that was identified by LR-GWAS, XP-EHH and FST that spanned among others the RAB6B and TF genes. This region was observed, following GWAS analysis, to be significantly associated with birth weight in sheep [76, 77]. 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 [78] and was found to be in a region under selection in broilers [79] and layers [80] suggesting that it may influence body weights. In several cases, T-TBI can destabilize host growth and development. To counter against such destabilization, we suggest that natural selection may be acting on the regions spanning RAB6B, TF and STXBP6 to enhance growth and development stability early in life as an adaptive strategy for survival in the context of high T-TBI challenge.

Conclusions

In this study, we investigated selection signatures in Tunisian indigenous sheep using 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. 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 driving tick resistance in sheep. We further hypothesized that STXBP6 and SLCO2A1 could be potential candidates for tick resistance in indigenous sheep and 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 may be enhancing growth and development stability as an adaptive strategy to tick infestation. For these genes to qualify as candidates for enhancing genetic resistance to ticks in sheep through precision breeding (genomic selection and/or gene editing), their potential effects should be quantified through gene expression studies involving resistance and susceptible animals and identify actionable variants encoding the trait. Such quantification would benefit from the inclusion of variables, such as skin and coat characteristics, to investigate their influence on tick infestation and resistance. Overall, our findings provide a genomic reference for understanding the genetic architecture of tick resistance and cryptic divergence in indigenous sheep.

Methods

Study animals, sampling and genotyping

The blood samples used in this study were collected by veterinary surgeons from farmers flocks. The sampling was done during routine animal health monitoring and surveillance following standard veterinary procedures that complied with animal welfare regulations as detailed in the guidelines for the care and use of animals of the Tunisian National Council of Veterinary Surgeons (TNCVS). Prior to blood sampling, verbal permission was sought and granted by animal owners who witnessed the procedure. Therefore, permission from the Ethics Committee of the TNCVS was not required.

Two breeds, Barbarine (B) and Queue Fine de l’Ouest (Q) from Tunisia were sampled. The animals graze natural communal pastures throughout the year except during summer when they are released to forage 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 selected for this study were sampled as part of a larger cross-sectional study on T-TBI carried out between 2018 and 2020 [44] in the northeast, northwest and southeast regions of Tunisia. The northeast and southeast regions are the homelands of the B breed although it has been dispersed throughout the country. It is the foundation of the “Tunis” and “Barbaresca” breeds found respectively, in the USA and Italy [81]. The B is also found in Libya and eastern Algeria [82]. The Q breed predominates in the northwest region although they have also 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 [83].

For the cross-sectional study, 439 mature ewes were tagged and monitored through eight rounds of sampling as detailed by [44]. 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, 165 (B = 105; Q = 60) 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). We used the infestation score of [84] to determine tick infestation loads as follows: 0 (no ticks), 1 (≤10), 2 (11–30), 3 (31–80), 4 (81–150), and 5 (>150), while piroplasms were detected with polymerase chain reaction [44].

From the results of population structure analysis, which showed no genetic differentiation between the two breeds, the 165 samples were assigned to two cohorts, irrespective of the breed based on extreme values of tick load and piroplasm infections viz.: (i) high tick-resistant (HR) - animals with no tick infestation(s) (score = 0) and/or no piroplasm infection(s) and (ii) low tick-resistant (LR) - animals that were highly infested with ticks (score > 81) and/or infected by piroplasms (Table 1). The DNA from the 165 samples was genotyped with the Illumina 600 K 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 [85].

The raw genotypes were assessed for quality with PLINK1.9 [86] 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 representing window size, step size and r2 threshold, respectively resulting in 338,180 SNPs that were used for genetic structure 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 [87], independent of 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 five replicate runs were performed 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 (HE) and observed (HO) heterozygosity, patterns of LD decay across the genome and effective population size (Ne) across generation time were investigated for each breed, the HR and LR cohorts and the genetic groups revealed by PCA and ADMIXTURE. HE and HO were calculated with PLINK v1.9. The extent of LD between pairwise SNPs was evaluated with the r2 statistic calculated with PLINK v1.9. The Ne over generation time was estimated with the equation Ne t = (1/4c) (1/r2 − 1) [88], where Ne t is the Ne t generations ago (t = 1/2c); r2 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 “FROH”. The latter was computed as the ratio of the total length of ROH to the length of autosomes (2.45 Gb) [89]. 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 from 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 after standardizing the sample sizes to 74 animals per 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 an ROH under selection. To identify the ROH streatches that are associated with tick resistance, we contrasted the HR and LR ROH regions and identified the ones that were specific to HR.

We used the population differentiation statistic, FST, to investigate regions of the genome showing divergence between HR and LR. The unbiased pairwise FST [90] was computed for each SNP with the R package “HIERFSTAT” [91] 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 FST values were standardized by Z transformation following [92]. To minimize the likelihood of false positive results, the windows that occurred within the top 0.001% of the normal distribution of the FST values in each chromosome and representing the most divergent regions between the two cohorts were considered putative candidates under divergent selection.

We used the software developed by [93] 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 [94,95,96]. We determined the candidate regions under positive selection by clustering the significant core SNPs (P-value < 0.05) within a distance of less than 100 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 the likelihood of false positive results. 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 top-most SNP found within the candidate regions identified. The Manhattan plots showing the genome-wide distribution of the SNPs were generated with the R package “qqman” v3.5.1.

Functional annotation of candidate regions

The candidate regions revealed by ROH were analysed and the ones that were specific to HR identified. These HR-specific ROH regions and the ones identified by LR-GWAS, XP-EHH and FST were analysed, and the ones that overlapped between at least two methods identified and merged using Bedtools v.2.28.0 [97]. 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. 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 [98]. We further determined the functions of the putative gene from the NCBI database (http://www.ncbi.nlm.nih.gov/gene/) and review of literature.

Availability of data and materials

The datasets generated and/or analysed during the current study are available in FIGSHARE repository, https://figshare.com/s/56deab19b63db8159c88; with the doi: https://doi.org/10.6084/m9.figshare.16915600.

Abbreviations

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

FST :

genetic differentiation statistic

OAR:

Ovine chromosome

References

  1. Wikel SK. Host immunity to ticks. Annu Rev Entomol. 1996;41:1–22.

    CAS  PubMed  Google Scholar 

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

    Google Scholar 

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

    CAS  PubMed  Google Scholar 

  4. Francischetti IMB, Sa-Nunes A, Mans BJ, Santos IM, Ribeiro JMC. The role of saliva in tick feeding. Front Biosci. 2009;14:2051–88. https://doi.org/10.2741/3363.

    Article  CAS  PubMed Central  Google Scholar 

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

    CAS  PubMed  Google Scholar 

  6. Jongejan F, Uilenberg G. The global importance of ticks. Parasitology. 2004;129 SUPPL.

    Google Scholar 

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

    PubMed  Google Scholar 

  8. Burrow HM, Mans BJ, Cardoso FF, Birkett MA, Kotze AC, Hayes BJ, et al. Towards a new phenotype for tick resistance in beef and dairy cattle: A review. Animal Production Science. 2019;59:1401–1427.

    CAS  Google Scholar 

  9. Ahmed J, Alp H, Aksin M, Seitzer U. Current status of ticks in Asia. In: Parasit. Rec. 2007. p. 159–62.

  10. Regitano L, Prayaga K. Ticks and tick-borne diseases in cattle. In: Bishop S, Axford R, Nicholas F, Owen J, editors. Breeding for Disease Resistance in Farm Animals. 3rd edition. London, UK: CAB International; 2011. p. 295–314.

  11. de La Fuente J, Kocan KM. Strategies for development of vaccines for control of ixodid tick species. Parasite Immunol. 2006;28:275–83.

    PubMed  Google Scholar 

  12. Ghosh S, Azhahianambi P, Borne MY-J of vector, 2007 U. Upcoming and future strategies of tick control: a review. J Vect Borne Dis. 2007;44:79–89. http://mrcindia.org/journal/issues/442079.pdf. Accessed 18 Oct 2021.

  13. Jonsson NN. The productivity effects of cattle tick (Boophilus microplus) infestation on cattle, with particular reference to Bos indicus cattle and their crosses. Vet. Parasitol. 2006;137:1–10.

    CAS  PubMed  Google Scholar 

  14. Yin H, Luo J. Ticks of small ruminants in China. In: Parasitology Research. 2007. p. 187–9.

  15. Mendes EC, Mendes MC, Sato ME. Diagnosis of amitraz resistance in Brazilian populations of Rhipicephalus (Boophilus) microplus (Acari: Ixodidae) with larval immersion test. Exp Appl Acarol. 2013;61:357–69.

  16. Abbas RZ, Zaman MA, Colwell DD, Gilleard J, Iqbal Z. Acaricide resistance in cattle ticks and approaches to its management: The state of play. Vet Parasitol. 2014;203:6-20.

  17. Willadsen P. Tick control: Thoughts on a research agenda. Vet Parasitol. 2006;138:161–8.

    PubMed  Google Scholar 

  18. Regitano LCA, Ibelli AMG, Gasparin G, Miyata M, Azevedo ALS, Coutinho LL, et al. On the Search for markers of tick resistance in bovines. Dev Biol. 2008;132:225–30.

  19. 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:280.

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

    PubMed  Google Scholar 

  21. de La Fuente J, Contreras M. Tick vaccines: Current status and future directions. Expert. Rev. Vaccines. 2015;14:1367–76.

    PubMed  Google Scholar 

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

  23. Kaiser MN, Sutherst RW, Bourne AS. Relationship between ticks and zebu cattle in southern Uganda. Trop Anim Health Prod. 1982;14.

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

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

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

    PubMed  Google Scholar 

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

    CAS  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

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

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

    PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    PubMed  Google Scholar 

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

    PubMed  Google Scholar 

  40. Gharbi M. Epidemiological Study of Sheep Anaplasmosis (Anaplasma ovis Infection) in Kairouan, Central Tunisia. Adv Parasitol. 2015;2:30–4.

    Google Scholar 

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

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

    CAS  PubMed  Google Scholar 

  43. Rjeibi MR, Darghouth MA, Gharbi M. Prevalence of Theileria and Babesia species in Tunisian sheep. Onderstepoort J Vet Res. 2016;83:a1040.

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

    Google Scholar 

  45. Alessandra T, Santo C. Tick-borne diseases in sheep and goats: Clinical and diagnostic aspects. Small Rumin Res. 2012;106 SUPPL.

    Google Scholar 

  46. Ogore PB, Baker RL, Kenyanjui M, Thorpe W. Assessment of natural ixodid tick infestations in sheep. Small Rumin Res. 1999;33:103–7.

    Google Scholar 

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

    CAS  PubMed  Google Scholar 

  48. 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:8.

  49. Grøva L. Tick-borne fever in sheep – production loss and preventive measures. Dept of Aquaculture and Animal Science. 2011;Ph.D.

  50. Moumène A, Meyer DF. Ehrlichia’s molecular tricks to manipulate their host cells. Microbes and Infection. 2016;18:172–9.

    PubMed  Google Scholar 

  51. Chen J, Ao L, Yang J. Long non-coding RNAs in diseases related to inflammation and immunity. Ann Transl Med. 2019;7:494–494.

    CAS  PubMed  PubMed Central  Google Scholar 

  52. 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:9250.

    CAS  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

  55. 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:106192.

  56. 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:29.

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

    CAS  PubMed  PubMed Central  Google Scholar 

  58. 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:73.

  59. 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:391-408.

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

    Google Scholar 

  61. Bedhiaf-Romdhani S, Djemali M, Zaklouta M, Iniguez L. Monitoring crossbreeding trends in native Tunisian sheep breeds. Small Rumin Res. 2008;74:274–8.

    Google Scholar 

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

    CAS  PubMed  Google Scholar 

  63. 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:506.

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  Google Scholar 

  66. 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:73.

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

    CAS  PubMed  Google Scholar 

  68. 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:374-82.e7.

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    PubMed  Google Scholar 

  71. Jonsson NN, Piper EK, Constantinoiu CC. Host resistance in cattle to infestation with the cattle tick Rhipicephalus microplus. Parasite Immunol. 2014;36:553–9.

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

  74. 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:84.

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

    Google Scholar 

  76. 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:e66569.

  77. Ghasemi M, Zamani P, Vatankhah M, Abdoli R. Genome-wide association study of birth weight in sheep. Animal. 2019;13:1797–803.

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  Google Scholar 

  79. Fu W, Lee WR, Abasht B. Detection of genomic signatures of recent selection in commercial broiler chickens. BMC Genetics. 2016;17.

  80. 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:e49525.

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

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

    PubMed  Google Scholar 

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

    Google Scholar 

  84. 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:413-25.

  85. Anderson R, McEwan J, Brauning R, Kijas J, Dalrymple J, Worley K, et al. Development of a high density (600 K) 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: ISGC; 2014.

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

    CAS  PubMed  PubMed Central  Google Scholar 

  87. Alexander DH, Novembre J, Lange K. Fast model-based estimation of ancestry in unrelated individuals. Genome Res. 2009;19:1655–64.

    CAS  PubMed  PubMed Central  Google Scholar 

  88. Sved JA. Linkage disequilibrium and homozygosity of chromosome segments in finite populations. Theor Popul Biol. 1971;2:125–41.

    CAS  PubMed  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

  90. Weir BS, Cockerham CC. Estimating F-statistics for the analysis of population structure. Evolution. 1984;38:1358–70.

    CAS  PubMed  Google Scholar 

  91. Goudet J. HIERFSTAT, a package for R to compute and test hierarchical F-statistics. Mol Ecol Notes. 2005;5:184–6.

    Google Scholar 

  92. 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:699.

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

    CAS  PubMed  PubMed Central  Google Scholar 

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

    CAS  PubMed  PubMed Central  Google Scholar 

  95. Ma. Y, Zhang. H, Zhang Q, Ding X. Identification of selection footprints on the X Chromosome in pig. PLoS One. 2014;9:e94911.

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

  97. Quinlan AR, Hall IM. BEDTools: A flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26:841–2.

    CAS  PubMed  PubMed Central  Google Scholar 

  98. Huang DW, Sherman BT, Lempicki RA. Systematic and integrative analysis of large gene lists using DAVID bioinformatics resources. Nat Protoc. 2009;4:44–57.

    CAS  Google Scholar 

Download references

Acknowledgements

We acknowledge the valuable contribution of veterinarians that facilitated access to sheep and to sheep owners for volunteering their flocks for 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 Dry 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 information

Authors and Affiliations

  1. Small Ruminant Genomics, International Centre for Agricultural Research in the Dry Areas (ICARDA), Addis Ababa, Ethiopia

    Abulgasim M. Ahbara, Aynalem Haile, Mourad Rekik, Barbara Rischkowsky & Joram M. Mwacharo

  2. Department of Zoology, Faculty of Sciences, Misurata University, Misurata, Libya

    Abulgasim M. Ahbara

  3. Lab. of Infectious Animal Diseases, Zoonosis and Sanitary Regulation, Univ. Manouba, Institution of Agricultural Research and Higher Education, National School of Veterinary Medicine of Sidi Thabet, 2020, Sidi Thabet, Tunisia

    Médiha Khamassi Khbou

  4. Laboratoire de Parasitologie, Institution de La Recherche et de L’Enseignement Supérieur Agricoles, Univ. Manouba, Ecole Nationale de Médecine Vétérinaire de Sidi Thabet, 2020, Sidi Thabet, Tunisia

    Médiha Khamassi Khbou, Rihab Rhomdhane, Limam Sassi & Mohamed Gharbi

  5. Animal and Veterinary Sciences Group, SRUC and Centre for Tropical Livestock Genetics and Health (CTLGH), Roslin Institute, Midlothian, Scotland

    Joram M. Mwacharo

Authors
  1. Abulgasim M. Ahbara
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Médiha Khamassi Khbou
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Rihab Rhomdhane
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Limam Sassi
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Mohamed Gharbi
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Aynalem Haile
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Mourad Rekik
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Barbara Rischkowsky
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. Joram M. Mwacharo
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

MG, MR, AH, BR and JMM conceived the study. MKK, RR and LS carried out the field work including sampling, DNA extraction and PCR under supervision of MG and MR. AMA analysed the data under supervision from JMM and both wrote the manuscript. All authors read and approved the final version of the manuscript.

Corresponding author

Correspondence to Joram M. Mwacharo.

Ethics declarations

Ethics approval and consent to participate

This study did not require any ethical approvals from the ethical committee of the National School of Veterinary Medicine for Animal Experimentation (NSVMAE) because the blood samples used for the study were collected by trained veterinary surgeons. Blood samples and ticks were collected using standard veterinary procedures in full respect of animal welfare and minimizing stress to the animals following the guidelines for the care and use of animals of the Tunisian National Council of Veterinary Surgeons. The blood was also sampled from the animals in the presence of their owners who provided verbal informed consent for their sampling according to the Charter of the Ethical Committee of the NSVMAE.

Consent for publication

Not applicable.

Competing interests

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

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Ahbara, A.M., Khbou, M.K., Rhomdhane, R. et al. Genome variation in tick infestation and cryptic divergence in Tunisian indigenous sheep. BMC Genomics 23, 167 (2022). https://doi.org/10.1186/s12864-022-08321-1

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-022-08321-1

BMC Genomics

ISSN: 1471-2164

Contact us