Confirming the identity of GX buffaloes
The breed of GX buffaloes was confirmed as a swamp-type based on a combination of morphological and karyotyping analysis, and sequencing of the mitochondria D-loop region. The morphological features of GX buffaloes(swamp-type) were compared with Murrah and NR buffaloes (river-type) as shown in Additional file 1: Figure S1a. The native GX buffaloes had phenotypic characteristics typical to the swamp-type; with a flat forehead, wide oronasal area, large abdominal perimeter, long sides, and backwardly swinging horns. The growth direction of the horns is horizontal to the forehead. The spine of GX buffaloes suddenly becomes flat after the last thoracic region. The Murrah and NR buffaloes exhibited typical river-type features; with a well-developed udder, and taller and wider body, wider pelvis, upwardly-curled horns, deeply pigmented skin and hair, and spine that extends further into the hindquarters (Cockrill 1976; Zhang et al. 2000; Presicce 2017).
Karyotyping analysis showed that GX buffaloes possess 2n = 48 chromosomes, of which 10 are meta- and submetacentric, 36 acrocentric somatic chromosomes, and 2 acrocentric sex chromosomes (see Additional file 1: Figure S1b and c) - features that match the karyotype characteristic of swamp-type buffaloes (Degrandi 2010). The nucleotide sequence of the mitochondria D-loop region of GX buffaloes was amplified and sequenced (NCBI accession number: MG491488). NJ tree showed a separation between swamp-type buffaloes and river-type buffaloes, but GX buffaloes were more closely related to swamp-buffaloes, further supporting the breed identity of GX buffaloes as a swamp-type. The similarity of D-loop nucleotide sequences was highest between GX, and WZ and XY buffaloes (99.8%). High similarity was also found with DC (99.7%), XL (99.7%) and FA (95.2%). There was less similarity with river-type buffaloes: Mediterranean (93.3%), Jaffarabadi (92.8%), Surti (92.5%), NR (91.5%), and Murrah (89.3%); similarity with the outgroup species of the genus Bison and genus Bos was 80.1–83.2% (see Additional file 1: Figure S1d). These results show that GX buffaloes are of a swamp-type, and are closely related to the Southwest Chinese cluster (DC, XL, WZ, XY)(Lei et al. 2007; Zhang et al. 2007).
TLRs 1–10 sequence analysis
We amplified and sequenced full-length TLRs 1–10 of GX buffaloes. The primer sequences, annealing temperature and expected PCR product length are shown in Table 1. The electrophoresis result shown in Fig. 1 indicates that all DNA fragments were amplified accorded to the expected size of each amplicon. The full-length sequences of TLRs 1–10 genes of GX buffaloes have been submitted to the NCBI GenBank, and the accession numbers are listed in Table 1. The coding sequences and deduced AA of these TLRs were compared with those of NR buffaloes, cattle, American bison, zebu, yak, wild yak, gayal, goat and sheep as “ingroups” (Table 3); then with other non-ruminant species including horse, pig, mouse, rat, rabbit, monkey and human as “outgroups” (see Additional file 2: Table S1).
Table 1
Primer sets, PCR annealing temperature and amplicon size, and gene accession number for swamp buffalo TLRs 1–10 genes
TLR gene
|
Primer sequence (5’-3’)a
|
Amplicon size (bp)
|
Annealing temperature (°C)
|
Accession number b
|
TLR1
|
F: AAAAACTTATCCCTGGAAACA
R: TATAAATCCCCGTGTTAATGT
|
2261
|
50
|
MG491497
|
TLR2
|
F: ATGCCACGTGCTTTGTGGA
R: TTATTGCAGCTCTCAAATTTAACCA
|
2347
|
50
|
MG491498
|
TLR3
|
F: AAAACGAACTGGATTTGGACTAA
R: TGGGTCAGATTATGTCAATTGGTA
|
2982
|
50
|
MG491489
|
TLR4
|
F: TCCTAAATTCACTAACCCTGAGGT
R:CTTAGCTCTAATTCAGTTTGATGGG
|
2986
|
52
|
MG491490
|
TLR5
|
F: CCTCAAGGCTGTTGCTCAGT
R: GCCTAGGAGATGGTGGTTACA
|
2651
|
52
|
MG491491
|
TLR6
|
F: AGAAGCTGTCATTTTGTTTACAT
R: AACCTTTCACATCATCATTTTC
|
2361
|
50
|
MG491492
|
TLR7
|
F: GGGTCCATCTCAGGCTGATC
R: TCCACTTCCTTTGTCATGTATTG
|
3393
|
55
|
MG491493
|
TLR8
|
F: GGTCAAAGCACTGCGAATTTCA
R: TTGCTGTGCATTTTGACTATGC
|
3238
|
50
|
MG491494
|
TLR9
|
F: CTAAGGGCTGCTCTCAGTGGTGGT
R: ACTCTGTGCTATTCGGCTGTCGTG
|
3845
|
52
|
MG491495
|
TLR10
|
F: AAAATTCATCTGAGGTGAACCAGTG
R: TAGGTCCCAATCTCCCTAGGAGTG
|
2517
|
55
|
MG491496
|
a Primer orientation: F, forward; R, reverse. |
b GenBank accession numbers for full-length sequences of TLRs obtained in the present study. |
Table 3
The similarity of the coding sequence (CDS), deduced amino acid (AA) sequence and Toll/IL-1 receptor (TIR) domain sequence of Toll-like receptor (TLR) genes between GX buffaloes and nine other ruminant species
Name
|
Species
|
CDS identity (%)
|
AA identity (%)
|
TIR identity (%)
|
GenBank No.
|
TLR 1
|
river buffalo
|
96.9
|
99
|
99.3
|
HQ327989.1
|
|
cattle
|
96
|
97.8
|
100
|
NM_001046504.1
|
|
wild yak
|
94.7
|
97.7
|
100
|
XM_005897339.2
|
|
zebu
|
94.5
|
97.2
|
99.3
|
XM_019962587.1
|
|
yak
|
97.1
|
97.2
|
100
|
KJ101605.1
|
|
gayal
|
95.7
|
96.1
|
97.9
|
KF905597.1
|
|
bison
|
94.9
|
97.8
|
100
|
XM_010833592.1
|
|
goat
|
96.1
|
94.6
|
97.2
|
HQ263209.1
|
|
sheep
|
95.3
|
96.6
|
99.3
|
NM_001135060.2
|
TLR 2
|
river buffalo
|
96.6
|
98.9
|
98.6
|
HM756162.1
|
|
cattle
|
92.7
|
97.1
|
97.9
|
NM_174197.2
|
|
wild yak
|
90.4
|
95.7
|
92.5
|
XM_014480820.1
|
|
zebu
|
91.5
|
96.8
|
97.2
|
XM_019978018.1
|
|
yak
|
96.8
|
96.9
|
97.9
|
KF977429.1
|
|
gayal
|
96.3
|
96.3
|
96.6
|
KF905578.1
|
|
bison
|
97.5
|
96.7
|
96.6
|
EU580539.1
|
|
goat
|
97.9
|
96.9
|
95.2
|
NM_001285603.1
|
|
sheep
|
95.2
|
92.3
|
91.7
|
NM_001048231.1
|
TLR 3
|
river buffalo
|
96.7
|
95.8
|
98.6
|
HQ343417.1
|
|
cattle
|
96.9
|
97
|
97.9
|
NM_001008664.1
|
|
wild yak
|
95.2
|
97.8
|
99.3
|
XM_005899004.2
|
|
zebu
|
95.1
|
97.6
|
99.3
|
XM_019953354.1
|
|
yak
|
98.2
|
97.8
|
99.3
|
KF990165.1
|
|
gayal
|
93.5
|
97.9
|
98.7
|
KF905599.1
|
|
bison
|
95.2
|
97.9
|
99.3
|
XM_010835803.1
|
|
goat
|
93.9
|
89.6
|
94.9
|
HQ263210.1
|
|
sheep
|
96.8
|
96.2
|
99.3
|
NM_001135928.1
|
TLR 4
|
river buffalo
|
93
|
97
|
97.9
|
HQ343416.1
|
|
cattle
|
91.8
|
96.7
|
97.9
|
NM_174198.6
|
|
wild yak
|
91.6
|
96.3
|
97.9
|
XM_005891938.1
|
|
zebu
|
95
|
96.2
|
98.6
|
KX138607.1
|
|
yak
|
97.4
|
96.3
|
97.9
|
KF977430.1
|
|
gayal
|
85.4
|
96.4
|
98.6
|
HQ327993.1
|
|
bison
|
86.8
|
96.6
|
98.6
|
AH015104.2
|
|
goat
|
93.6
|
93.2
|
97.3
|
NM_001285574.1
|
|
sheep
|
95.9
|
93.7
|
97.3
|
NM_001135930.1
|
TLR 5
|
river buffalo
|
98.2
|
98.7
|
100
|
NM_001290918.1
|
|
cattle
|
94.9
|
94.9
|
95.2
|
NM_001040501.1
|
|
wild yak
|
94.8
|
97.9
|
97.9
|
XM_005903731.1
|
|
zebu
|
97.9
|
97.6
|
97.3
|
GQ866979.1
|
|
yak
|
95.7
|
97.3
|
95.9
|
GU647093.1
|
|
gayal
|
92.5
|
92.4
|
91.4
|
HQ392514.1
|
|
bison
|
95.2
|
97.7
|
97.9
|
JN615236.1
|
|
goat
|
96.3
|
92.4
|
95.7
|
NM_001285699.1
|
|
sheep
|
96.8
|
95.1
|
97.9
|
NM_001135926.1
|
TLR 6
|
river buffalo
|
96.6
|
94.7
|
96.1
|
HQ327992.1
|
|
cattle
|
97.7
|
96.9
|
96.1
|
NM_001001159.1
|
|
wild yak
|
98.3
|
97.8
|
99.2
|
XM_005897340.1
|
|
zebu
|
98.1
|
97.7
|
98.4
|
XM_019962585.1
|
|
yak
|
98.4
|
98
|
99.2
|
KM359137.1
|
|
gayal
|
98.7
|
98
|
N/A
|
KC161554.1
|
|
bison
|
98.4
|
98.1
|
98.4
|
XM_010833587.1
|
|
goat
|
95.1
|
91
|
90.7
|
NM_001285540.1
|
|
sheep
|
96.6
|
94.3
|
99.2
|
NM_001135927.1
|
TLR 7
|
river buffalo
|
98.5
|
99.9
|
100
|
GU903503.1
|
|
cattle
|
98.5
|
99
|
100
|
NM_001033761.1
|
|
wild yak
|
97.4
|
99
|
100
|
XM_005891026.2
|
|
zebu
|
97.4
|
99
|
98
|
XM_019956405.1
|
|
yak
|
94.8
|
97.4
|
100
|
GQ438259.1
|
|
gayal
|
95.3
|
98.2
|
N/A
|
KF905583.1
|
|
bison
|
91.9
|
97.1
|
100
|
XM_010830634.1
|
|
goat
|
94.7
|
92.4
|
86.1
|
NM_001285543.1
|
|
sheep
|
96.9
|
98.5
|
100
|
NM_001135059.1
|
TLR 8
|
river buffalo
|
99.7
|
99.4
|
99.3
|
NM_001290928.1
|
|
cattle
|
97.2
|
96.6
|
99.3
|
NM_001033937.1
|
|
wild yak
|
97.1
|
98.1
|
100
|
XM_005891027.2
|
|
zebu
|
94.3
|
98
|
100
|
XM_019956418.1
|
|
yak
|
98.4
|
97.9
|
100
|
KM359140.1
|
|
gayal
|
98.1
|
97.3
|
98
|
KF905575.1
|
|
bison
|
96.3
|
98.2
|
100
|
XM_010830635.1
|
|
goat
|
99.3
|
98.8
|
99.3
|
NM_001285757.1
|
|
sheep
|
97.7
|
97.6
|
99.3
|
NM_001135929.1
|
TLR 9
|
river buffalo
|
96
|
97.2
|
100
|
HQ242779.1
|
|
cattle
|
94.7
|
98.2
|
98.6
|
NM_183081.1
|
|
wild yak
|
97.1
|
96.5
|
99.2
|
XM_014478899.1
|
|
zebu
|
96.8
|
98.3
|
99.3
|
KX138608.1
|
|
yak
|
97.8
|
97.9
|
98
|
KM359141.1
|
|
gayal
|
95.9
|
96.3
|
98.6
|
KF905585.1
|
|
bison
|
83.5
|
98.4
|
99.3
|
XM_010851751.1
|
|
goat
|
95.7
|
95.2
|
98
|
EU747825.1
|
|
sheep
|
94.2
|
95.3
|
98
|
NM_001011555.1
|
TLR 10
|
river buffalo
|
94.3
|
95.2
|
95.9
|
HQ327991.1
|
|
cattle
|
94.8
|
96.9
|
90.4
|
NM_001076918.2
|
|
wild yak
|
94.9
|
96.9
|
98.6
|
XM_005897338.1
|
|
zebu
|
95.8
|
96.7
|
97.9
|
EU006659.1
|
|
yak
|
95.7
|
96.1
|
98.6
|
KF776527.1
|
|
gayal
|
96
|
96.7
|
98.6
|
KF905606.1
|
|
bison
|
94.7
|
97
|
98.6
|
XM_010833593.1
|
|
goat
|
93.1
|
87.8
|
90.4
|
NM_001285541.1
|
|
sheep
|
96.6
|
95.4
|
95.9
|
NM_001135925.1
|
The coding sequence, deduced AA sequence and the predicted TIR domain of most TLRs were conserved across the nine ruminant species (especially river-type buffaloes), but were divergent from non-ruminant mammals. There was a high degree of similarity (> 90%) with coding sequences of ruminant (river-type buffalo, cattle, sheep and goat) TLRs, whereas a 67–89% similarity was detected to the sequences of horse, pig, mouse, rat, rabbit, monkey, and human. At the AA level, the similarity of TLRs between GX buffaloes and river-type buffaloes was 94.7–99.9%, whereas a89.6-99% similarity with ruminant species and 48–91% similarity to non-ruminant species were detected. TLR7 of GX buffaloes showed the highest amino acid homology to that of river-type buffaloes (99.9%). Although the coding sequences of TLR 6 of GX buffaloes had 96.6% identity with river-type buffaloes, TLR6 AA sequence of GX buffaloes was more similar to yak, gayal or bison (~ 98%), cattle (96.9%), wild yak (97.8%) and zebu (97.7%); than to NR buffaloes (94.7%). AA sequences of GX buffalo TLR9 and TLR10 were more similar to domestic cattle, zebu or bison (~ 98% and 97%, respectively) compared to NR buffaloes (97.2% and 95.2%, respectively). The lower AA sequence similarities in TLRs 6, 9 and 10 between the swamp-type and river-type buffaloes point out to breed-specific differences in the biological function of these TLRs. Compared with non-ruminant mammals, the similarities between TLRs of GX buffaloes and ruminant species were high, suggesting that TLRs might have evolutionary conserved functions across ruminants.
Prediction of functional domains
A simple modular architecture tool (SMART) was used to predict the functional domains of TLRs, including the signal peptide (SP), transmembrane(TM), TIR and LRR domains (Fig. 2). Similar to river-type buffaloes, TLRs 1–10 of swamp-type buffaloes possessed LRR and TIR domains, however the TM domains varied across different TLRs. No TM domain was found in TLR7 nor TLR9, of both GX buffaloes and NR buffaloes. This was consistent with an earlier report on river-type buffaloes (Dubey et al. 2013). Surprisingly, only one TM domain was found at the C-terminus of TLR5 and TLR6in GX swamp-type buffaloes, rather than the two TM domains found in the homologous TLRs of river-type buffaloes, which was inconsistent with a previous study(Dubey et al. 2013). TMHMM was employed to confirm the finding, and the result completely matched the SMART prediction. Instead of the N-terminus TM domain, a SP domain was found. This feature suggests that TLR5 and TLR6 of swamp-type buffaloes are secretory proteins.
The cytoplasmic TIR domain plays a central role in cell signal transduction cascade downstream of TLRs through adaptors, such as myeloid differentiation factor 88 (MyD88), TIR domain-containing adaptor inducing interferon (TRIF), TRIF-related adaptor molecule(TRAM) and TIR domain containing adaptor protein(TIRAP)/Mal(Yamamoto et al. 2004; Kawai and Akira 2010). Sequences of the intracellular signal-transducing TIR domains were subjected to multiple alignments across species. TIR domains were found to be conserved across ruminant TLRs, with a similarity that ranged from 86.1–100%. TIR identity of TLR5, 7 and 9 between swamp-type and river-type buffaloes was 100%, whereas that of TLR6 (96.1%) and TLR10 (95.9%) was relatively lower, when compared to other TLRs. For TLR1 and TLR8, the TIR similarities were up to 100% between GX buffaloes and Bos species, but were slightly lower (99.3%)between GX buffaloes and NR buffaloes. Although GX buffaloes are more closely related to NR buffaloes than other ruminant species in the phylogenetic tree, the result indicated that the TIR domain of TLR6 was more conserved among GX buffaloes, wild yak, yak and sheep (99.2%) compared to river-type buffaloes (96.1%). There was also a 98.6% similarity in TIR10 between GX buffaloes and wild yak, yak, gayal or bison; but only a 95.9% similarity between GX buffaloes and NR buffaloes, and a 90.4% similarity between GX buffaloes and domestic cattle. These results indicate that coding sequences, deduced AA sequence and predicted TIR domains of most TLRs, were relatively more conserved among ruminants - but were divergent from human, mouse and chicken, as expected.
The highly-conserved LRR regions, located extracellularly, are associated with pathogen ligand binding(Kobe and Kajava 2001). To predict the module number of LRR motifs present in TLRs of GX buffaloes, the LRR finder web application based on a large LRR database (tLRR db) was used (Matsushima et al. 2007; Offord et al. 2010). The number of LRRs in TLRs 1–10 across various mammal species was compared, and the results are presented in Table 4. In swamp-type buffaloes, the number of LRR modules varied from 18 (TLR10) to 26 (TLRs 7, 8 and 9), which was similar to the range reported in river-type buffaloes(Dubey et al. 2013). Most TLRs of swamp-type buffaloes shared an equal number of LRR modules with river-type, except for TLRs 1, 6 and 10. TLRs 7, 8 and 9 were found to have an equal number of LRRs among sequences analyzed from all ruminant species (26 LRRs), except gayal TLR7 (25 LRRs), bison TLR7 (uncertain) and cattle TLR8 (25 LRRs), suggesting their conserved nature among ruminant species. GX buffalo TLRs 1, 2, 3, 4, and 5 had an equal or similar (+/- one LRR) number of LRR modules to cattle, wild yak, zebu, yak, bison and sheep. TLR6 and TLR10 of swamp-type buffaloes and zebu had equal number of LRRs, whereas the number of LRRs in these two TLRs in swamp-type buffaloes was similar to other ruminant species. The mammalian LRR motif has been shown to play a role in protein-protein interaction(Kobe and Deisenhofer 1994). Duplication or deletion of individual LRR motifs may alter the protein’s stability and affinity, and subsequently change the function of ligand recognition in some TLRs (Björklund et al. 2006; Peter et al. 2009; Vieux and Barrick 2011). Our result indicates a possible evolutionary divergence in most of the TLR1 subfamily (TLRs 1, 6 and 10) between swamp-type and river-type buffaloes, which might lead to differential recognition of TLR ligands, despite their highly homologous sequences.
Table 4
Number of Leucine-rich repeats (LRRs) identified based on the amino acid sequences of the 10 Toll-like receptor (TLR) genes in GX buffaloes and nine other ruminant species. Accession numbers of TLRs used for the prediction analysis are provided in Table 3
TLRs
|
GX buffalo (this study)
|
NR buffalo
|
Cattle
|
Wild yak
|
Zebu
|
Yak
|
Gayal
|
Bison
|
Goat
|
Sheep
|
TLR1
|
19
|
18
|
19
|
19
|
19
|
20 a
|
17
|
20
|
19
|
19
|
TLR2
|
20
|
20
|
20
|
20
|
19
|
20
|
17
|
19
|
20
|
20
|
TLR3
|
24
|
24
|
24
|
24
|
24
|
24
|
17a
|
24
|
21
|
24
|
TLR4
|
22
|
22
|
22
|
22
|
21
|
21
|
22
|
22
|
22
|
22
|
TLR5
|
23
|
23
|
23
|
23
|
23
|
23
|
22
|
23
|
23
|
23
|
TLR6
|
19
|
20
|
20
|
20
|
19
|
20
|
12 a
|
20
|
18
|
20
|
TLR7
|
26
|
26
|
26
|
26
|
26
|
26
|
25
|
14 a
|
26
|
26
|
TLR8
|
26
|
26
|
25
|
26
|
26
|
26
|
26
|
26 a
|
26
|
26
|
TLR9
|
26
|
26
|
26
|
26
|
26
|
26
|
26
|
26
|
26
|
26
|
TLR10
|
18
|
20
|
20
|
20
|
18
|
19
|
20
|
20
|
20
|
20
|
a denotes partial AA sequence. |
SMART analysis showed that TLRs 2, 4, 5, 6, 7 and 8 of GX buffaloes contained a SP domain at their N-terminus, suggesting their role as secretory proteins, however, TLRs 1, 3, 9 and 10 were non-secretory. As for NR buffaloes, the SP domain was absent in TLRs 1, 3, 5, 6 and 10, but present in the rest of the TLRs. The absence of a SP domain in TLRs 1, 6 and 10 had been noted in other Indian river-type buffalo breeds (Banerjee et al. 2012). It is clear that the occurrence of SP domains in TLRs can vary between breeds during their evolutionary divergence. The result of the SMART prediction concurred with the result of the SignalP4.1 prediction(Nielsen 2017)(see Additional file 3: Figure S2). SignalP-4.1 also provided the cleavage sites of GX buffalo TLRs 2, 4, 5, 6, 7 and 8, which were located at 20–21 aa, 25–26 aa, 20–21 aa, 22–23 aa, 26–27 aa and 18–19 aa, respectively. SPs are critical for protein translocation and export(von Heijne 1990). A TLR carrying a SP domain determines the N-terminus of the mature TM protein on either the non-cytoplasmic side or cytoplasmic side (Käll et al. 2004). Combined with the results from the TM prediction, it is possible to suggest that TLRs 2, 4, 5, and 6 of GX buffaloes are located on the plasma membrane. The only exception was the presence of a SP domain in TLR3 of both buffalo types based on our SignalP-4.1 analysis, but this was not revealed by SMART prediction. This result might be attributed to differences in the algorithms used in the two approaches; however, the presence or absence of a signal peptide does not completely control the localization of a protein.
The LRRCT and LRRNT domains (cysteine flanking cap regions of LRR at the C- and N-terminus, respectively) were also predicted in the GX buffalo TLRs 1–10. LRRCT was found to be present in all GX swamp-type buffalo TLRs except TLR9, as was predicted in river-type buffalo TLRs; whereas LRRNT was absent in all TLRs, except TLR3 and TLR7, as expected (Dubey et al. 2013). The high similarity in the structural features of LRRCT and LRRNT, suggest a similar function in receptor activation, or ligand-binding events, among TLRs 1–10 in these two buffalo breeds.
Prediction of protein structure
Analysis of the predicted GX buffalo TLRs 1–10 amino acid sequences indicated that the putative AAs encoded by ORFs were 727 aa, 877 aa, 904 aa, 891 aa, 871 aa, 736 aa, 1050 aa, 1033 aa, 1081 aa and 812 aa, respectively. The results of the prediction of the physical and chemical properties of these proteins can be found in Additional file 4: Table S2. The GX buffalo TLR 2 and 9 were positively charged, according to the predicted charge characteristics and theoretical pI; whereas the other TLRs were negatively charged. A comparative analysis of the AA composition of the putative product encoded by the GX buffalo TLRs 1–10 genes, showed that Leu had the highest frequency (14.4%-18.2%), followed by Ser (8%-10.5%); both are very similar to the content of Leu and Ser in NR buffalo TLRs 1–10 (14.4%-18.7% and 8%-10.5%, respectively) (see Additional file 5: Table S3). Analysis of the hydrophobicity of TLR proteins, in this study, was based on the machine-calculating grand average of hydropathicity (GRAVY) values. The GRAVY values of GX buffalo TLRs 1–10 were − 0.046, -0.203, -0.125, 0.031, 0.026, -0.057, -0.13, -0.135, -0.035 and − 0.053, respectively. These values revealed that TLR4 and TLR5 were hydrophobic proteins, and the rest of TLRs were all hydrophilic proteins. The hydrophobicity of TLRs 1–10 in NR buffaloes showed consistent results with GX buffaloes (the comparative results are shown in Additional file 6: Figure S3).
Combined analysis using the ProtComp, TargetP and PredictProtein programs, indicated that TLRs 1, 2, 6, 10 are more likely to be located in both the cytoplasm and membrane, whereas TLRs 4, 5, 6 are located on the plasma membrane (extracellular), and TLRs 3, 7, 8, 9 are in the cytoplasm or intracellular organelles. Our analysis also showed a high possibility for the existence of a secretory pathway - characterized by containing a SP in TLR3 and TLR4. Prediction of transmembrane helices in TLR proteins was performed by TMHMM 2.0 (see Additional file 7: Figure S4); and showed that the classical type I TM domain was present in TLRs 1, 2, 3, 4, 5, 6, 8 and 10, whereas TM did not appear in either TLR7 and TLR9 of GX buffaloes, which was similar to NR river-type buffaloes, and agreed with the prediction result obtained by SMART analysis. The biological relevance of the absence of a TM domain in bubaline TLR7 compared to other ruminant species remains to be elucidated.
Additionally, to further characterize the tertiary structure, protein homology models for these GX buffalo TLRs were generated by applying the SWISS-MODEL Homology Modelling program (see Additional file 8: Figure S5). The 3D graphics may provide information for further predictions on the active domain of proteins to help design appropriate ligands or reagents. The QMEAN4 ranged between − 2.22 and − 7.28, while GMQE score ranged between 0.55 and 0.77 (if excluding TLR10 which found unmatched template of TLR10 from the database), indicating an overall good quality for most models (see Additional file 9: Table S4).Although this is a rough analysis of homologues for 3D structure for these bubaline TLRs, it has however been revealed that all 3D models for these TLR proteins were closely matched to the homologous TLR templates (with the exception of TLR10), what preliminarily suggests their conservation of tertiary structure therefore the possible biology to those of template TLRs that have been commonly studied.
Evolution of ruminant TLR genes, TIR and LRR domains
The phylogenetic tree revealed the relative genetic distance between different ruminant species, especially between the two types of buffaloes examined in our study. The AA sequences of TLRs 1–10 from GX buffaloes, NR buffaloes and the examined ruminant species, were subjected to phylogenetic analysis using the NJ algorithm with 1,000 bootstrap resampling (Fig. 3a). The same analysis was used for TIRs and LRRs (Fig. 3b and c). It is clear that clustering of major TLR gene subfamilies with high bootstrap values, could be divided into three clades: TLR1 subfamily (TLRs 1, 2, 6 and 10), TLR7 subfamily (TLRs 7, 8 and 9), and the singletons (TLRs 3, 4 and 5) together as a single clade originating from the same branch.
Regarding the singleton TLRs (TLRs 3, 4 and 5), GX buffalo TLR4 and TLR5, but not TLR3, were highly conserved between the two buffalo breeds - compared with species of genus Bos or Bison. There was a higher similarity between TLR2 or TLR8 of goats with GX buffaloes, but not with river-type buffaloes or any Bos species. An additional phylogenetic analysis of TLR AA sequences - comparing the evolutionary relationship of GX buffaloes, major ruminant species and non-ruminant mammals -revealed a closer relationship of GX buffalo TLRs to those of ruminants than to other mammalian TLRs (see Additional file 10: Figure S6), a finding consistent with previous studies(Banerjee et al. 2012; Dubey et al. 2013).
mRNA expression in selected tissues from GX buffaloes
Transcripts of TLRs 1–10 genes and the reference gene GAPDH were successfully quantified in all examined tissues and blood leukocytes of GX buffaloes using a more sensitive SYBR Green dye-based RT-qPCR assay. The expected size of the amplified products was verified by electrophoresis (Fig. 1b). The results of tissue-specific TLRs 1–10 expression levels are shown in Fig. 4. Our expression analysis of leukocytes showed that relative TLR4 expression was significantly higher than any other TLR in leukocytes (p < 0.0001), and then the expression of TLR4 in all other examined tissues (p < 0.0001). Other TLRs in leukocytes showed low expression levels, among which TLRs 5, 6 and 9 had the lowest expression levels.
Among all examined organs and tissues, TLR1 expression was highest in HLNs and lowest in the lungs. The highest expression level of TLR2 was in the heart (p < 0.0001) and the lung (p < 0.01), and lowest expression was in the liver. TLR2 expression in SILN (p < 0.05), heart (p < 0.0001), lung (p < 0.01), and kidney (p < 0.05) was significantly higher than any other TLR in the same group, except TLR1 and TLR10 in the kidney. TLR3 had relatively mild or moderate expression levels in all tissues. TLR4 expression in spleen was significantly higher than TLR5 (p < 0.01), TLR6 (p < 0.05), TLR9 (p < 0.01), and TLR10 (p < 0.05). In MLN and spleen, higher expression values of TLR4 (p < 0.01) were found, as well as relatively higher expression of TLR7 and TLR8 mRNAs, compared to all other tested tissues.