Evolution of the chicken Toll-like receptor gene family: A story of gene gain and gene loss

Aims

The major aims of this investigation were to determine what TLRs are present in chicken and determine where unmapped TLR genes are located. These data then allowed us to define orthologs in other species and examine TLR gene family evolution using phylogenetic methods. These analyses allowed us to characterise the variation within and between vertebrate lineages, as well as the dynamics of TLR gene gain/gene loss. In addition, TLR gene nomenclature issues were also addressed.

Definition of orthologs and TLR gene nomenclature

TLR gene discovery in chicken and other vertebrates

The chicken TLR repertoire contains tandem duplications of the genes for TLRs 1L and 2 [31, 17]. A partial sequence of chicken TLR1LB has been described previously by Yilmaz et al. [17]. Using RT-PCR we have sequenced and cloned the full length mRNA [Genbank: DQ518918]. This sequence includes two further exons which, although transcribed, are not translated. Using synonymous substitution rates (see Materials and Methods) we estimated the gene duplication events that gave rise to TLR1LA, 1LB and 2A, 2B to have taken place 147 and 66 Mya respectively. The chicken is also found to be missing a number of TLRs, which are present in most mammals. The TLR7, 8 and 9 subfamily is present in fish and mammals but is only represented by TLR7 in the chicken [14, 16, 17]. TLR8 is present as a pseudogene [32] and TLR9 has been deleted.

Genetic Mapping of chicken TLR1LA and B genes

The chromosomal locations of most chicken TLRs are known [17], with the exception of TLR1LA and TLR1LB which, based on a single sequence contig [Genbank: NW_001471687] map close together on an unmapped chromosome [17]. Sequencing of a genomic fragment from TLR1LB locus from parents of the East Lansing reference mapping cross, identified two SNPs in TLR1LB. Genetic linkage analysis of TLR1LB showed linkage to chromosome 4 (see Table 3).

Structural analysis of proteins and selective constraints on functional domains

Comparisons of structural predictions for the chicken TLRs show them to have a varied structure [17]. The newly discovered TLR15 has the unusual feature of many LRRs clustered towards the C-terminus of the molecule and few at the N-terminus (Figure 1). TLR7 has the uncharacteristic feature of a splice variant with two predicted trans-membrane (TM) domains only previously found in mammalian TLR6 (see Additional file 2). The predicted structure of the TLR1LB fragment described by Yilmaz et al. [17] is shown to have a further predicted TM domain towards the N-terminus whereas the majority of TLRs have a single TM domain between the extracellular and intracellular domains (Additional File 2). The predicted structure of this TLR is otherwise very similar to TLR1LA only missing the two LRRs at the N-terminus (see Figure 1). By contrast the duplications of TLR2 are marked by the absence of an LRR from TLR2A, which is present in all other known TLR2s (Additional File 2).

Phylogeny of vertebrate TLRs

Molecular phylogenies based on sequence alignments are only as accurate as the alignment data from which they are produced; consequently it is important that the alignment data quality is determined [33–35]. In order to assess the alignment's tree like structure, likelihood mapping was carried out and showed the sequence alignment data to have strong phylogenetic signal (Figure 2). Both phylogenetic methods used produced near identical topologies, which were strongly supported with bootstrapping (see Figures 3 and 4). Close up images of individual clades in Fig. 3 are shown in Additional files 4, 5 and 6. Close up images of individual clades in Fig. 4 are shown in Additional files 7, 8 and 9.

The majority of TLRs form specific clades, with each species containing one copy of each receptor (e.g. TLRs 3, 4 and 5). The main exception to this trend is the assemblage of receptors named TLR1, which have independently evolved in fish, Xenopus, chicken, opossum and eutherian mammals (see Figures 3 and 4). These genes are not orthologs but have often been given the same name, which creates some confusion (see Table 1). The same is true for TLR5S in fish, where a soluble TLR5 has independently evolved in different fish species [36] (see Figures 3 and 4). The presence of most TLR gene families across so great a range of vertebrate taxa indicates that the majority of TLR families must have been established prior to the evolutionary radiation of vertebrates in the Cambrian Period 540-505 Mya [37].

Phylogenetic analyses show that chicken TLR21 is an ortholog of fugu and Xenopus TLR21. Chicken TLR15 shares a clade with TLRs 1, 2, 6, 10 and 14 (Figures 3 and 4), but the precise nature of this relationship remains unresolved due to the small number sequences available. TLRs 3, 5, 7, 8 and 9 form a clade that mainly follow the expected species tree; the exceptions to this being the loss of chicken TLRs 8 and 9 and the gain of soluble versions of TLR5 in teleost fish. Gene gain has been comparatively common in the clade consisting of TLRs 2, 14 and the 1, 6, 10 assemblage. The TLRs in this clade are found to respond to bacterial pathogens, and in mammals and chicken have been shown to function heterodimerically [3, 38].

No specific ligand has been identified for the clade containing TLR 13, 21 and 22 (Table 2). The neighbour joining tree (Figure 3) attaches the TLR 11–12 clade to this group. The positioning of this clade is not as well supported by neighbour joining as it is using maximum parsimony, which places the TLR11–12 clade basal to the TLR3, 5, 7, 8 and 9 group (Figure 4). Both the TLR11–12 clade and TLR13-21-22 clade share a history of multiple gene loss events. TLRs 14, 22 and 23 have been found to have been lost from all land animals (see Figures 3 and 4).

Comparison of the predicted structure of Toll-like receptors

From protein alignments as well as SMART predictions it is evident that the TIR domain is highly conserved across many different receptor families whereas the extracellular LRR region shows great variation in the positioning of LRRs between receptors although a strong level of structural conservation across taxa for a given receptor (see Additional File 2).

Chicken TLR7 has a predicted TM domain at the C-terminus as have the majority of TLRs, but also has an extra TM domain at the N-terminus. Similarly an extra TM domain at the N-terminus tip is predicted for Chicken TLR1LB as well as TLR6, a receptor present only in eutherian mammals. A TM domain is predicted at the N-terminal in fugu TLR1 but in this case the usual TM domain is not found at the C-terminal (Additional File 2). As the predicted N-terminal TM domain of TLR 1LB is also found in a number of related TLRs it is possible that this is a conserved feature. However, the absence of such a predicted structure from all examples of TLR7 other than chicken may mean that this feature is a result of the SMART program falsely interpreting a hydrophobic region. A pseudogene caused by an insertion event for TLR8 in chicken has been identified but TLR9 in chicken has been lost leaving only flanking regions where the TLR9 locus is predicted [32].

Phylogenetic analyses of the Toll-like receptor gene family

There have been previous phylogenetic analyses of TLR data [4, 10, 15, 17]. Since the inclusion of more sequences improves the data quality of the alignment and consequently the phylogenies produced [19, 39], this investigation has been able to improve on phylogenetic analyses carried out by previous studies owing to the increased number of sequences available. This is corroborated by the likelihood mapping analysis (Figure 2) and the high bootstrap values by both neighbour joining (Figure 3) and parsimony (Figure 4) methods for the major TLR family groups.

Also very similar trees are produced from neighbour joining and maximum parsimony analyses (Figures 3 and 4).

Gain of Toll-like receptor genes

The vast majority of duplications of whole TLR genes occurring since the Cambrian period have occurred in the clade containing TLRs 1, 2, 6 and 10. The first of these took place 300 Mya leading to the establishment of TLR 10 in the lineage that gave rise to the mammals. This was followed by TLRs 1 and 6 in the eutherian mammals 130 Mya [4]; data which are supported by Figures 3, 4 and Roach et al. [15]. A further duplication has occurred 147 Mya in the lineage that gave rise to modern birds, leading to avian TLR1LA and TLR1LB followed by the latest duplication 65 Mya producing TLR2A and TLR2B in the chicken. The duplication of TLR2 in the chicken lineage is not a unique event, other independent duplication events of TLR2 have occurred in X. tropicalis (Figures 3 and 4), the American alligator [4] and the ancestor of marsupial and eutherian mammals [15] the timing of these duplications, however, remains unknown.

All of the LRRs of TLR1LA are in positions shared by closely related TLRs (Additional File 2). The much shorter TLR1LB has an almost identical 3' sequence and structure, but the two 5' exons and the start of the third and final exon are disrupted by stop codons and frame shifts. A start codon is found a short way into the third exon producing a product lacking the two N-terminal LRRs and including an extra predicted TM domain. The predicted structure of this receptor is the same as that reported in Yilmaz et al. [17] although the complete coding sequence of this gene is now complete (this paper) [GenBank: DQ518918].

Loss of Toll-like receptor genes

Compared to the paucity of gene gain in the evolution of the TLR family there are many instances of gene loss. TLRs 14, 22 and 23 are not found in any land dwelling vertebrates and TLR4 has been lost in the majority of studied fish species, yet is preserved in land animals. TLRs 11, 12 and 13 found in a small number of animals are lost in birds and many mammals, however a pseudogene for human TLR12 has been reported [30]. Although TLR19 of the zebrafish and TLR16 of X. tropicalis clearly represent ancestral sequences of TLRs 11 and 12 (see Figures 3 and 4), there is no trace of a receptor of this kind to be found either in other fish, mammals or the chicken. This indicates that members of this group must have been lost on more that one occasion. Functional genes for TLRs 8 and 9 are absent from the chicken; TLR8 is present as a pseudogene [32] and TLR9 has been deleted from the genome.

The inclusion of certain taxonomic groups would assist in our understanding of these gene losses. The inclusion of fish species (other than the ray-finned fish) and more basal tetrapods (squamata, anapsida etc.) would help to give a more precise idea of the point at which TLRs 22, 14 and 19 were lost.

Sequence databases

Candidate chicken TLR genes were identified using extensive searches for sequence homology using BLAST started by using the mouse and human sequences identified in Smith et al. [16] and the TIR protein consensus sequence from Meijer et al. [10] in the NCBI database [40] and the Ensembl genome browser [41]. Full coding region nucleotide and corresponding protein sequences were collected for as great a taxonomic range as possible as this improves the validity of the phylogenetic analysis [39]. Collected sequences were used in BLAST searches both into specific genomes and the general sequence database with the intention of identifying all available TLR sequences.

The full sequence for chicken TLR1LB was determined using the original version of partial sequence Genbank: AY633573, the contig Genbank: NW_001471687 (the details described in Yilmaz et al. [17]) and EST Genbank: CD762233. The complete sequence was confirmed by sequencing cDNA fragments produced by RT-PCR (for list of primers see Additional File 3). To ensure that all chicken TLRs had been identified, Hidden Markov models (HMMs) were used to analyse the chicken genome [42]. When expected TLR homologues were found to be absent from the chicken (e.g. TLR9), the sequence contig for the expected chromosomal region was identified to ensure that the absence of a TLR was not the result of a gap in the genome sequence.

Bioinformatics

Complete protein sequences were assembled and labelled in BioEdit [43] before being aligned using ClustalX [44]. This program automatically produces a phylogeny that can be viewed in TreeView [45], which can be used as a rough guide although it is not suitable for in depth phylogenetic investigation [33]. The alignment was checked by visual inspection. To ensure that the alignment dataset was suitable for phylogenetic analyses likelihood mapping was carried out using the program TreePuzzle [46]. This method, although little used, is of great value for confirming the suitability of a dataset for use in tree building [35]. Likelihood mapping works by breaking down the alignment into quartets of sequences and assessing the likelihood of their three possible phylogenetic topologies. This information is displayed as points in an equilateral triangle where points located in the apexes represent strong 'tree like' phylogenetic signal, points between apexes give equal likelihood for two different tree topologies and points in the centre indicating 'star like' evolution [35, 46, 47].

A convenient way to infer events of gene duplication and gene loss is to note to what extent the gene phylogenies match the expected species phylogeny [20, 21, 48]. Incongruence in the close association between the two can usually be explained by episodes of gene duplication and/or gene loss in a given lineage [18, 20, 21, 34, 48]. Gene duplication followed by positive natural selection is seen as the principle means by which new gene functions arise [26, 49]. A study relating to the evolution of gene families on the chordates found the majority of duplications to have occurred early in chordate evolution. This can be easily maintained by the similarity between the orthologous genes in different species being greater than that found between members of the same gene family in the same species [50, 51].

Phylogenies were created in PAUP [52] using both parsimony and neighbour-joining methods. These two methods are used as each one analyses the data in different ways [19, 33]. Consequently, if tree topology is not shared by both methods it shows that these relationships may not represent evolutionary heritage. To determine the support for individual clades, 1000 bootstrap replicates were carried out on the trees produced by each method [53].

The structure of chicken TLRs was predicted using the SMART program [54] to analyse full length amino acid sequences.

RNA extraction and tissue sources

Tissues were collected and stored at -80°C from a Gallus gallus broiler presumed to be healthy. All tissues were ground under liquid nitrogen with a pestle and mortar. 100 mg samples of each tissue were trizol extracted and RNA quantified using a bench top spectrophotometer. To ensure RNA extractions were free from genomic DNA contamination, RNase free DNase (Promega) was used in accordance with the manufacturer's instructions. RNA samples were reverse transcribed using SuperScript II reverse transcriptase (Invitrogen Carlsbad CA) following the manufacturers guidelines, including the use of RNase H to remove complimentary RNA from the newly synthesised cDNA. All cDNAs were standardised using RT-PCR against the reference gene, GAPDH and the absence of contamination determined by the absence of bands in samples that had not undergone reverse transcription.

Reverse-transctriptase PCR (RT-PCR) and PCR

PCR was carried out in 10 μl volumes consisting of 10 pg template, 10× dNTPs (2 mM each), 10× PCR buffer, 5× GC solution, 0.4 U Taq and 5 pM each primer. List of primer sets used is shown in Additional file 3. PCR conditions for all reactions were as follows: 15 min 95°C followed by 30 cycles of 95°C for 30 seconds, 60°C for 30 seconds and 72°C for one minute followed by an extension step of 72°C for 5 minutes. Amplified PCR products were run by electrophoresis in a 1% agarose gel stained with ethidium bromide and visualised using an ultra violet transilluminator.

DNA sequencing

As the sequence of TLR1LB was unclear from sequence files and EST data, the full sequence was confirmed using primers shown in Additional File 3. To improve sequence quality of TLR1LB TOP sequence (see Additional File 3) we collected PCR products from agarose gels using the QIAEXII gel extraction kit (Qiagen Hilden). Inserts were grown in bacterial colonies on agar plates before colonies were inoculated into LB medium and plasmids purified using the QIAprep Miniprep kit (Qiagen Hilden). For all other sequences (those used in the SNP analysis) PCR products were purified using ExoSap and sequenced in both directions using the BigDye^® Terminator v3.1 Cycle sequencing kit.

Genetic linkage mapping of Toll-like receptor genes

Two single nucleotide polymorphisms (SNPs) were identified in the TLR1LB expression fragment (for primer list see Additional File 3) from Red Jungle Fowl and White Leghorn parent birds [55]. Sequences from offspring were sequenced and genotyped and these data used for genetic mapping. TLR1LA is known to map closely to TLR1LB [Genbank: NW_001471687]. Consequently this information can map both genes. The map locations of other chicken TLRs was determined using cDNA sequences in a BLAT search on the UCSC genome browser [56]. The sequence for chicken TLR2B was shown to contain a gap of unidentified nucleotides. This region was sequenced in both versions of TLR2 and the gap found to be an assembly artefact on both the February 2004 and May 2006 genome builds [Genbank: EF595650].

Calculation of time of divergence of duplicated genes

Nucleotide alignments of chicken TLR duplications were created using DAMBE [57] which aligns nucleotide sequences to the protein alignment one codon at a time, thus allowing functionally similar areas of each gene to remain aligned and not become victim of codon slippage. Non-synonymous (dN) and synonymous sequence divergence (dS) were estimated with the codeml program in PAML [58] in Runmode = -2 with codonfreq = 2 using the Goldman and Yang [59] substitution model. The substitution rates were estimated using the nucleotide frequencies at each codon position using the F3 × 4 model [58] to estimate the expected codon frequencies.

In order to get an estimate of the age of these duplications (T) the equation, r = K/(2T) [24] was used, where r is the neutral substitution rate per site per year and K is the total synonymous sequence divergence (d_S) for the total aligned sequence (0.4119 and 0.184 for TLRs 1 like and 2 respectively). We used an overall estimate of r = 1.4 × 10^-9 that we obtained by dividing the autosomal sequence divergence between chicken and turkey (0.1008) from Axelsson et al. [60] with the time of divergence between chicken and turkey of 28 × 10⁶ – 45 × 10⁶[61, 62].