Determination of Genetic Variations of Toll-Like Receptor (TLR) 2, 4, and 6 with Next-Generation Sequencing in Native Cattle Breeds of Anatolia and Holstein Friesian ^†

Abstract

In recent years, the focus of disease resistance and susceptibility studies in cattle have been on determining patterns in the innate immune response of key proteins, such as Toll-like receptors (TLR). In the bovine genome, there are 10 TLR family members and, of these, TLR2, TLR4, and TLR6 are specialized in the recognition of bacterial ligands. Indigenous cattle breeds of Anatolia have been reported to show fewer signs of clinical bacterial infections, such as bovine tuberculosis and mastitis, and it is hypothesized that this might be due to a less stringent genetic selection during breeding. In contrast, Holstein-Friesian cattle have been under strong selection for milk production, which may have resulted in greater susceptibility to diseases. To test this hypothesis, we have compared the TLR2, TLR4, and TLR6 genes of Anatolian Black (AB), East Anatolian Red (EAR), South Anatolian Red (SAR), Turkish Grey (TG), and Holstein (HOL) cattle using next-generation sequencing. The SAR breed had the most variations overall, followed by EAR, AB, TG, and HOL. TG had the most variations for TLR2, whereas SAR had the most variations in TLR4 and TLR6. We compared these variants with those associated with disease and susceptibility traits. We used exon variants to construct haplotypes, investigated shared haplotypes within breeds, and proposed candidate haplotypes for a disease resistance phenotype in Anatolian cattle breeds.

1. Introduction

Toll-like receptors (TLRs) recognize conserved patterns in diverse microbial molecules called microbial-associated molecular patterns (MAMPs). These include lipopolysaccharide recognized by TLR4 and lipopeptides recognized by the heterodimer formed by TLR2 with either TLR1 or TLR6 [1]. In addition, TLRs also react to damage-associated molecular patterns (DAMPs) [2,3] released after cellular damage and, thus, play a crucial role in initiating the innate immune response [4]. Genetic variations found in the genes encoding TLRs have been associated with disease susceptibility and resistance in a variety of animal species [1,5,6]. In the bovine genome, 10 members of the TLR family (TLR1–10) have been identified and mapped to specific chromosomes [7]. Bovine TLR2, TLR4, and TLR6 genes are located on BTA17, BTA8, and BTA6, respectively, and have defined amino acid (aa) lengths [8,9]. Independent of their aa length, all TLRs identified so far contain three domains: an extra-cellular ligand binding domain (ECD) consisting of various numbers of leucine-rich repeats (LRR), a trans-membrane domain (TMD), and an intra-cellular domain, which is also known as the Toll/Interleukin-1 Receptor (TIR) domain [10], of which the TMD and TIR domain seem to be highly conserved between species [11].

Several studies have suggested that disease resistance, susceptibility, and severity of clinical signs in individuals or breeds may be attributed to altered ligand binding caused by single nucleotide polymorphisms (SNPs) specifically within the ECD of TLR genes [11,12]. Given the increased use of antibiotics in food-producing animals, native cattle breeds could provide an important genetic resource to breeders. Many of these indigenous breeds show an enhanced disease resistance due to a co-evolution with specific pathogens over decades, potentially resulting in the development of genetic resistance. Indeed, due to low selective pressure, local breeds are generally not as productive as high-yielding breeds selected for specific QTLs and, thus, may have preserved their genetic makeup over the years.

Anatolia has five native cattle breeds; Anatolian Black (AB), South Anatolian Red (SAR), East Anatolian Red (EAR), Turkish Grey (TG), and Native Southern Yellow. In general, these breeds grow slowly, are adapted to extensive breeding conditions, and present the morphologically look of semi-wild cattle [13]. Milk production in these breeds is between 700–1000 kg per lactation period per cow. The breeds are adapted to challenging environmental conditions and are known to be resistant to specific diseases. For example, SAR and EAR breeds are known to be resistant to blood parasite infections [14], whereas cows from the AB breed showed an extremely low bovine tuberculosis incidence rate when kept with Mycobacterium bovis-infected cattle [15], as well as a low mastitis-incidence rate [15]. As the corresponding disease-causing bacteria have been described to bind to TLR2, TLR4, and TLR6, the aim of the study was to determine the presence of SNP in these genes in the AB, EAR, SAR, and TG indigenous breeds compared to those present in the corresponding genes of HOL cattle.

2. Materials and Methods

2.1.DNA Isolation and Amplicon Sequencing

To determine variations in the genomic sequences for TLR2, TLR4, and TLR6, blood samples from AB (n = 20), EAR (n = 20), SAR (n = 20), TG (n = 10), and HOL (n = 10) were used. Blood samples were collected from each breed across Turkey (Figure 1).

Genomic DNA (gDNA) extraction was performed using a commercially available gDNA extraction kit (Qiagen Blood and Tissue, Hilden, Germany). Extracted DNA was measured by spectrophotometry (Nanodrop 2000, Thermo, Wilmington, DE, USA) and integrity was visualized by agarose gel electrophoresis. To obtain reliable variant information, paired end sequencing was performed at 50× coverage to analyze TLR2, TLR4, and TLR6, which spanned 13, 11, and 19 kb regions, respectively. Different primer pair combinations were used for each gene (for details see: Table S1). Primers were designed using the Primer 3 software package [16] and were spaced out over 2500–3500 bp intervals. Hot start taq DNA polymerase (Phire Hot Start II, Thermo Fisher, Bremen, Germany) was used in PCR applications. PCR conditions were as followed; 95 °C 60 s, for 45 cycles; 95 °C 10 s, 60 °C 120 s, 72 °C 20 s, and final elongation at 72 °C 60 s. For GC-rich regions 5% DMSO was added to the reaction to enhance PCR results (for details see: Table S2). PCR amplicons were visualized using SybrSafe (Invitrogen, Paisley, Scotland, UK) stained agarose gels. NGS library preparation was performed following the manufacturer’s instructions (Nextera Library Preparation Kit, Illumina Inc., San Diego, CA, USA). Amplicons were sequenced using paired end sequencing on the MiSeq platform (Illumina in Intergen Lab Inc., Ankara, Turkey), and reads were aligned to the present bovine genome version available Btau 4.6.1 (Btau7) using the MiSeq software (Illumina, San Diego, CA, USA). The resulting alignment was used to construct binary-aligned map (bam) files.

2.2. Sequence Analysis, Variant Verification, and Protein Modelling

TLR2, TLR4, and TLR6 sequence datasets from individuals were analyzed using Picard [17], BAM tools, SAM tools [18], and GATK [19] to generate variation call files (VCF) for each individual [20]. SnpSift [21] and SnpEff analysis tools [22] were used to annotate variants with the SNP138 variant collection. Thereafter, VCFs were aggregated at the breed level to determine novel SNP and InDel variants, and the results were stored in distinct .vcf files. Positions of both SNPs and InDels were lifted over to a newer assembly version (Btau8) in UCSC [8] and new positions were uploaded to the Ensembl Variant Effect Predictor (VEP) [23] database. The accuracy and the annotation of the identified genetic variations were assessed using VEP [23] and SnpSift [21], respectively.

Although VEP [23] and SnpSift [21] are useful tools to assess accuracy and to annotate novel candidate SNPs, we also screened variants using their coverage and frequency to eliminate false variant calls. Potentially important non-synonymous and novel SNP variants located in exons were validated subsequently by using Sanger sequencing. New sets of primers covering the regions including variants were designed by using Primer3 [16] (Table S1). The PCR products were sequenced using a BigDye Terminator v3.1 cycle sequencing kit and ABI310 automatic sequencer (Applied Biosystems, Foster City, CA, USA). Amplification primers were used for bidirectional sequencing. Obtained data was analyzed using BioEdit software [24].

After analyzing the variants at the amino acid level, a protein model was constructed for the most important variant identified in TLR2. Models were constructed using the Modeller [25] software package, and validated using ProCheck [26], Verify 3D [27], ERRAT [28], and ProQ [29].

3. Results

Obtained .bam files were visualized using Integrative Genomics Viewer (IGV) [30]. The average read length was determined to be 125 bp. Some .bam files contained ambiguous alignments due to the presence of nonspecific amplicons produced during the amplification of gene regions. Thus, all ambiguous alignments and reads with less than 125 bp were removed and filtered according to their mapping quality before calling variants using MQ > 50. Across all analyzed sequences from all breeds, five intronic regions in 13 individuals could not be amplified and, therefore, NGS results are missing for these individuals. Furthermore, one individual from the HOL breed appeared to have variations found only in Anatolian breeds and was subsequently excluded from the analysis.

Within the three TLR genes analyzed, a total of 360, 463, 520, 423, and 274 SNP variants were determined at the breed level in AB, EAR, SAR, TG, and HOL individuals used in this study, including 26, 23, 68, 22, and four novel SNPs (

Table 1. SNP variants identified in TLR2, TLR4, and TLR6 genes. In each cell of the table, the total/novel number of SNPs is shown, or only the total if no novel SNPs were found.

Gene	Breed
	AB	EAR	SAR	TG	HOL
TLR2	133/11	153/17	155/13	187/16	45/3
Exon	23	27	28	26	6
missense	13	16	17/1	15/1	4
synonymous	10	11/1	11/1	11	2
Intron	93/10	111/15	110/10	146/15	34/3
3′UTR	15/1	13/1	15/1	14	4
5′UTR	2	2	2	1	1
TLR6	209/14	230/5	285/52	179/5	155
Exon	26	23	26	24	9
missense	10/2	9/1	9/1	8/1	3
synonymous	16/1	14	17/2	16/1	6
Intron	177/8	203/4	252/46	153/2	103
3′UTR	2	1	3/2	2	1
5′UTR	4/3	4	5/1	1/1	2
TLR4	18/1	80/1	80/3	57/1	74/1
Exon	2	16	13	10	15
missense	1	6	2	3/1	5
synonymous	1	10	11	7	10
Intron	14	59/1	62/3	45	55/1
3′UTR	1	2	2	1	1
5′UTR	1	3	5	1	5

). According to average data on determined SNPs and genes, the highest SNP variation was found in the TLR6 gene, whereas the lowest one was found in TLR4 (Figure 2).

The SNP variants spread in the whole genes; however, InDel variants were identified only in intronic regions, suggesting no effect on the gene function (

Table 2. InDel variants identified in TLR2, TLR4, and TLR6 genes. In each cell of the table, the total/novel number of variants is shown.

Gene	Breed
	AB	EAR	SAR	TG	HOL
TLR2	12/12	16/13	15/12	16/13	6/6
TLR6	23/9	28/14	38/24	21/13	17/3
TLR4	3/1	5/3	8/3	6/4	5/1

). Novel variants and variants potentially impacting on molecule structures identified by NGS results were subsequently confirmed using bidirectional Sanger sequencing (Figure 3).

By just analyzing SNPs occurring in the exons of TLR2, TLR4, and TLR6, a total of 33, 24, and 46 SNP variants covering the privates and the shared among breeds, respectively, were determined and according to these SNPs 36, 25, and 98 different haplotypes per corresponding TLR were constructed using ShapeIT software [31]. Obtained phased haplotypes (Table S3) were visualized in a median joining (MJ) tree using the Network program [32] (Figure 4 and Figure 5), except TLR6, due to a high number of haplotypes. Within all haplotypes, two main haplotypes were identified for TLR2 and TLR6, whereas three major haplotypes were identified for TLR4. Analyzing haplotypes by breed, we identified 36 haplotypes in TLR2 for Anatolian breeds, but only four within HOL. Interestingly however, we were unable to identify any breed specific haplotypes for TLR4 and TLR6, but we were able to identify eight and 12 haplotypes, respectively, that were shared between Anatolian breeds and HOL (Table S3).

Missense and synonymous variations were analyzed according to reference proteins for all genes and several changes were found in the analyzed genes affecting the amino acids characteristics (

Table 3. Missense and synonymous variations on protein domain level of TLR2 (according to reference sequence NP_776622.1).

Domain	TLR2 (aa)	AB	EAR	SAR	TG	HOL	SNP ID
LRR1	54..77	62 N/N	62 N/N	62 N/N	62 N/N		rs68268249
		63 E/D	63 E/D	63 E/D	63 E/D	63 E/D	rs55617172
		68 G/S	68 G/S	68 G/S	68 G/S		rs68268250
LRR3	102..125		119 W/L	119 W/L	119 W/L		rs211243949
LRR5	151..175		152 R/Q	152 R/Q		152 R/Q	rs43706434
LRR7	200..223		201 S/N			201 S/N	rs110491977
		211 I/V	211 I/V	211 I/V	211 I/V	211 I/V	rs43706433
LRR8	224..250	227 F/L	227 F/L	227 F/L	227 F/L		rs68268251
LRR9-10	251..308
LRR11	309..337				315 R/R		rs68268253
		326 H/Q	326 H/Q	326 H/Q	326 H/Q		rs68343167
LRR12	338..361	337 R/Q	337 R/Q	337 R/Q	337 R/Q		rs68343168
LRR13	362..388
LRR14	389..414	405 T/M	405 T/M	405 T/M	405 T/M		rs68268255
LRR15	415..437	417 N/S	417 N/S	417 N/S	417 N/S		rs68268256
		436 G/G	436 G/G	436 G/G	436 G/G		rs68268257
LRR16-18	438..500
LRR19	501..524	502 S/A	502 S/A	502 S/A	502 S/A		rs68268258
				530 A/A			novel
				531 G/S			novel
LRR20	533..586	544 F/F	545 F/F	546 F/F	547 F/F	544 F/F	rs68268259
		563 R/H	563 R/H	563 R/H	563 R/H		rs68268260
		569 H/H	569 H/H	569 H/H	569 H/H	569 H/H	rs41830058
Trans Membrane	588..608			574 R/W			novel
		593 A/A	593 A/A	593 A/A	593 A/A		rs68268261
		594 A/A	594 A/A	594 A/A	594 A/A		rs68343169
		605 T/M	605 T/M	605 T/M	605 T/M		rs68343170
					634 A/V		novel
TIR	640..784		644 F/F				novel
					655 V/A
		665 H/Q	665 H/Q	665 H/Q	665 H/Q		rs68268263
		675 H/H	675 H/H	675 H/H	675 H/H		rs68343171
		685 I/I	685 I/I	685 I/I	685 I/I		rs68268264
		738 E/Q	738 E/Q	738 E/Q			rs207552166
		738 E/E	738 E/E	738 E/E	738 E/E		rs68268266
ATG16Lmotif	761..778	765 P/P	765 P/P	765 P/P	765 P/P		rs68268267

,

Table 4. Missense and synonymous variations on protein domain level of TLR4 (according to reference sequence NP_776623.5).

Domain	TLR4 (aa)	AB	EAR	SAR	TG	HOL	SNP ID
LRR1	55..76
LRR2	79..100
LRR3	103..124
LRR4	127..148
LRR5	151..172		151 N/T			151 N/T	rs8193049
LRR6	176..197
LRR7	205..225
			238 N/K				rs8193050
			276 F/F				rs8193051
			347 A/E			347 A/E	rs8193053
LRR8	352..373
LRR9	374..394		374 P/P	374 P/P	374 P/P	374 P/P	rs8193054
				381 K/R			rs8193055
				385 L/L			rs8193056
				389 G/G	389 G/G		rs8193057
LRR10	400..422
LRR11	423..444
LRR12	448..469
LRR13	472..495				482 S/Y		novel
LRR14	497..518		507 Q/Q	507 Q/Q	507 Q/Q	507 Q/Q	rs8193059
LRR15	521..542
LRR16	545..568	552 S/S	552 S/S	552 S/S	552 S/S	552 S/S	rs8193060
LRR17	-
LRR18	-
LRR19	-
LRRCT	579..626			589 S/S			rs8193061
			609 C/C			609 C/C	rs8193062
					622 S/S		rs8193063
			625 N/N			625 N/N	rs8193064
			649 G/G		649 G/G	649 G/G	rs8193065
			640 V/I			640 V/I	rs8193066
Trans membrane	633..653
			664 G/G			664 G/G	rs8193067
		674 T/I	674 T/I	674 T/I	674 T/I	674 T/I	rs8193069
				676 D/D			rs8193070
TIR	677..815

and

Table 5. Missense and synonymous variations on protein domain level of TLR6 (according to reference sequence NP_001001159.1).

Domain	TLR6 (aa)	AB	EAR	SAR	TG	HOL	SNP ID
LRR_RI	<43..164	37 D/N		37 D/N	37 D/N
LRR1	54..77	61 Q/Q	61 Q/Q	61 Q/Q	61 Q/Q		rs68268271
LRR2	78..101	87 R/G	87 R/G	87 R/G	87 R/G		rs68268272
LRR3	102..122	116 S/P	116 S/P
LRR4	123..147	135 D/H	135 D/H	135 D/H	135 D/H		rs520121582
LRR5-LRR6	148..196
LRR7	197..219	214 D/N	214 D/N	214 D/N	214 D/N	214 D/N	rs43702941
		217 A/A	217 A/A	217 A/A	217 A/A		rs68268273
LRR8-LRR12	220.354
LRR13	355..378	374 D/D			374 D/D	374 D/D	rs68268274
LRR14	379..404	395 T/A	395 T/A	395 T/A	395 T/A		rs68268275
		400 K/K	400 K/K	400 K/K	400 K/K		rs211657505
LRR15	405..428	425 S/S	425 S/S	425 S/S	425 S/S		rs68268276
LRR16	429..449
LRR17	450..473		458 H/H	458 H/H			rs68268277
LRR18	474..495
LRR19	496..519	505 N/N	505 N/N	505 N/N	505 N/N		rs55617146
		526 V/A	526 V/A	526 V/A	526 V/A	526 V/A	rs68343174, rs133754378
		526 V	526 V	526 V	526 V	526 V	rs68343175, rs136574510
LRRCT	529..582	539 D/D	539 D/D	539 D/D			rs68343176
			544 V/I	544 V/I			rs55617465, rs68268279
		573 K/K	573 K/K	573 K/K	573 K/K		rs55617193
Trans membrane	585..605	589 V/I	589 V/I	589 V/I	589 V/I	589 V/I	rs55617317, rs207882984
		605 L/L	605 L/L	605 L/L	605 L/L	605 L/L	rs68268280, rs378853146
TIR	641..784	642 F/F	642 F/F	642 F/F	642 F/F		rs438448894
		669 I/V	669 I/V	669 I/V	669 I/V		rs210580164
		674 H/H	674 H/H	674 H/H	674 H/H	674 H/H	rs209572763
		676 R/R	676 R/R	676 R/R	676 R/R		rs68343178
		680 A/A		680 A/A	680 A/A		novel
		684 I/I	684 I/I	684 I/I	684 I/I		rs68343179
		700 F/F	700 F/F	700 F/F	700 F/F		rs55617339, rs211454671
		701 V/V	701 V/V	701 V/V	701 V/V	701 V/V	rs207586910
			709 S/S	709 S/S			rs55617335
				E/E 710			rs68268282

).

Impact of Identified SNPs on Protein Model

Within the identified SNPs, we identified in LRR11 of TLR2 an amino acid change, which resulted in the change from a charged residue into an uncharged residue. This area is part of the TLR2 ligand binding domain, spanning LRRS 9–12. We visualized the importance of this change by constructing a hybrid protein model for one AB individual, carrying the H326Q change, and other variants with one HOL individual. The resulting model was visualized in PyMol [33] (Figure 6).

4. Discussion

The goal of the present study was to assess whether genetic variations in the sequences of innate immune receptors of indigenous cattle breeds compared to those sequences seen for these genes in HF cattle may potentially explain some of the observed genetic resistance to specific bacterial pathogens.

Several association studies have been conducted to find susceptibility related alleles on TLR2, TLR4, and TLR6 genes [34,35,36,37]. None of the previously identified susceptibility associated alleles or haplotypes were found in the present study in indigenous cattle breeds, however, all of them were identified in HOL.

Anatolian breeds seem to be more genetically resistant to infection with bacteria, such as Mycobacterium bovis and mastitis-causing bacteria and, when analyzing TLRs involved in the recognition of these bacterial pathogens, we indeed identified breed specific SNPs within the genes for the receptors investigated which differ significantly from those present in HOL cows.

Interestingly, when grouping the identified SNPs into haplotypes, Anatolian breeds grouped into different haplotypes for TLR2, but not for the other TLRs. This is an important observation, as those bacterial diseases, to which Anatolian breeds have been described to be more resistant, are mainly caused by bacteria binding to TLR2 [38]. Analyzing the identified SNPs in TLR2 in more detail, one aa change was identified in LRR1–10, four between LRR11–LRR20, one in the TM, and four in the TIR domains, respectively, in native cattle breeds, whereas only one was determined in LRR5 in HOL. Considering that the ligand-binding region of TLR2 encompasses LRR9–12, the most important change causing amino acid characteristic changes (H326Q) was found in LRR11. In addition, the aa changes identified in the TIR domain (H665Q and E738Q) need to be further investigated, as these might impact on subsequent intracellular signaling events, similar as described recently for bovine TLR5 [39].

In previous studies, aa changes L227P, H305P, and H326Q in the bovine TLR2 gene have been described to be under positive selection [10]. In our study, the aa at position 227 was L, whereas the aa at position 326 was Q, and no change was determined for the aa at position 305. Furthermore, aa changes L227P, H326Q, N417S, and H665Q have been identified as being specific to Bos indicus cattle breeds [10], similar to T405M [10], which we also identified in our study. It is currently assumed that these variations may represent geographical differences, being driven by a different microbiological environment. It has also been suggested that the Bos indicus-specific aa changes H326Q and R563H might also be found in cattle breeds that originated in a similar geographical and microbial environment [10].

Indeed, blood parasites have been described to cause substantial economic losses in terms of production [40]. Genetic variations in TLR2 of Bos indicus cattle breeds are assumed to impact on blood parasite infections [41]. When the EAR hybrid cattle population breed in Diyarbakir were screened for Theileria blood parasite, only 24 out of 100 samples taken from clinically healthy cows tested positive for Theileria [42]. Similarly, while clinically healthy 24 EAR hybrid cattle (n = 111) and 21 Brown Swiss cows (n = 177) from Erzurum tested positive for the presence of Babesia spp. [43]. Taken further into account that previous phylogenetic studies concluded that Anatolian cattle breeds were Bos indicus and Bos taurus hybrids [44,45], we believe that either similar selective pressure may exists in Anatolian cattle breeds due to the similar geographic/microbial environment compared to pure Bos indicus breeds, or that Bos indicus breeds may have been crossed in at an earlier time due to an increase in the resistance of the local breeds to various infectious diseases.

In comparison with previous studies [34,46,47], the highest number of genetic variation in the analyzed genes were found in Anatolian breeds, except for TLR4, for which we identified more variations in the HOL breed, compared to AB and TG breeds. In these studies, researchers sampled a minimum one individual of each breed that went under artificial selection and analyzed innate immunity-related genes, partially. It is a known fact that the selection pressure for the quantitative traits associated with productions has a negative effect on immunity traits [48,49]. The cited studies showed that cattle breeds were under strong selection pressure, which might lead to a decrease in the variation on the gene regions. Nevertheless, with natural selection animals which cannot resist diseases and environmental conditions cannot find a chance for reproduction; thus, they are eliminated from the population [50] leading to accumulation of resistance-related variations. In the present study Anatolian cattle breeds were indigenous breeds that evolved under natural selection over the years. It can be assumed that the determined polymorphism and haplotype have a potential positive effect on immunity traits. In this context, when taking into account Figure 4 and Figure 5, EAR and SAR, which were raised closer to the center of domestication, were seen as most divergent from the other breeds, suggesting more potential for disease resistance. Determination of the lowest haplotype number in the HOL breed might be due to a low sampling size and/or inbreeding for years. However the TG breed has the same sample size with the HOL breed, and 12 haplotypes were determined for TLR2 (Figure 4).

In addition to TLR2, we also assessed the occurrence of SNPs in TLR4 and TLR6, known to form heterodimers with TLR2. In the TLR6 gene, synonymous SNPs were positively associated with the susceptibility for bovine tuberculosis [51]. In the analyzed individuals of the presented study there were no breed-specific differences for the four synonymous variations. In addition to this, an SNP array study associated protein tyrosine phosphatase receptor T (PTPRT) and myosin IIIB (MYO3B) with bovine tuberculosis resistance [51]. The highest variation was identified in TLR6, but non-synonymous variations were seen in the LRR1–LRR14 in native cattle, whereas only one variation was observed in HOL cattle. One variation found in LRR14 may impact resistance and subsequently might be under positive selection as it was only detected in native cattle breeds [47].

With regards to TLR4, we identified two non-synonymous SNPs in LRR7 (N238K and A347E), which does not contribute to the MD2 binding region. The remaining variants identified all represented synonymous variations. The TLR4 region was analyzed for somatic cell score (SCS) and three SNPs were associated with high and low SCS values (intron1 rs8193046 A/G, exon3 rs8193060 C/T, and 5′UTR rs29017188 C/G) [52]. The ACC haplotype was associated with low SCS, whereas the GTG haplotype was associated with high SCS. High allele frequencies were found for the ACC haplotype in the TG cattle breed. TLR4 is also located in one of the QTL loci for milk production and mastitis [5]. Within the analyzed breeds, the AB breed has the lowest milk production and variation for TLR4, whereas highest variation number was determined in the high milk yield breeds, SAR, EAR and HOL as shown in Figure 5.

5. Conclusions

It is worth mentioning that the number of identified novel InDel variants is significantly higher compared to the identified number of SNPs. We believe that this can be attributed to the fact that InDels were mainly identified in intronic regions, whereas exonic regions in TLRs are more conserved, as well as the fact that the genetic composition Anatolian breeds, representing a hybrid of Bos taurus and Bos indicus breeds, has not been studied before. Given global warming due to climate change and increased anti-microbial resistance due to the overuse of antibiotics, new approaches are needed to manage infectious diseases in farm animals. A key step towards these is the identification of genotypes conferring resistance to both disease and adverse environmental conditions. Our identification of non-synonymous SNPs and novel variants of TLR2, TLR4, and TLR6 genes can provide one such set of variants for future association studies and the validation of the candidate haplotypes at a cellular level.