The following strains were used for all assays as indicated: M. bovis strain AF2122/97 and M. bovis BCG Tokyo; M. tuberculosis H37Rv. Strains were cultured in BD Middlebrook 7H9 media (BD Difco™, USA), supplemented with OADC and 0.05% Tween 80 (Sigma-Aldrich, UK). For M. bovis culture only, growth media was additionally supplemented with sodium pyruvate; for M. tuberculosis culture only, media was additionally supplemented with glycerol (Sigma-Aldrich, UK).

Four lanes of 10x Genomics data were generated using an Illumina HiSeq. This provided ~1200 million reads with a 150bp read length. The data was assembled de novo using the 10x Genomics assembly tool Supernova v2.0.1, using the option to output a phased genome, which generates two homologous assemblies where each scaffold in the assembly has one of the paternal or maternal haplotypes (48). Using one of the haplotypes, we used RepeatMasker, using ‘bovidae’ as the repeat library to soft-mask repeat regions in the assembly (49). We tested the assembly for completeness using several quality control metrics. First, BUSCO (v5.0.0) was used to identify the number of mammalian single-copy orthologs present in the genome (using the database mammalia_odb10) (50). Using the Kmer Analysis Toolkit, we examined the content and distribution of kmers in the assembly to that of the reads using the tool ‘kat comp’, ignoring the first 16 bases of each R1 read in order to omit 10x barcodes from the analyses. We then visualised the distribution of kmers using the tool ‘kat plot’ (51).

The LiftOff tool was used to lift over the B. taurus genome annotation (ARS-UCD1.2.103) to that of the Brown Swiss, noting any genes that were unmapped in Brown Swiss (52). We then identified genes in the Brown Swiss genome annotation that were present in multiple copies (compared to B. taurus ARS-UCD1.2.103) and downloaded the functional information for these genes using ENSEMBL BioMart database (53). For genes annotated as being duplicated three or more times and where a gene name or description was missing, we downloaded the coding sequence for that gene and used megablast to search the NCBI nr database and noted the highest scoring hit (with 100% query coverage and an E-value equal to 0), that contained functional information about the gene (i.e., disregarding hits for clones, BACs, isolates, etc.).

Complementary DNA (cDNA) from TLR2 mRNA sequences isolated from PBMCs from six clinically healthy pedigree BS and four HF were investigated for sequence comparison. All cows were female and housed under Home Office License PPL7009059. Age and lactation cycle of sampled cows were recorded. After identification and selection of candidate SNP H326Q within the coding sequence (CDS) of TLR2 in the BS samples, further cDNA was generated from nine BS and 13 HF PBMC samples from the above-mentioned farms to be assessed for the presence of this SNP. Additionally, 17 Sahiwal and 17 Boran cDNA samples were kindly provided by Drs Thomas Tzelos and Tim Connelley (both The Roslin Institute, University of Edinburgh, UK) for TLR2 SNP analysis.

Blood for peripheral blood mononuclear cells (PBMC) isolation and subsequent macrophage (MØ) generation was collected from clinically healthy pedigree HF and BS cows. All procedures were carried out under the Home Office license (PPL7009059) approved by the RVC’s AWERB Committee. Blood was drawn into sterile glass vacuum bottles containing 10% v/v acid citrate and centrifuged to separate the buffy coat, which was then diluted in phosphate buffered saline (PBS, Thermo Fisher Scientific, UK) with 2% foetal calf serum (FCS, Sigma-Aldrich, UK). Lymphoprep (d = 1.077 g mL^-1, StemCell Technologies ™, Canada) was underlaid and PBMC were isolated by density gradient centrifugation. PBMC were washed and red blood cells lysed, resultant cells were cultured in RPMI 1640 with GlutaMAX™ (Gibco, UK), supplemented with FCS (Sigma-Aldrich, UK) and penicillin/streptomycin (Thermo Fisher Scientific, UK) at 37 °C with 5% CO₂. To derive MØ, PBMC were incubated additionally with 10% filtered L929 supernatant. Media was replaced after three days and adherent cells were harvested after 6 days. Matured MØ were phenotypically assessed by flow cytometry using antibodies specific for bovine cell surface expressed markers including CD14, CD16, CD32, CD11b, CD80, CD163 and MHC class II. All antibodies and isotype controls were purchased from Bio-Rad, UK ( Table 1 ). Data were acquired using a FACS Calibur (BD Biosciences, UK) with Cell Quest Pro software (BD Biosciences, UK), counting 10,000 events. Data was exported as FSC files and analysed using FlowJo V10 (FlowJo LLC, USA).

For sequencing and cloning of the tlr2 gene, total RNA was isolated from bovine PBMC using a RNeasy Mini kit (Qiagen, UK) with on-column gDNA digestion (DNAse from Ambion, UK) according to the manufacturer’s instructions. Total RNA was reverse transcribed to cDNA using iScript cDNA Synthesis Kit (Bio-Rad, UK) according to the manufacturer’s instructions. cDNA samples were checked for concentration and purity via the A260/280 ratio using a Nanodrop ND-1000 Spectrophotometer (Thermo Fisher Scientific, UK). A RT-negative control of each RNA sample was analysed by polymerase chain reaction (PCR) for expression of β-actin to confirm lack of gDNA contamination of RNA.

Full length tlr2 sequences were amplified using primers designed using Primer-BLAST (NCBI, USA) and Primer3Plus ( Table 2 ) (187) based on bovine genome Bos_taurus_UMD_3.1.1 (NCBI, USA) as reference. PCR reactions were performed using Easy-A High-Fidelity PCR master mix (Agilent Technologies, UK) on a Mastercycler Pro Thermal Cycler (Eppendorf, UK). Following amplification, the PCR products were purified using the MinElute PCR purification kit (Qiagen, UK) according to the manufacturer’s protocol. For TA cloning of full length bovine tlr2 into the pDrive plasmid (Qiagen, UK) or pGEM-T plasmid (Promega, UK), manufacturer’s instructions were followed. Plasmids were transformed into NEB^® 10-beta competent E. coli (New England Biolabs, UK) and colony PCR performed to verify the presence of the tlr2 insert. Sequences were validated using Sanger sequencing; DNA sequencing was performed by DNA Sequencing & Services (MRC I PPU, School of Life Sciences, University of Dundee, Scotland, www.dnaseq.co.uk) using the corresponding vector-specific sequencing primers. A minimum of 4 different plasmids derived from unique colonies were sequenced for each animal. Bovine TLR2 contigs were assembled using the bovine TLR2 sequence Bos taurus Hereford breed NM_174197.2 (NCBI, USA) as a reference. TLR2 nucleotide sequences were translated into corresponding amino acid sequences and aligned in CLC Main Workbench V7.6.4 (CLCbio, Denmark). Detected SNPs were classified according to the dbSNP database (NCBI, USA). Analysis of SNP data was curated and as at present, the NCBI dbSNP database contains only human data, allele changes, CDS position and function of each SNP was manually cross-checked and referenced against Ensemble and the European Variation Archive (EVA).

The Human Embryonic Kidney (HEK) 293T cell line, stably transfected with the gene encoding NF-κB inducible secreted embryonic alkaline phosphatase (SEAP), was cultured SEAP media as described recently (54). Bovine or human TLR2 sequences of interest were cloned into pcDNA3.3 TOPO TA mammalian expression vector (Thermo Fisher Scientific, UK), and transfected into SEAP HEK cells using TurboFect transfection reagent (Thermo Fisher Scientific, UK) according to the manufacturer’s protocol. After 24 h, cells were counted and re-seeded in triplicates in a 96-well plate to rest for another 24 h before either evaluation of transfection efficiency or stimulation with various agonists. Surface expression of TLR2 was confirmed by flow cytometry using either human anti-bovine CD282 (anti-TLR2) antibody or mouse anti-human CD282 antibody, both conjugated with Alexa Flour 647 (Bio-Rad, USA). HuCAL Fab-dHLX-MH Ab (Bio-Rad, USA) was used as a negative control as recommended by the manufacturer.

We recently described the SEAP assay and the correlation between SEAP and CXCL8 production in detail (54), see also Supplementary Figure 2 . In brief, SEAP HEK cells expressing TLR2 constructs, or PBMC-derived bovine MØ were specifically stimulated using 100 ng mL^-1 diacylated FSL-1 (a TLR2/TLR6 antagonist; In vivogen, USA); 1 mg mL^-1 of Pam₃CSK₄ (a TLR2/TLR1 antagonist, In vivogen); 19 kDa lipoprotein antigen (represents Rv3763 or Mb3789) (EMC microcollections, Germany), M. bovis BCG (MOI of 10) or MTB H37Rv at (MOI of 5). The direct NF-κB activator phorbol 12-myristate 13-acetate (PMA, Sigma-Aldrich, UK) was used at 200 ng mL^-1 as a positive control. Sterile culture media and mock transfection with an empty vector plasmid were used as negative controls. Supernatants were harvested for CXCL8 ELISA as described by Cronin et al. (55) for bovine MØ, or Quantikine Human CXCL8 ELISA kit (R&D Systems, USA) for human MØ, or colorimetric SEAP assay from HEK cells. SEAP activity as induced by NF-kB activation was assessed by the addition of pre-warmed Quantiblue reagent (In vivogen, USA) followed by incubation for 24 h at 37°C. For both assays, optical density was measured using a SpectraMax M2 spectrometer (Molecular Devices, UK).

We used the 10x Genomics tool Supernova to assemble a haplotype-phased Brown Swiss cow genome. After trimming, the mean read length was 138.5 bases, with a raw read coverage of 59.64. The genome size was calculated to be 2.66Gb (of which 98.08% of the genome was present in scaffolds >10Kb) with a contig N50 of 523Kb and a scaffold N50 of 26Mb ( Supplementary Table 1 ). We used BUSCO to identify reconstructed single copy orthologs. From 9226 BUSCO orthologs 95.2% were recovered, 93.1% were single-copy, 2.1% were duplicated, 1.5% were fragmented, and 3.3% were missing. We used KAT to examine the distribution of kmers in the assembly to that in the reads ( Figure 1 , red distribution). The genome shows a normal distribution of kmers found once in the assembly and once in the reads. It also shows a number of low frequencies kmers (black peak between 0-10) present in the reads, but not in the assembly that most likely represent sequencing errors removed from the assembly. This continues into a shallow distribution of kmers missing from the assembly (black distribution between 10-40), which most likely represents haplotypes found only in the alternate phased assembly.

We noted genes that were missing and duplicated in the Brown Swiss genome, ordering the duplicated genes in order of the most numerous duplications ( Supplementary Table 2 ). We note that all genes annotated as having five or more copies belonged to one of five genes. Three of these genes were repeat/transposon related genes. The remaining two genes were a predicted Bos indicus x Bos taurus elongation factor 1-alpha 1 pseudogene (6 copies), and the Bos taurus T cell receptor gamma cluster 1 (TCRG1) gene (128 copies).

TLR2 is an important PRR involved in the recognition and uptake of MTBC family members. Key SNPs in TLR2 have been identified that affect susceptibility to mycobacterial infection for both humans and cattle (56, 57). The CDS of the bovine tlr2 gene of six BS and four HF cattle were cloned and compared to the bovine NCBI reference sequence genome (RefSeq) comprised of that of a B. taurus Hereford breed bull (NM_174197.2, NCBI, USA, accessed April 2017) in order to identify potential functional SNP sites. A total of 19 SNP variants were detected across the breeds, comprising both synonymous and missense mutations ( Table 3 ). Of the eight missense mutations, seven occur in the ECD region of TLR2 and one in the TIR domain. One synonymous SNP at mRNA position 202 could not be classified with any available database, was described by Jann et al. (47) and is potentially novel. This position was called as ‘A’ in Water Buffalo, ‘G’ in BS and B. taurus, and ‘R’ (purine = G/A) in B. indicus, thus potentially indicating the presence of a heterozygous allele in B. indicus. Of all SNPs identified across the breeds, two of the non-synonymous SNPs have been already reported in the literature (44, 47, 56) and were only observed in the BS animals. One identified SNP variant, rs55617172, leading to an amino acid change of aspartic acid (D) to glutamic acid (E), both negatively charged, at amino acid position 63, has been significantly associated with TB resistance in a TB case- control cattle study (44). However, SNP rs68343167, was only identified in BS isolates in this study, and results in an amino acid change of the positively charged histidine (H) to the uncharged glutamine (Q) at amino acid position 326 (H326Q). This amino acid position has been suggested to be of relevance for ligand binding and functionality of human TLR2 (47, 58), bovine TLR2 and was additionally identified in the Anatolian Black Bos taurus breed (56). In addition, position 326 resides in leucine rich repeat 11 (LRR11) which together with LLR12 are the key domains involved in ligand biding and heterodimerisation of TLR2 with TLR1 and TLR6 (58).

Having identified multiple SNPs within the coding sequence of tlr2 gene across breeds compared to the RefSeq genome, SNP rs68343167 (H326Q) was selected for further investigation using a larger sample size, including cDNA from two additional breeds, Boran and Sahiwal, from the B. indicus species. A total of 66 cDNA samples (Brown Swiss n=15, HF n=17, Sahiwal n=17, Boran n=17) were assessed for H326Q SNP analysis and frequency analysis ( Table 4 ). The reference codon at the site of interest is encoded by the CAT sequence (TT genotype). Individual animals carrying the SNP had either heterozygote codons (TA genotype), encoded by either CAT or CAA, which is responsible for the H326Q variant; or they were homozygote for the CAA sequence (AA genotype).

Interestingly, all HF samples tested as part of the present study were homozygous for the B. taurus Hereford RefSeq (TT genotype). Of the BS samples, six samples were homozygous for the RefSeq sequence, whereas the majority (n=9) were found to be heterozygous (TA genotype), containing both, the RefSeq sequence and the H326Q SNP. Of the Bos indicus breeds, the Boran samples revealed two samples homozygous for the RefSeq sequence, eight heterozygous and seven homozygous for the SNP. Thus, the heterozygote genotype was most abundant in the BS and Boran population. All Sahiwal samples (n=17) were homozygous for the SNP variant (AA genotype).

Alignment of TLR2 transcripts revealed additional SNP sites unique to BS, BS and B. taurus as well as B. indicus. Of note are positions H326Q as reported here and by others, as well as positions 417, 563 and 665; these latter 3 amino acid positions are particularly interesting as like H326Q these SNPs were found to introduce identical residue changes to our full length TLR2 analysis. Further, these SNPs were found to be present in BS cattle investigated ( Table 3 ) but not present for HF cattle suggesting a role for genotype impact as found for H326Q. Position 417 and 665 are noteworthy, belonging to LRR15 and the TIR domain respectively. As LRR15 is close to the LRR regions thought to be involved in dimerization of TLR2 with co-receptors this SNP may play a role in effective PRR dimers forming. As the TIR domain is crucial for receptor dimerization of TLR cytoplasmic domains as well as signal potentiation; SNP H665Q may also modulate functional responses of TLR2.

On establishing distinct genotypes for the H326Q SNP in a range of cattle species (i.e., HF = TT, BS and Boran = TA and Sahiwal = AA), we assessed in the next step this association with TLR2 gene models and global sequence alignments. From the genotype analysis BS and Boran cattle represent an intermediate between the genotype extremes, perhaps representing genomic heritage; for example, Boran cattle are predominantly B. indicus with underlying influences of European and African B. taurus (59, 60). TLR2 gene models comparing the BS genome to B. taurus, B. indicus and B. bubalis (water buffalo) highlights an alignment of BS with B. taurus as TLR2 is derived from a single transcript in both breeds ( Supplementary Figure S3 ). Highlighting the heterozygous genotype of BS cattle, BS resembled B. taurus at H326Q within LRR11 of the extracellular domain of TLR2 ( Supplementary Figure S4 ) when the CDS of each possible TLR2 transcript was aligned to the three cattle species, Interestingly, additional sites where the BS and B. taurus sequences differed from B. indicus where I211V, F227L and R337K. Positions 211 and 337 were identified from the initial TLR2 sequencing described above.

Having identified the H326Q candidate SNP and its representation in the present study population of different cattle breeds, specifically in the BS population, the potential impact of this SNP on boTLR2 function in this breed was assessed. We utilised a previously published HEK cell assay to express the different TLR2 constructs in the absence of additional PRR (54). As well as the BS H326Q TLR2 variant, we included the Hereford reference variant and a L306P SNP variant that has been identified by others (47, 56) but not identified in the present study. We were also interested to assess the functionality of TLR2 across species, and therefore included the huTLR2 reference sequence (NCBI RefSeq Accession Number NG_016229.1). Furthermore, novel chimeric TLR2 constructs featuring the ECD of huTLR2, with the TIR domain of boTLR2 were created to delineate the contribution of each component to TLR2 functionality.

To confirm that the different TLR2 constructs were similarly expressed on the cell surface of the SEAP HEK cells, immunostaining was performed using bovine and human TLR2 specific antibodies. The mean fluorescence intensity (MFI) values of the positively stained population were similar for all constructs ( Figure S1 ). Stimulation with PMA, a non-specific activator of NF-κB, was used as a positive control and resulted in consistently high SEAP secretion for all transfectants. SEAP activity in response to ligands was normalised to PMA responses as an internal control.

The H326Q variant showed significantly higher TLR2-dependent NF-кB response in comparison to RefSeq Hereford TLR2 (p < 0.001 Figure 2 ). Constructs containing the L306P SNP did not respond significantly differently in comparison to the RefSeq Hereford boTLR2 version. The response of SEAP HEK cells transfected with huTLR2 was significantly higher in response to both FSL-1 and Pam₃CSK₄ in comparison to the bovine constructs. Interestingly, the chimeric bov:hu TLR2 receptor showed higher NF-κB activity than the RefSeq Hereford TLR2 HEK cell line, but significantly lower activity than the huTLR2 receptor in response to both FSL-1 and Pam₃CSK₄.

Using cell culture supernatant recovered from the above SEAP assay experiments, CXCL8 production in response to FSL-1 and Pam₃CSK₄ stimulation were measured ( Figures 3B, D ). As previously, there was an increase in CXCL8 production in the H326Q variant in comparison to the bovine RefSeq Hereford TLR2, however this was found only to be statistically significant in response to FSL-1 stimulation (p<0.01). CXCL8 production measured in huTLR2 transfectants was significantly higher than in boTLR2 transfectants in response to both ligands. CXCL8 production by chimeric bo:hu TLR2 HEK cells was significantly higher than any of the boTLR2 variants to FSL-1 (p<0.05), but lower than those from huTLR2, reflecting the observations made using the SEAP reporter assay of NF-κB activity.

Having assessed the activity of TLR2 variant HEK cells using TLR2-specific agonists and determining good correlation between SEAP activity and CXCL8 secretion ( Supplementary Figure S2 ), we next challenged the TLR2-HEK cells with live mycobacteria and measured the NF-κB activation responses. Interestingly, the response of SEAP HEK cells expressing huTLR2 to challenge with M. tuberculosis was markedly lower than when they were challenged with TLR2-specific lipopeptides i.e., Pam₃CSK₄ and FSL-1 in the previous experiment ( Figure 3 ). Only the huTLR2 variant produced a significant NF-κB response upon challenge with M. bovis BCG compared to mock simulations. NF-κB responses for H326Q containing TLR2 were significantly lower than Hereford reference TLR2 and those constructs containing L306P. All boTLR2 variants responded more strongly to M. tuberculosis than to M. bovis BCG, even though the MOI was lower (MOI of 5) than with the M. bovis BCG infection (MOI of 10).

Differences in functionality at the whole cell level conferred by boTLR2 breed variants were further assessed using PBMC-derived MØ. Cells were isolated from BS, containing the heterozygote H326Q variant, or HF cattle, which have the TLR2 sequence identical to the RefSeq Hereford breed. Due to Cat3 access restrictions, we measured the CXCL8 response to the TLR2-specfic agonist, FSL-1 and M.bovis BCG. In agreement with TLR2-HEK cell responses ( Figure 4 ), MØ generated from BS carrying the H326Q variant showed a significantly higher CXCL8 responses than HF MØ upon FSL-1 stimulation (p<0.05, Figure 5A ) and the 19 kDa protein, but not M.bovis BCG ( Figure 5B ). We further verified a high degree of correlation between the SEAP NF-kB and CXCL8 assays in a control experiment, by collecting data for both from the same stimulated cells and performing regression analysis ( Figure S2 ). Interestingly, unlike the SEAP HEK cell experiments, the PBMC-derived MØ responded similarly to MTB-derived 19 kDa lipoprotein and M. bovis BCG compared with FSL-1. We also detected a significantly greater response by BS derived MØ to M. bovis BCG compared to untreated cells, while the response by HF-derived MØ was not statistically greater than untreated cells.