Discovery of host genetic factors through multi-locus GWAS against toxoplasmosis in sheep: addressing one health perspectives

Abstract

Toxoplasma gondii stands as one of the most successful pathogens, capable of infecting nearly all warm-blooded species. It is estimated that up to 50% of human population might harbor Toxoplasmosis infections. One of the primary transmission routes is the consumption of tissue cysts from infected farm animals used for food production. Thus, controlling Toxoplasmosis in farm animals is of vital importance for human health and food safety. Selective breeding in farm animals, where available, could complement classical control measures like biosecurity measures, vaccination, and test-and-cull methods. This multidisciplinary approach will make the eradication of Toxoplasmosis more effective. For this purpose, we conducted four multi-locus genome-wide association (GWA) approaches to identify the polygenic factors underlying innate resistance to Toxoplasma gondii in naturally infected sheep. Our findings indicate that 16 single nucleotide polymorphisms (SNPs), exhibiting varying degrees of statistical power, play a significant role in host immunity against T. gondii infection. We propose the genes containing these SNPs or located within 100 ± Kb of them (PLSCR5, EPHA3, DGKB, IL12B, CGA, WDR64, TMEM158, CLMP, and SIAE) as potential candidate genes. This study represents the first exploration of host genetic factors against Toxoplasmosis in livestock, utilizing the ovine paradigm as its foundation.

Introduction

The One Health approach acknowledges the intricate link between human, animal, and environmental health. It seeks optimal health outcomes through interdisciplinary collaboration. Approximately 60% of known and 75% of newly emerging infectious diseases in humans originate from animals [1]. Toxoplasmosis is a major One Health concern because it affects human, domestic animal, and wildlife health. It also disrupts ecosystems and poses risks to people who depend on animal resources [2, 3].

The etiological agent, Toxoplasma gondii, is an apicomplexan parasite with a complex life cycle. It replicates sexually in felids and asexually in various intermediate hosts, infecting nearly all warm-blooded species [4, 5]. T. gondii induces epigenetic modifications in rodent amygdalae, reducing their innate fear of felids [6, 7]. Transmission occurs via three infective stages: oocysts, tachyzoites, and bradyzoites. Felids shed oocysts in feces, which become infectious after sporulation. Bradyzoites, the dormant form, persist in tissue cysts and spread through ingestion. Both oocysts and bradyzoites transform into tachyzoites once inside a host, rapidly invading and multiplying within cells [8]. Transmission risk factors include consumption of undercooked meat, unpasteurized milk, and contaminated soil, with multiple cat ownership possibly increasing infection likelihood, although this remains debated [9, 10]. Congenital transmission occurs most frequently in late gestation, with early pregnancy infections leading to severe fetal complications [11, 12]. In healthy individuals, Toxoplasmosis is usually asymptomatic but can be life-threatening in immunocompromised individuals [13]. In the U.S., it is the second leading foodborne zoonotic cause of death [14]. Some studies associate T. gondii with behavioral disorders, although a causal relationship is yet to be established [15,16,17]. Global prevalence ranges from 10 to 80%, with lower rates in Northern Europe and North America and higher rates in tropical regions [5, 12, 18].

Toxoplasmosis is widespread in farm animals, with infected livestock serving as a major transmission source for humans. T. gondii seroprevalence in pigs varies widely across regions. Reported rates range from 0.8 to 64.4% in Europe, 0.6–82.7% in North America, 1.3–90.4% in South America, 0–27.7% in Asia, and 9.3–40.6% in Africa [19]. These variations likely result from differences in production systems, management practices, geography, and diagnostic methods, such as Enzyme-Linked ImmunoSorbent Assay (ELISA) and Modified agglutination test. In goats, seroprevalence ranges from 1.7 to 75% in Asia, 6.8–53.4% in North America, 3.7–81.8% in South America, 3.7–73.8% in Europe (excluding Turkey), and 4.6–90% in Africa [20]. For sheep, rates vary between 4 and 84.5% in Europe (excluding Turkey), 27.1–73.8% in North America, 7–66.4% in South America, 3–72.6% in Asia, and 5.6–68% in Africa [21].

In Turkey, despite these influencing factors, a review of seven studies reported an average seroprevalence of 47.2% (22.3–93.5%) in 1416 tested cattle. Six studies involving 1048 goats found a mean seroprevalence of 56.2% (12.9–95.2%). Meanwhile, 23 studies on 6522 sheep showed an average seroprevalence of 31.5% (9.5–100%) [22]. Toxoplasmosis in livestock primarily affects reproduction, leading to abortions, stillbirths, and neonatal deaths. A meta-analysis estimated T. gondii prevalence in aborted fetuses and stillbirths at 42% in sheep, 31% in goats, and 19.3% in cattle [23].

A comprehensive strategy for controlling this foodborne zoonosis should integrate traditional disease control measures with host genetic profiling. Where feasible, selective breeding for T. gondii resistance could be a valuable strategy. This approach, aligned with the One Health framework, would support sustainable livestock breeding, enhance food safety, and reduce drug residues. Ultimately, it would improve the health of humans, animals, and the environment.

Few studies on humans and model organisms have explored host genetic mechanisms against T. gondii infection, but they have yielded key insights. For example, IFN-γ-dependent immunity plays a crucial role in controlling infection, with distinct genetic factors influencing resistance to acute Toxoplasmosis and Toxoplasmic encephalitis. Genetic variation affects susceptibility to ocular Toxoplasmosis in mice, while rats lacking inducible nitric oxide synthase resist infection through ROS-mediated pyroptosis. The Toxo1 locus regulates resistance in LEW rats, and co-evolutionary dynamics have shaped specialized immune responses in rodents, involving TLR11, TLR12, and Immunity-Related GTPases—elements absent in humans [24, 25]. However, no studies have examined host genetic factors in T. gondii infection in livestock species.

This study aimed to identify host genetic factors influencing T. gondii resistance in naturally infected sheep. To achieve this, sheep were serologically tested for T. gondii infection and genotyped using 50 K single nucleotide polymorphisms (SNPs) distributed across the ovine genome, enabling genome-wide association (GWA) studies. The most commonly used statistical approach in GWA studies is the Mixed Linear Model (MLM), which effectively controls for population stratification and inflation from polygenic backgrounds However, MLM-based single-marker GWA studies may fail to detect multiple loci that simultaneously influence the traits of interest due to their one-dimensional genome scanning, making them less effective for identifying polygenic effects. Another limitation of MLM-based GWA studies is the challenge of multiple test corrections for statistical significance thresholds. The widely used Bonferroni correction is often overly conservative, potentially excluding important markers by setting excessively stringent thresholds [26,27,28].

To overcome these limitations, various multi-locus approaches have been developed, including the Multi-Locus Random-SNP-Effect Mixed Linear Model (mrMLM) [27], Fast Multi-Locus Random-SNP-Effect Mixed Linear Model (FASTmrMLM) [29], Iterative Sure Independence Screening EM-Bayesian LASSO (ISIS EM-BLASSO) [30], fast multi-locus random-SNP-effect EMMA (FASTmrEMMA) [26], polygenic-background-control-based least angle regression plus empirical Bayes (pLARmEB) [31], and polygenic-background-control-based Kruskal-Wallis test plus empirical Bayes (pKWmEB) [32]. These multi-locus approaches rely on multidimensional scanning and employ a two-stage detection process. In the first stage, markers passing quality control are screened using a single-locus approach with a less stringent significance threshold (typically p < 0.005 or 0.001). Putative markers that meet this threshold are then selected for further analysis. In the second stage, these selected markers are simultaneously evaluated using multi-locus GWA models. Due to their multi-locus nature, these approaches do not require Bonferroni correction for multiple testing [26, 27, 30,31,32].

Since most complex traits have polygenic effects, this study applied four multi-locus GWA methods. These included three parametric models (mrMLM, FASTmrMLM, and ISIS EM-BLASSO) and one non-parametric model (pKWmEB) to explore genetic influences on individual antibody responses to T. gondii. The ultimate goal was to identify genetic markers associated with innate resistance/susceptibility to Toxoplasmosis and provide genomic insights for marker-assisted selection (MAS) as part of integrated control strategies against Toxoplasmosis in sheep.

This report presents the first investigation into host genetic factors influencing Toxoplasmosis in livestock, using sheep as a model. It evaluates the potential of MAS application while considering broader strategies from public health and food safety perspectives.

Methods

Animals

The animal material used in this study, Karacabey Merino lambs, belongs to the Sheep Breeding Research Institute (SBRI), an official institution under the Turkish Ministry of Agriculture and Forestry. All necessary official permissions were obtained from the relevant institution prior to the study. The study was conducted during the summer of 2023. Blood specimens were obtained via aseptic venipuncture of the jugular vein from each selected animal. Specifically, two blood samples, measuring 5–6 mL in volume, were collected for DNA isolation and serological testing. DNA extraction was performed using the commercial kit (QIAamp^® DNA Blood Mini Kit, QIAGEN, Hilden, Germany). Following blood procurement, the animals were maintained under standardized husbandry protocols, ensuring continuity of routine care and nutritional management, and no euthanasia procedures were performed.

Karacabey Merino, recognized as one of the most popular commercial sheep breeds in Turkey, underwent improvement initiatives in the 1940s through the strategic backcrossing of German Mutton Merino rams with native Kivircik ewes. Following these breeding efforts, the Karacabey Merino population has remained closed to backcrossing for over four decades. Serving as the central breeding facility, SBRI has been instrumental in the dissemination of Karacabey Merino throughout Turkey since the 1940s. While specific data regarding Toxoplasmosis prevalence within SBRI flocks was unavailable, the uncontrolled ingress of numerous stray cats into the premises suggests the potential for a certain prevalence. Consequently, the observation of a notable Toxoplasmosis prevalence would not be unexpected in such circumstances.

The present study was meticulously designed to ensure the acquisition of robust genotype-phenotype association results, beginning with the selection of samples. Specifically, all sheep were chosen from identical shelters and subjected to uniform management practices, reflective of their birth and rearing locations. Furthermore, to standardize exposure duration, all selected sheep were of similar age, ranging between 18 and 20 months. This approach aimed to ensure uniformity in the presumed loads of T. gondii oocysts within the shared environment, encompassing factors such as soil, manure, and drinking water. Given the homogeneity in age across the selected sheep, it was presumed that the duration of exposure to T. gondii oocysts was comparable among individuals. Consequently, each sheep selected for the study was deemed to bear an equal risk of infection, considering both the intensity and duration of exposure.

Serological test

To determine antibody response to T. gondii, a total of 760 serum samples were tested with indirect Enzyme-Linked ImmunoSorbent Assay (ELISA). A commercial test kit, ID Screen^® Toxoplasmosis Indirect Multi-species (IDvet, Grabels, France), was used according to the manufacturer’s instructions. The commercial ELISA kit was prepared using the P30 (SAG1) antigen targeting T. gondii-specific IgG antibodies, which appear and persist for years. Following the prescribed protocols, ELISA microplates underwent incubation and subsequent washing steps, after which optical densities (OD) for each sample were measured at a wavelength of 450 nm.

The Sample-to-Positive ratio (S/P) was computed using the following formula:

$$\displaylines{ S/P = \cr \frac{{OD\left( {sample} \right) - OD\left( {mean\:of\:the\:negative\:controls} \right)}}{{OD\left( {mean\:of\:the\:positive\:controls} \right) - OD\left( {mean\:of\:the\:negative\:control} \right)}} \cr \times 100 \cr} $$

As per the manufacturer’s manual, serum samples with S/P % values ≤ 40% were classified as negative, while those with 40% < S/P% < 50% were considered doubtful, and samples with S/P% ≥ 50% were deemed positive. To ensure the attainment of stringent genotype-phenotype association results, samples exhibiting weak antibody responses, which could potentially be interpreted as suspiciously positive, were excluded. Thus, samples falling within the 40–100% S/P range were eliminated from the analysis. Subsequently, samples with S/P % < 40% were categorized as seronegative, while those with S/P% > 100% were regarded as seropositive. Consequently, a total of 683 sheep (comprising 321 rams and 362 ewes) were retained for further analysis.

Genotyping and quality controls

The animals underwent genotyping using the Ovine 50 K BeadChip on the Illumina platform at a private laboratory (Gen Era Diagnostic, İstanbul, Türkiye), achieving a total genotyping rate of 0.995. Initially, genotype data for sex chromosomes were excluded from the analysis. Subsequently, data quality control procedures were conducted utilizing Plink 1.07 [33]. The quality control (QC) parameters employed were as follows: Minor Allele Frequency (MAF) > 0.05, animal missing genotype rates (mind) < 0.1, and Hardy-Weinberg equilibrium threshold > 0.00001. Upon completion of QC, a total of 36,529 SNPs and 701 sheep met the predetermined criterion thresholds. Following the application of the aforementioned ELISA S/P criteria (seronegative group; S/P % < 40% and seropositive group; S/P% > 100% ), 683 sheep were selected for further genotype-phenotype association analyses.

Association study

Four multi-locus GWA approaches (mrMLM, FASTmrMLM, ISIS EM-BLASSO, and pKWmEB) were employed to identify markers associated with the antibody response against T. gondii. For parametric multi-locus methods (mrMLM, FASTmrMLM, ISIS EM-BLASSO), the individual sample-to-positive (S/P) ratio was utilized as the trait. The normality test of S/P values failed; therefore, Box-Cox and logarithmic transformations were separately tested to fit the data to a normal distribution. The Box-Cox transformation provided the best fit for normalizing the S/P data. Conversely, for the nonparametric multi-locus method (pKWmEB), binary formatted (case-control) S/P values were used as the dichotomous trait [32]. A multi-locus GWA analysis was performed using the default parameters with mrMLM R package [34]. Kinship matrices were generated via the mrMLM software and incorporated into all four multi-locus models to control for cryptic relatedness. To address population stratification, Principal Component Analysis (PCA) was conducted using TASSEL software [35] with the first five principal components included in the statistical model. A Chi-square test was employed to assess the potential effect of animal gender on the immune response against T. gondii; consequently, sex type was added to the model as a covariate. In multi-locus GWA studies, a generally accepted threshold for the logarithm of odds (LOD) score is ≥ 3 [26,27,28, 30,31,32]. To obtain the most robust association results, we employed the LOD criterion (LOD ≥ 5) for statistical significance. Manhattan plots were generated separately for parametric and non-parametric multi-locus statistics to visualize significant SNPs (LOD ≥ 5) across chromosomes. To detect any inflation or deflation of test statistics due to systemic biases under the null model, such as unnoticed population stratification or cryptic relatedness, The Quantile-Quantile (QQ) plots were utilized.

Genomic annotation, gene ontology, and pathway analysis

Genomic annotation was performed using the Ensembl Genome Browser online database (https://www.ensembl.org/) using the Oar_v3.1 assembly. The genes located within ± 100 Kb proximity of each associated SNP underwent further analysis to determine the pathways they are involved in. Gene ontology and pathway analysis were carried out using the popular online tool ShinyGO [36] with the Ovis aries database, employing a false discovery rate (FDR) cutoff of 0.05 for statistical significance.

Statistical power analysis

Statistical power analysis was performed using Cohen’s f² effect size metric with the “pwr” R package (v4.3.3), based on the parameters of MAF and the proportion of variance explained (r²). The significance threshold for statistical power calculation was set at p < 2.72E-4 by applying a Bonferroni correction to the 184 SNPs that passed the initial p ≤ 0.005 criterion and underwent multi-locus analysis in the second phase of the association analysis. GWAS-based genotype-phenotype association analyses typically require a statistical power of ≥ 0.8, which is widely accepted as a threshold for ensuring the reliability of results. Therefore, in this study, the threshold for statistical power of ≥ 0.8 was also applied [37, 38].

Results

The mean S/P value for the seronegative group was 17.13 (min, 1.32; max, 39.89), while for the seropositive group, it was 259.7 (min, 109.3; max, 342.2). The mean antibody titer of the seropositive group was found to be 15.2 times greater than that of the seronegatives. A detailed distribution of S/P values across male, female, and overall animals is presented as a box plot in Fig. 1. The overall seropositivity rate was determined to be 24.16%, with the distribution of seropositivity observed as 30.8% (99 out of 321) in rams and 18.2% (66 out of 362) in ewes. The Chi-square test revealed that the seropositivity rate was significantly higher in rams compared to ewes (X² = 14.8, p = 0.00012).

Fig. 1

Box plot illustrating the distribution of S/P values according to positive and negative groups, as well as male and female categories. Four multi-locus methods (mrMLM, FASTmrMLM, ISIS EM-BLASSO, and pKWmEB) were employed to identify markers associated with the host immune response to T. gondii infection. A total of 16 Toxoplasmosis-associated SNPs meeting the threshold of LOD 5 ≥ were identified. Among them, 13 SNPs were detected using parametric methods (mrMLM, FASTmrMLM, and ISIS EM-BLASSO; Table 1; Fig. 2A), while 3 SNPs were identified using nonparametric multi-locus methods (Table 1; Fig. 2B). QQ plot indicated no obvious evidence of inflation or deflation factors in the test statistics, confirming that possible population stratification and/or cryptic relatedness were effectively controlled during the GWA study (Fig. 3)

Table 1 Identified Toxoplasmosis-associated markers that Meet the LOD ≥ 5 threshold

Fig. 2

Manhattan plots illustrating parametric and nonparametric GWAS results

Fig. 3

QQ plot for GWAS results

The marker OAR1_259333832.1 was identified using three multi-locus approaches (mrMLM, FASTmrMLM, and ISIS EM-BLASSO). Additionally, markers OAR3_178922031.1 and s24197.1 were identified using two multi-locus approaches (mrMLM and FASTmrMLM), while markers oar3_OAR4_22415561 and OAR12_34325493.1 were also identified using two multi-locus approaches (mrMLM and ISIS EM-BLASSO; Fig. 4).

Fig. 4

The figure 4 illustrated along with the associated SNPs and the multi-locus GWAS methods used to detect them. The size of the nodes representing the methods is directly proportional to the number of SNPs detected by that method. Similarly, the size of the nodes representing the detected SNPs is directly proportional to the number of different methods that identified the relevant SNP

Within ± 500 Kb proximity of the associated markers, 11 known and 3 novel genes (ENSOARG00020035061, ENSOARG00020029346, and ENSOARG00020030227) were identified (Table 2).

Table 2 Genes identified within ± 500 Kb proximity to associated markers

Statistical power analysis demonstrates that the associated SNPs exhibit varying levels of statistical power. Specifically, 12 out of the 16 SNPs surpass the classic threshold of 0.8, indicating that the majority of the SNPs possess sufficient power to detect significant associations (Fig. 5).

Fig. 5

Forest plot illustrating statistical power and r² of associated SNPs. Statistical power analysis showing achieved statistical powers with current sample sizes. The x-axis represents statistical power ranging from 0.35 to 1.00. Blue dots indicate the actual power values, while red horizontal lines represent the r² size for each marker. Higher power values (closer to 1.0) suggest better capability to detect true effects, with most samples showing power values between 0.80 and 1.0, indicating adequate statistical strength

The utilization of gene ontology and pathway analysis has revealed that the identified candidate genes, (PLSCR5, EPHA3, DGKB, IL12B, CGA, WDR64, TMEM158, CLMP, and SIAE), located within a ± 100 Kb range of associated markers, predominantly participate in pathways crucial for modulating the innate immune system’s response against pathogens. These pathways include the RIG-I-like receptor signaling pathway, Th1 and Th2 cell differentiation pathways, C-type lectin receptor signaling pathway, and Phosphatidylinositol (PI) signaling system. Besides bacterial pathogens (i.e., Legionellosis and Pertussis), these genes are particularly implicated in host defense mechanisms against protozoal diseases, such as African trypanosomiasis, Leishmaniasis, Chagas disease, Amoebiasis, and Toxoplasmosis (Fig. 6).

Fig. 6

Pathways involving candidate genes located within a ± 100 kb range of associated markers. Enrichment analysis showing differentially regulated pathways and diseases. The x-axis represents fold enrichment scores (0–90), while each row represents a distinct biological pathway or disease condition. Pink dots indicate the significance level (-log10 FDR = 1.31) of the enrichment.

Discussion

This report describes the underlying host genetic factors against T. gondii in naturally infected Turkish sheep. This study used a two-step multidimensional genome scanning approach. Four multi-locus statistical methods, mrMLM, FASTmrMLM, ISIS EM-BLASSO, and pKWmEB, were applied to identify polygenic effects on immunity to Toxoplasmosis. Our findings suggest that after standardizing exposure intensity, duration, and other confounding factors, multigenic effects can partially explain variations in the immune response to Toxoplasmosis.

Serosurvey results indicate that Toxoplasmosis is common in Karacabey Merino sheep, with a significantly higher prevalence in rams compared to ewes. While limited studies have compared toxoplasmosis seroprevalence between rams and ewes, some studies reported no significant difference [22, 39]. This discrepancy may be due to differences in the management systems for rearing rams and ewes. Additionally, these studies do not clearly outline methods for standardizing confounding factors like exposure intensity and duration.

Statistical power analysis showed that 12 of the identified SNPs have high statistical power (> 0.8), suggesting a strong likelihood of detecting true associations. Four SNPs fall below the 0.8 threshold, with two specific SNPs, OAR8_53470280.1 and Chr3:30289994, showing particularly low power. This indicates that these markers may be less reliable for detecting true associations. While most SNPs associated with T. gondii infection have satisfactory statistical power, those with power below 0.8 should be carefully evaluated for further investigation, as their ability to detect significant associations may be compromised.

Among the identified candidate genes, those located within a ± 100 Kb distance from the associated SNPs were further investigated for potential links to innate immunity. The IL12B gene encodes the IL-12p40 subunit, a component of two heterodimeric cytokines, IL12 and IL23, both crucial to immune responses [40]. IL12 regulates T and natural killer (NK) cells by stimulating the production of γ-interferon, thereby enhancing their proliferation and cytotoxicity. Additionally, IL12 may counteract pro-angiogenic signals observed in the progression of malignancies [41]. IL12B has been extensively studied in immune responses and has been associated with several diseases, including psoriasis, Crohn’s disease, ulcerative colitis, mycobacterial diseases, juvenile idiopathic arthritis, Takayasu arteritis, and congenital cytomegalovirus infection [42,43,44,45,46,47,48]. IL12B has also been implicated in the pathogenesis of Crohn’s disease and ulcerative colitis [44]. An experimental T. gondii challenge in both feline immunodeficiency virus (FIV)-infected and healthy control cats demonstrated increased IL12 expression in healthy cats, suggesting a significant role for IL12 in the immune response to T. gondii [49]. Furthermore, T. gondii dense granule proteins, GRA15 and GRA24, induce the secretion of IL12, IL18, and IL1 by macrophages. This release of cytokines stimulates IFNγ secretion from peripheral blood mononuclear cells (PBMCs), which in turn controls T. gondii proliferation [50].

The PLSCR5 gene is part of the phospholipid scramblase (PLSCR) family, which is involved in the translocation of phospholipids between the two leaflets of a lipid bilayer. While direct evidence of PLSCR5’s molecular functions and its association with immunity is limited, other members of the PLSCR family have been implicated in various immune responses. For example, PLSCR1, another member of the family, has been linked to antiviral and antibacterial responses, as well as autoimmune diseases [51]. Furthermore, studies have shown a significant upregulation of PLSCR in response to pathogen-associated molecular patterns in the planarian species Dugesia japonica, indicating a potential role for PLSCR in the immune response triggered during pathogen invasion [52]. Although the biological function of the PLSCR gene family has not yet been fully understood, its expression during the gametocyte stage of Plasmodium spp., which are related to T. gondii, suggests a potential role in other stages of the parasite’s lifecycle or transitions between them. PLSCR may contribute to lipid regulation during the sexual stages. Given the diverse lipid-related functions of human PLSCR1, PLSCR in Plasmodium spp. could operate through similar mechanisms [53].

The EPHA3 gene encodes the EphA3 receptor tyrosine kinase, which has been recognized as a target for antibody immunotherapy in acute lymphoblastic leukemia. The EPHA3-specific monoclonal antibody IIIA4 binds and activates both human and mouse EPHA3 [54]. The interaction between EPHA3 and Ring Finger Protein 5 (RNF5), a protein involved in downregulating viral gene expression in Kaposi sarcoma-associated herpesvirus infection, further suggests EPHA3’s role in immune system regulation [55]. EPHA3 expression has also been observed in regions of tumor neovasculature, tumor-associated immune cells, and tumor-infiltrating cells in glioblastoma, highlighting its potential involvement in the tumor microenvironment and immune response [56]. In a gene function study, it was reported that certain genes, including EPHA3, are among the host genetic factors involved in Plasmodium sporozoite infection, the causative agent of malaria [57]. Additionally, a recent study has suggested a link between variations in the expression levels of the EPHA3 gene and tolerance to Ascaridia infections in chickens, proposing that differential gene expression could play a critical role in the host’s ability to resist parasitic infections [58].

Genetic studies have provided preliminary support for the influence of DGKB on cognitive complexity [59]. Genetic variants at specific genes, including DGKB, have been reported to be associated with a decreased glucose-stimulated insulin response. This underscores the importance of pancreatic beta cell dysfunction in the genetic susceptibility to hyperglycemia and type 2 diabetes [60, 61]. Polymorphisms within the DGKB gene have been suggested to affect various cellular processes in human CD34-positive cells [62]. These processes include differentiation, proliferation, transport, and migration ability induced by stromal cell-derived factor 1. Diacylglycerol kinases family (DGKs) modulate intracellular concentrations of two critical signaling molecules, diacylglycerol (DAG) and phosphatidic acid (PA), by catalyzing the conversion of DAG to PA. Among DGK isoforms, DGKzeta is notably detected in macrophages and dendritic cells. Studies reveal that DGKzeta deficiency disrupts immune responses, leading to reduced interleukin-12 (IL-12) and tumor necrosis factor-alpha (TNF-α) production upon Toll-like receptor (TLR) activation in both experimental and live models. This deficiency also correlates with heightened protection against endotoxin-induced shock but paradoxically increases susceptibility to T. gondii infection [63]. Moreover, the DGKB gene has been implicated in genetic resistance to Paratuberculosis in Italian goat breeds, indicating a potential role for DGKB in innate immunity [64].

The glycoprotein hormone alpha polypeptide gene (CGA) is a component of various glycoprotein hormones. These include follicle-stimulating hormone (FSH) and luteinizing hormone (LH) [65,66,67]. These hormones are crucial for reproductive processes, such as ovarian folliculogenesis and the regulation of steroidogenesis in gonadal cells [65, 66]. CGA is also necessary for the production of human chorionic gonadotropin (hCG), a vital hormone produced by trophoblast cells during early pregnancy and commonly used as a biomarker in pregnancy tests [68]. A systematic review has reported that certain hormones modulate T. gondii infection in various animal models. The CGA gene encodes components involved in the structure of various glycoprotein hormones. While it is not directly associated with innate immunity, it can be speculated that its expression levels may indirectly influence the progression of T. gondii infection in the host [69].

In several studies, the WDR64 gene has been associated with obesity [70] and hematological malignancies [71]. Other members of the WDR gene family have been linked to neurodevelopmental disorders and brain malformations [72, 73]. The relationship between the central nervous system and the immune system is profound, influencing various physiological processes and disease states. Given this interconnection, the potential association between the WDR64 gene and the progression of T. gondii infection presents a particularly compelling avenue for scientific investigation.

The TMEM158 gene has been associated with various aspects of cancer, including tumor progression, chemotherapy sensitivity, and immune infiltration. It has been shown that TMEM158 expression is downregulated in prostate cancer and is linked to disease progression and anti-tumor immune infiltration, particularly with NK cell and mast cell enrichment [74]. Its association with cancer and immune responses has led to the proposition that TMEM158 is a tumor-suppressor gene and a target gene in the mutator pathway, indicating a potential role for TMEM158 in regulating cellular processes related to tumor suppression and mutagenesis [75]. Additionally, higher expression levels of TMEM158 have been observed in individuals with severe West Nile virus infection, suggesting a potential role for TMEM158 in the immune response to viral infections [76, 77]. Moreover, studies have demonstrated the role of the TMEM gene family in some bacterial and viral diseases. For example, the role of TMEM208 in the immune response of teleost fish has been reported, with studies highlighting its critical importance in the defense mechanisms of large yellow croaker (Larimichthys crocea) against visceral white spot disease, caused by Pseudomonas plecoglossicida [78]. Additionally, several studies have consistently highlighted the link between genetic variations in the TMEM154 gene and resistance to Visna/Maedi, a progressive and fatal viral disease in small ruminants [79, 80].

The CLMP gene encodes the Coxsackie- and adenovirus receptor-like membrane protein. This gene has been implicated in various physiological and pathological processes. It is known to be associated with congenital short-bowel syndrome (CSBS) [81,82,83]. The CLMP gene plays a crucial role in cell adhesion. It is essential for maintaining normal intestinal homeostasis and development. This is achieved by stabilizing the gut vascular barrier positioned between endothelial and epithelial cells within the junctional complex. Moreover, CLMP has been identified as having potential prognostic value in colon cancer, with its expression correlating with immune infiltration and inflammatory response [84]. Additionally, CLMP has been recognized as an adipocyte adhesion molecule (ACAM) and is upregulated in obesity [85]. Furthermore, CLMP has been linked to immune cell migration in the context of multiple sclerosis [86].

The SIAE gene, encoding sialic acid acetylesterase, has been implicated in immune function and autoimmune diseases. Studies have demonstrated that functional rare genetic variations in SIAE are associated with the regulation of human immune tolerance [87]. SIAE is involved in regulating B-cell immune responses and serves as a negative regulator of B lymphocyte signaling, thereby contributing to the maintenance of tolerance in germinal centers [87, 88]. Additionally, SIAE has been linked to the regulation of immune cell signaling, as evidenced by its involvement in recruiting inhibitory sialic-acid-binding immunoglobulin-like lectins (SIGLEC) receptors to limit immune cell signaling [89].

Pathway analysis using the KEGG database based on the Ovis aries genome revealed that the identified pathways predominantly involve innate immunity, particularly in response to pathogens (Fig. 7). For instance, the RIG-I-like receptors (RLRs) are pivotal components in the surveillance of viral infections, orchestrating the transcriptional upregulation of type I interferons and various antiviral genes essential for host defense. Recent investigations have elucidated that RLR activation can be prompted by both viral and endogenous RNA, a phenomenon integral to mounting an effective antiviral response. However, dysregulated RLR activation poses the risk of immunopathological outcomes, underscoring the necessity for stringent regulation of RLR-mediated signaling cascades [90]. Furthermore, the RLR signaling pathway interacts with other immune signaling pathways, such as the Toll-like receptor signaling pathway, to enhance immune function and combat infections [91, 92].

Fig. 7

Network illustration of detected pathways related to resistance/susceptibility to Toxoplasmosis

The phosphatidylinositol (PI) signaling system, particularly through phosphatidylinositol 3-kinase (PI3K), has been implicated in innate immune cells, influencing inflammation mediators, cell adhesion, and cytokines [93]. The involvement of the PI signaling system has been associated with immune dysregulation and immunodeficiency [94]. Moreover, it is suggested that the PI signaling system may be linked to immune-related pathways such as WNT, MAPK, TGF-β, and JAK-STAT [95].

Th1 and Th2 cell differentiation constitute another cascade orchestrating the immune response. Upon encountering antigens, naive CD4 + T cells undergo differentiation into distinct subsets of effector T helper (Th) cells, a process regulated by the cytokine milieu during activation. These Th cell subsets, including Th1 and Th2, elicit diverse immune responses. Th1 cells, characterized by their interferon-gamma (IFN-γ) production, play a critical role in immunity against intracellular pathogens. Conversely, Th2 cells, which produce interleukin 4 (IL-4), IL-5, and IL-13, are pivotal in mediating immunity against parasitic worms and are also implicated in allergic inflammation [96].

The C-type lectin receptor (CLR) signaling pathway is a crucial component of the immune system, playing a significant role in the recognition of a wide range of pathogens and the initiation of immune responses. CLRs constitute a diverse group of receptors that recognize carbohydrate structures on various pathogens, including viruses, bacteria, parasites, and fungi [97, 98]. More than fifty CLRs have been identified, and they are functionally diverse, playing essential roles in processes such as phagocytosis, endocytosis, pathogen recognition, and complement activation [99, 100].

Numerous studies have investigated host genetic factors in Toxoplasmosis, focusing on humans or other model organisms rather than farm animals. Studies have recognized for over two decades that IFN-γ-dependent cell-mediated immunity is crucial. It plays an important role in controlling T. gondii infection and preventing Toxoplasmic Encephalitis (TE). IFN-γ is produced by various T and non-T cells, and its production is enhanced by IL-12, IL-18, Bcl-3, NF-κB [2], and the CD40-CD40L interaction. Regulation of IFN-γ-mediated immune responses is essential to prevent immunopathology resulting from overstimulation. Host genes determine resistance/susceptibility to infection, with distinct genes involved in acute infection resistance versus TE development. In mice, five genes govern acute stage survival. The L^d gene within the H-2 complex confers TE resistance. Additionally, in AIDS patients with TE, HLA-DQ3 is associated with susceptibility, while HLA-DQ1 is linked to resistance [24]. In a study, immune mechanisms in Ocular Toxoplasmosis (OT) were investigated by comparing C57BL/6 (B6) and BALB/c mice to assess genetic susceptibility. After T. gondii infection, B6 mice exhibited more severe ocular pathology, higher parasite loads, and upregulated immune response genes, including triggering receptor expressed on myeloid cells 1 (TREM-1), interleukin (IL)-1β, IL-33/ST2, and Toll-like receptors (TLRs), alongside increased neutrophil infiltration in ocular tissues. The study highlights how these immune factors correlate with OT pathogenesis across genetic backgrounds, emphasizing the critical influence of host genetics on disease susceptibility [101].

It has been demonstrated that the absence of inducible nitric oxide synthase (iNOS) confers complete resistance to T. gondii infection in iNOS knockout (iNOS−/−) rats. When iNOS−/− rat peritoneal macrophages (PMs) are infected with T. gondii, elevated levels of reactive oxygen species (ROS) are generated. These ROS are primarily induced by T. gondii-secreted GRA43. Furthermore, it has also been found that the increased ROS levels lead to the depletion of mitochondrial superoxide dismutase 2 and trigger PM pyroptosis and subsequent cell death [102].

Cavailles et al. (2014) reported that natural resistance to Toxoplasmosis in the LEW rat model is governed by the Toxo1 locus on chromosome 10, which directs parasite control via macrophage-dependent mechanisms. Toxo1 is located within a conserved haplotype block in resistant rat strains. It harbors polymorphisms, particularly in the NLRP1A gene. These polymorphisms are associated with the activation of pyroptosis as a defensive mechanism. Macrophages bearing Toxo1 alleles exhibit pyroptosis-like features upon T. gondii infection, involving caspase-1 activation and IL-1β secretion [103].

Gazzinelli et al. (2014) highlighted that the co-evolutionary dynamics between T. gondii and rodents have led to complex innate immune responses, characterized by the emergence of specialized molecular components such as Toll-like receptors (TLR11, TLR12), and an expanded repertoire of Immunity-Related GTPases (IRGs), which play pivotal roles in host resistance. Murine hosts rely on TLR11-, TLR12-, and IRG-dependent pathways to control T. gondii infection. Conversely, these elements are absent in humans, raising questions about human innate immunity to the pathogen. The lack of functional TLR11, TLR12, and IRG genes within the human genome suggests that humans, and possibly other vertebrates, lack co-evolutionary pressure with T. gondii, preventing the development of these specialized defenses [25].

Genetic factors influencing resistance or susceptibility to T. gondii infection exhibit significant variation across species. Previous studies in humans and model organisms have primarily focused on candidate gene approaches to investigate immune responses. The absence of GWAS in other species limits our ability to compare gene functions across taxa, restricting deeper insights into genetic mechanisms and their relevance across species. Our study presents the first GWAS in livestock, identifying 16 Toxoplasmosis assaciated SNPs. The observed divergences between findings in humans, model organisms, and sheep likely arise from differences in genotyping techniques and species-specific adaptations to the parasite. These findings underscore the need for expanded cross-species genetic research. Future GWAS across diverse farm animals may uncover unique mechanisms of resistance or susceptibility, advancing our understanding of evolutionary adaptations to T. gondii and informing species-specific control strategies within the One Health framework.

Limitations and concerns

Confounding factors

Despite controlling for confounding factors such as breed, age, location, and husbandry practices in the sampled material, the absence of data on T. gondii oocyst concentration in soil, animal bedding, and consumed water can be considered as one of the limitations of this study.

Stringent significance threshold

In this study, a LOD score ≥ 5 was chosen instead of the generally accepted minimum significance threshold of LOD score ≥ 3. While this stricter threshold significantly reduces the false-positive rate, it also increases the potential for missing markers with smaller effects. Additionally, each multi-locus model tested inherently detects associated markers with varying LOD scores. This naturally leads to a reduction in the co-detection rate of markers that have LOD scores within the 3–5 range across different multi-locus methods.

Statistical power of detected markers

Among the 16 markers identified, four of them had statistical power values below the 0.8 threshold. Caution should therefore be taken when interpreting these markers, as their ability to detect true associations may be compromised.

Exclusion of samples with intermediate ELISA S/P values

Concerns may arise regarding potential bias due to the exclusion of samples with ELISA Sample-to-Positive (S/P) values between 40 and 100. By removing samples close to the positivity threshold (S/P = 50), we ensured a minimum 2.5-fold difference in antibody titers between the positive and negative groups. This decision was based on a scientific cost-benefit analysis, as including false-positive or false-negative samples in the statistical analyses was anticipated to pose greater risks to the reliability of the results.

Need for validation and functional studies

The results reported in this study need to be retested in different populations before being applied in Marker-Assisted Selection programs. Furthermore, functional studies, such as gene expression analysis or CRISPR-based gene function studies, are required to confirm the roles of the candidate genes identified here.

In conclusion, the management of Toxoplasmosis in farm animals is crucial for safeguarding public health, ensuring food safety, and upholding animal welfare standards. Selective breeding practices may complement existing control measures, offering an additional strategy for managing the disease. Our comprehensive multi-locus genome-wide association analysis identified 16 significant SNPs, 12 of which demonstrated sufficient statistical power. These SNPs appear to play a notable role in influencing the host immune response to T. gondii infection. The findings from this study provide new insights into the immune pathways related to resistance and susceptibility to Toxoplasmosis in sheep. These insights could aid future MAS application efforts. This research underscores the interconnectedness of human, animal, and environmental health in the context of infectious disease management, aligning with the One Health approach.