A Genome-Wide Association Study for Resistance to Tropical Theileriosis in Two Bovine Portuguese Autochthonous Breeds

1. Introduction

Tropical Theileriosis is a tick-borne disease caused by a haemoprotozoan, Theileria annulata, which affects cattle production in Europe, Asia, and Africa [1]. This parasite can affect the host’s immune system and cause damage to red blood cells, leading to varying degrees of anemia [2]. This will result in production shortfalls and increased mortality and morbidity rates of parasitized animals. Thus, tropical theileriosis causes economic losses associated with reduced production capacity and the need to carry out control, for example, through the use of acaricides [3].

Theileriosis control strategies may include curing affected animals and prevention [4]. Curing affected animals can be difficult and involves the administration of substances such as imidocarb, erythromycin, oxytetracycline, quinozoline, and naphthoquinone derivatives (parvaquone and buparvaquone), whose efficacy is questionable [5,6,7]. In addition, buparvaquone, for example, can leave residues in animal products such as meat and milk, implying long safety intervals, and its use is therefore limited in several countries, including European countries, where its use is not approved by the EMA (European Medicines Agency) [5,8]. In an attempt at control, where it is difficult to achieve a cure, it is possible to resort to culling affected animals, although this results in significant economic losses for the farm [4]. As for prevention, this can be based on strategies that include control of ticks (e.g. use of acaricides), the environment (e.g. biosecurity measures and avoiding movement of animals from non-endemic to endemic areas), and the animal (e.g. selection of resistant animals and use of vaccines) [4,5]. The use of acaricides in the control of ticks and tick-borne diseases is one of the main strategies used, but it is costly [4]. In addition, increasing acaricide resistance is also a concern in the implementation of this control strategy [4,7,9]. Acaricide resistance results from the selection of specific heritable traits in a tick population, resulting from exposure of the population to an acaricide. As a result, there will be an increase in the number of ticks that survive after administration of the recommended dose of the same acaricide to which resistance already exists [9]. In turn, the use of vaccines, an animal-based control strategy, also has some limitations. The only vaccines that appear to be effective are live attenuated vaccines, but more research is needed to fully understand their mechanism of action [5,7]. Subunit vaccines are difficult to obtain, considering the genetic diversity of these parasites [5]. Finally, the use of resistant animals could be another, more sustainable, strategy to minimize the impact of the disease [6,8].

The use of animals that are genetically tolerant/resistant to ticks and tick-borne diseases is thus considered the most natural control alternative [10]. Natural disease resistance refers to the inherent ability of an animal to resist diseases when exposed to pathogens without prior exposure or immunization [4]. Several authors have reported that cattle of autochthonous breeds from endemic regions are more resistant than those of exotic breeds [3,11,12,13]. For example, different responses to T. annulata infection have been identified in the Sahiwal breed (B. indicus) and the Holstein breed (B. taurus) [11,14].

An important prerequisite in breeding strategies for disease tolerance/resistance is the attempt to understand the genetic control mechanisms of diseases. Nowadays, it is known that tick loads in cattle have moderate to high levels of heritability, which can range from 0.40 to 0.54 [15,16]. In addition, some haemoparasitoses such as those caused by Anaplasma marginale and Ehrlichia ruminantium appear to have significant heritabilities, with values of 0.16 and 0.19, respectively [17]. Genome-wide association studies have been conducted for some vector-borne diseases [10]. In addition, single nucleotide polymorphisms (SNPs) have been tested to identify genetic variants associated with complex traits. SNPs distributed throughout the genome can detect and map mutations underlying variation in target traits by a process called genome-wide association analysis (GWAS) [18]. GWAS allows to identification of genetic markers, candidate genes, and QTLs for individual traits. The discovery of Quantitative Trait Locus (QTL) is an important step in identifying and understanding genetic variants associated with economically important phenotypes, and genome-wide association study (GWAS) has become a widely used approach for identifying QTLs and genome regions associated with phenotypes [19]. High-density arrays capable of genotyping thousands of SNPs are already available for cattle, allowing for increased genomic coverage and statistical power [20].

In this study, we performed a GWAS using a high-density array to identify SNP genetic markers for resistance to Tropical Theileriosis in cattle of the Portuguese autochthonous breeds Alentejana and Mertolenga and to identify the genes or genomic regions that are most likely involved in tolerance/resistance to Tropical Theileriosis. Both breeds are officially at risk of extinction and have conservation and improvement programs. The data obtained in this work could be important information to be included in these breed improvement programs.

2. Materials and Methods

2.1. Sample collection and selection of phenotypically extreme animals

The blood samples used in this study were collected in the Alentejo region of Portugal, the country’s leading beef producing region, with more than 65% of national production. This region is the place of origin of the Portuguese autochthonous breeds under study, namely the Alentejana and the Mertolenga breeds, and is also an endemic region for T. annulata [21]. The Alentejana and Mertolenga breeds belong to a group of 15 autochthonous Portuguese cattle breeds and are at risk of extinction [22]. The Alentejana cattle breed has around 22 000 breeding females registered in the herd book, spread over around 174 farms, although only around 8 000 are purebred, and the rest are used for crossbreeding with males from exotic international breeds [23]. On the other hand, the Mertolenga cattle breed, which currently has 27 000 breeding females, is the largest of Portugal’s 15 autochthonous breeds. As with the Alentejana breed, only around 8 000 females are purebred, the rest being used to crossbreed with males from exotic international breeds [24].

Initially, 843 blood samples were randomly collected from purebred Alentejana and Mertolenga cattle, which did not show any clinical signs of T. annulata infection. The collections were carried out between November 2018 and December 2019, in 420 Alentejana breed animals and 423 Mertolenga breed animals. These animals all belong to farms with extensive or semi-extensive production regimes, so these animals live outdoors and on pasture.

Between 3 and 5 mL of blood were collected from the jugular vein of each animal and stored individually in tubes containing EDTA (ethylenediaminetetraacetic acid). This blood sampling was carried out by the technicians of the cattle associations of the autochthonous breeds under study. The tubes containing the blood sample were frozen at -20°C before being sent to the laboratory. During transport and storage, the temperature conditions were maintained.

Initially, these samples were analyzed at the Molecular Genetics Laboratory of the National Institute for Agricultural and Veterinary Research (INIAV). In the laboratory, 300 µL of each of the 843 blood samples were used to perform DNA extraction using the Cytogene® Blood Kit (India, Cytogene), following the manufacturer’s instructions. Thereafter, DNA from each sample was subjected to amplification with Polymerase Chain Reaction (PCR) of a fragment, with about 319 base pairs (bp), of a T. annulata merozoite-piroplasm surface antigen gene, Tams 1 [25]. In all reactions, were used a negative sample of T. annulata (negative control), without DNA, and a positive sample (positive control), property of the Parasitology Laboratory of the INIAV. The amplified samples were analyzed by 1.5% agarose gel electrophoresis, and the gel was visualized with an ultraviolet (UV) transilluminator. Positive and negative PCR controls were introduced in each PCR and a molecular weight marker (NZYDNA VI) was placed on each gel. Samples were classified as positive (samples infected with T. annulata, with the presence of Tams 1 gene) and negative (samples not infected with T. annulata, without the presence of the Tams 1 gene). From the 843 blood samples analyzed, 96 samples from animals of the Alentejana breed and 96 samples from animals of the Mertolenga breed were selected. In the case of the animals of the Alentejana breed, all the samples positive for T. annulata (30 samples) were used and the farms to which they belonged were identified. The remaining 66 negative samples were selected from the same farms as the previous ones. In the case of the Mertolenga breed, we used 48 samples from infected animals and 48 samples from non-infected animals (n = 96), selecting the latter based on the farm to which the infected animals belonged. The presence of the Tams 1 gene (infected animals) or its absence (non-infected animals) were the phenotypic extremes considered in our study. Thus, infected animals, with the presence of the Tams 1 gene in the blood sample, were considered susceptible to Tropical Theileriosis, while non-infected animals, without the presence of the Tams 1 gene in the blood sample, were considered resistant to Tropical Theileriosis.

2.2. SNP genotyping and GWAS analysis

All 192 samples were genotyped in two plates of 96 samples of the Axiom™ Bovine Genotyping 100K Array, one plate for each breed, in a laboratory providing animal health and food safety services (Segalab, Portugal), following the best practices workflow from the manufacturer. The “BestandRecommended” SNPs were selected for the analysis of both breeds. PLINK v1.9 was used for quality control [26]. SNPs were filtered out whenever the minor allele frequency (--maf) was below 5% and missing data (--geno) was higher than 10%. Similarly, samples were filtered out when there was 10% or more missing data (--mind). To account for population structure, samples were filtered based on cryptic relatedness (--min 0.2), and their distribution was inferred by means of a principal component analysis (--pca var-wts). Whenever a principal component was responsible for a significative separation of samples into groups, which was not expected from sample records, it was used as a covariate for the association analysis. Genome-wide association analysis was performed in PLINK under the logistic model for disease traits (--logistic). Case/control phenotypes for T. annulata were included in the phenotypes files (--pheno). The origin of each sample, namely its breeder or the district of provenance, was used as a covariate (--covar). The association was evaluated using quantile-quantile plots. The p-values for all SNPs were adjusted for false discovery rate (FDR). Also, the threshold of significance was calculated following the Bonferroni correction.

2.3. Data analysis

Data analysis was supported by querying the Bovine Mine and National Center for Biotechnology Information databases [27,28]. When analyzing the data using these databases, the aim was to identify the overlap of suggestive SNPs with possible genes and QTLs or to identify genes in their vicinity.

3. Results

3.1. Descriptive statistics

In the sample initially used (n = 843), it was possible to find positivity to T. annulata of 7.10% (30/420) in Alentejana animals and 14.4% (61/423) in Mertolenga animals [21]. As mentioned above, we selected 96 Alentejana animals and 96 Mertolenga animals, whose characteristics are described in the table below (Table 1).

3.2. Genome-wide associations

To identify SNP genetic markers and QTLs for resistance to Tropical Theileriosis, we performed a GWAS using a genotyping approach. After implementing data quality control measures, out of the 100 000 SNPs on the array, 75 126 SNPs were used in the association test for the Alentejana breed and 81 357 SNPs for the Mertolenga breed, and the analysis was performed for each breed separately.

Considering an initial p-value ≤ 0.05, it was possible to find 7 833 significant SNPs in the Alentejana breed, and 7 157 significant SNPs in the Mertolenga breed. All SNPs were classified as non-significant based on the adjusted p-values after FDR. Furthermore, the Bonferroni correction also established a threshold of significance lower than the lowest p-value found. Despite this, based on other GWAS works available in the literature, it was decided that the suggestive value of significance to be used in this work would be p-value ≤ 10^-4 (Figure 1) [29,30,31]. Thus, it was possible to find 24 significant SNPs for Alentejana cattle and 20 significant SNPs for Mertolenga cattle (Table 2). In this table, we can see that for both breeds, the protective allele of the identified SNPs appears in a higher percentage in the animals under study. Thus, we found that the protective allele appears more frequently in 75% of the SNPs studied in the Alentejana breed and 80% of the SNPs studied in the Mertolenga breed.

For the Alentejana breed, considering the 24 SNPs suggestive of genomic significance for resistance to Tropical Theileriosis, it was not possible to find any overlapping QTLs. On the other hand, in the case of the Mertolenga breed, considering the 20 SNPs suggestive of genomic significance, 7 QTLs were found in the vicinity of the suggestive SNPs, distributed over 4 chromosomes (chromosomes 3, 7, 8, and 29) (Table 3).

In addition to the identification of QTLs based on SNPs suggestive of genomic significance, the genes where these SNPs overlapped with or the genes upstream and downstream of them, in the case of intergenic variants, were also considered. Thus, it was possible to identify different genes on chromosomes 2, 4, 6, 9, 12, 15, 16, 17, 20, 22, 25, and 26 for the Alentejana breed (Table 4), and on chromosomes 3, 4, 7, 8, 10, 11, 12, 13, 21, and 29 for the Mertolenga breed (Table 5).

4. Discussion

In this work, GWAS was performed using a high-density bovine SNP’s array that allowed the identification of 24 significant SNP’s for Alentejana breed cattle and 20 significant SNP’s for Mertolenga breed cattle (p ≤ 10^-4). For these SNPs, in both breeds, we found that the protective allele is the most prevalent. In addition, only in the case of the Mertolenga breed it was possible to identify 7 QTLs already described, associated with SNPs suggestive of resistance/susceptibility to Theileriosis. These QTLs are associated with traits such as average daily gain, muscle carnosine, creatine and creatinine content, body weight at birth, larger muscle area, and percentage of glycated kappa-casein in milk. These data are extremely important because the genetic selection of animals may consider important health traits but cannot neglect the productive traits of these animals. Thus, it is important to use these QTLs to enhance genetic and economic gains by incorporating the information obtained into breeding programs. Thus, these programs may allow the selection animals with higher combined economic value for the next generation by combining productive and non-productive traits, such as resistance to Theileriosis [36].

It was also possible to identify several annotated genes that overlap with the suggestive SNPs from GWAS, or in their vicinity. For the Alentejana breed, it was possible to find genes associated with the regulation of cell proliferation, differentiation and survival, such as EGFR (epidermal growth factor receptor), cell growth such as DCHS2 (dachsous cadherin-related 2), and the formation of plasma cation channels such as TRPC3 (transient receptor potential cation channel subfamily C member 3) [37,38,39]. In addition to these functions at the cellular level, of the genes described, the MAP3K1 (mitogen-activated protein kinase kinase 1) gene was identified, which acts in the MAPK pathway, a signal transduction pathway that modulates physiological and pathophysiological cellular responses. In this pathway, mitogen-activated protein kinases regulate important cellular processes such as proliferation, stress response, apoptosis, and immune defense, regulating the production of T helper 1 (Th1) and T helper 2 (Th2) lymphocyte responses. Currently, the ability of protozoan parasites such as Trypanosoma cruzi, Trypanosoma congolense and Leishmania spp. to modulate the host immune response by intervening in the MAPK pathway to favor their replication and survival has been described [40,41]. Another gene recognized was CMTM7 (KLF like MARVEL transmembrane domain containing 7), which belongs to the superfamily encoding chemokine like factors. Lack of CMTM7 has already been shown to cause a reduction in the innate B-cell population and lead to natural IgM and IL-10 deficiency [42]. IL 10 is one of the interleukins present in the highest concentration when cattle are infected with T. annulata and is essential in the development of the immune response against this agent [43]. In turn, the SSFA2 (sperm specific antigen 2) gene was identified, whose presence has already been reported to be associated with greater resistance to the development of clinical mastitis in cattle [44]. In addition to all these, the ATG13 (autophagy related 13) gene was also identified, as responsible for cellular autophagy, i.e. programmed cell death when cells are aged, degenerated, or non-functional. The action of this gene in the autophagy process in cells infected with the Bovine Viral Diarrhea Virus (BVDv) has already been described [45]. BVDv is an intracellular pathogen, at a certain stage of its cycle, like T. annulata [46].

Interestingly, two genes were identified in the group of Alentejana animals that seem to be associated with neurotransmitter concentration. It is known that the abnormal concentration of neurotransmitters is one of the factors that affect the health status, temperament, and welfare of animals, but the genetic basis of this abnormality is still unknown [47]. Despite this, there are already references that the PCDH15 gene (protocadherin related 15), which was identified in this work, may be associated with the regulation of this concentration. On the other hand, it was also possible to identify the LOC107133268 gene (protocadherin-16-like), whose functions are not yet fully described in cattle, but the homolog, in humans, seems to be involved in the modulation of synaptic transmission and the generation of specific synaptic connections [48]. In addition, the genes LOC107133203 (olfactory receptor 1S1-like) and LOC104974344 (olfactory receptor 10Q1), which encode olfactory receptor proteins responsible for the recognition and transduction of olfactory signals, have also been identified and can trigger a neuronal response that triggers the perception of an odor [49].

In addition, were identified genes potentially associated with resistance to Theileriosis in Alentejana breed animals that appear to have significant productive importance, such as the ATP12A (ATPase H/K transporting non-gastric alpha2 subunit), SHISA9 (shisa family member 9) and FSTL5 (follistatin like 5) genes. The ATP12A gene encodes the H+/K+ ATPase Type 2 protein, a membrane protein involved in transmembrane cation transport, which is present in sperm and acts in the acrosome reaction at fertilization [50]. On the other hand, the SHISA9 gene is associated with pre-weaning growth, while the FSTL5 gene appears to be associated with muscle hypertrophy in cattle, inducing follistin to increase insulin action at the skeletal muscle level [51,52,53,54]. Finally, the BRINP3 (BMP/retinoic acid inducible neural specific 3) gene was also identified, for which its association with meat quality and fertility in heifers is reported [55,56].

Regarding the genes identified, there is an overlap of two SNPs with the SHISA9 and EGFR genes, which affect the animals’ productive capacity and cell regulation, respectively. In addition, one SNP overlaps with the DCHS2, MAP3K1, CMTM7, PCDH15, SSFA2, and ATG13 genes, most of which act in cell regulation and protective cell response. All other genes described are near the SNPs identified.

In the case of the Mertolenga breed, it was also possible to identify genes with different functions. Thus, the UOX (Bos taurus urate oxidase) gene was identified, which encodes the urate oxidase enzyme, which has the function of degrading uric acid, a hematological parameter that is increased in animals infected with T. annulata [57,58]. In addition, the TMC2 (transmembrane channel like 2) gene, which encodes proteins responsible for the formation of mechanosensitive ion channels at the tips of sensory cells in the inner ear, has been identified in mammals [59]. In humans, it is reported to be associated with hearing loss and epidermodysplasia verruciformis [60,61]. Cumulatively, the ZBTB44 (zinc finger and BTB domain containing 44) gene has been identified that appears to be associated with the macrophage-dependent immune response in Mycobacterium avium subsp. paratuberculosis infection in cattle [62].

Regarding genes associated with production, it was possible to identify the EED (embryonic ectoderm development) gene, whose association with the milk production capacity of cattle is already been reported [63]. In addition, the LOC107131273 gene (multidrug resistance-associated protein 4-like) was identified, for which there is an indication of differential expression in cows pregnant with large fetus, which lead to the development of dystocia, with significant impacts on production [64]. Furthermore, the MTNR1B (melatonin receptor 1B) gene has been identified, which is expressed in mammalian oocytes, and the TRNAC-GCA (tRNA-Cys) gene appears to be associated with sperm quality [65]. In addition to these, the DRD1 (dopamine receptor D1) gene, which encodes the dopamine D1 receptor protein, has been identified. These receptors are prime candidates in the regulation of energy for the maintenance of homeostasis and are implicated in the regulation of feeding behavior in cattle [66]. Finally, the MEF2C (myocyte enhancer factor 2C) gene encoding a myocyte enhancer protein was identified and is important for skeletal, cardiac, and smooth muscle development [67,68].

For the Mertolenga breed, two SNPs were found to overlap with the MEF2C gene and one SNP with the LOC170131273, UOX, and ZBTB44 genes. In addition, the UOX gene is also found to overlap with a QTL (associated with average daily gain). This gene is a strong candidate for resistance to Tropical Theileriosis due to the apparent link with the pathogenesis of the disease, warranting further investigation. In addition, the fact that it is related to a QTL associated with average daily gain could also be a good indicator, since one of the losses associated with infection by T. annulata is reduced productivity [1]. All other genes recorded are near the significant SNPs identified.

In beef cattle, it is known that it is of utmost importance to increase the ability to resist diseases, which is strongly associated with their immune performance. However, the productive traits of these animals cannot be disregarded, meaning that genes associated with the production of more muscle, better quality meat, or even milk in the case of suckler females are of utmost importance.

To the best of the authors’ knowledge, there has been no other published work analyzing GWAS data from the autochthonous Portuguese Alentejana and Mertolenga breeds. As such, this is a study that contributes to deepening genomic knowledge of Portuguese genetic resources. On the other hand, the authors are also unaware of any similar work aimed at associating resistance/tolerance to Tropical Theileriosis using GWAS. Despite this, it was found that some functions of genes associated with resistance to other parasitic diseases, such as Amoebiasis, Trypanosomosis, Toxoplasmosis, and Leishmaniasis, are also common to candidate genes for resistance to Theileriosis, identified in this study. Thus, we highlighted the action on the cell membrane (CMTM7, TMC2, and ATP12A genes), ATP binding (ATP12A gene) and immune response (MAP3K1, CMTM7, SSFA2, ATG13, and ZBTB44 genes) [69]. This work then made it possible to identify some SNPs suggestive of an association with resistance/tolerance to Tropical Theileriosis and the genes and QTLs that overlap or are in the vicinity. In the future, it would be important to replicate the results presented here on a larger sample of animals, focusing on the significant SNPs.

5. Conclusions

In this work, 24 candidate SNPs for resistance to T. annulata infection were identified in the Portuguese Alentejana autochthonous breed and 20 candidate SNPs in the Portuguese Mertolenga autochthonous breed. For both breeds, the protective allele of the identified SNPs appears in a higher percentage in the animals under study. Also, 7 QTLs were found in the Mertolenga breed, of which one overlaps with the candidate gene UOX. This gene appears to be associated with the pathogenesis of Tropical Theileriosis. In the case of the Alentejana breed, the MAP3K1, CMTM7, SSFA2, and ATG13 genes will be good candidate genes for resistance to Tropical Theileriosis, due to their importance in regulating the immune response or their already described impact on resistance to other diseases. Thus, due to the importance that these genes seem to have in Tropical Theileriosis, further studies will be required, focusing on the SNPs identified in this work. On the other hand, the data found in this study could be used to define markers to be applied in breeding programs using marker-assisted selection, for both breeds under study.