Toll-Like Receptor Evolution in Birds: Gene Duplication, Pseudogenization, and Diversifying Selection.

Toll-like receptors (TLRs) are key sensor molecules in vertebrates triggering initial phases of immune responses to pathogens. The avian TLR family typically consists of ten receptors, each adapted to distinct ligands. To understand the complex evolutionary history of each avian TLR, we analyzed all members of the TLR family in the whole genome assemblies and target sequence data of 63 bird species covering all major avian clades. Our results indicate that gene duplication events most probably occurred in TLR1 before synapsids diversified from sauropsids. Unlike mammals, ssRNA-recognizing TLR7 has duplicated independently in several avian taxa, while flagellin-sensing TLR5 has pseudogenized multiple times in bird phylogeny. Our analysis revealed stronger positive, diversifying selection acting in TLR5 and the three-domain TLRs (TLR10 [TLR1A], TLR1 [TLR1B], TLR2A, TLR2B, TLR4) that face the extracellular space and bind complex ligands than in single-domain TLR15 and endosomal TLRs (TLR3, TLR7, TLR21). In total, 84 out of 306 positively selected sites were predicted to harbor substitutions dramatically changing the amino acid physicochemical properties. Furthermore, 105 positively selected sites were located in the known functionally relevant TLR regions. We found evidence for convergent evolution acting between birds and mammals at 54 of these sites. Our comparative study provides a comprehensive insight into the evolution of avian TLR genetic variability. Besides describing the history of avian TLR gene gain and gene loss, we also identified candidate positions in the receptors that have been likely shaped by direct molecular host-pathogen coevolutionary interactions and most probably play key functional roles in birds.

Introduction

The species richness and wide range of ecological adaptations are remarkable in birds (Jetz et al. 2012). Despite being similar to mammals in many aspects of their biology, birds evolved many similar traits (including homoiothermy and parental care) convergently through their analogous ecological life strategies (Farmer 2000; Emery and Clayton 2004; Olkowicz et al. 2016; Lovegrove 2017). Recent advance in avian genomic research, which started with the Avian Phylogenomics Project (http://avian.genomics.cn/en/; Zhang, Jarvis, et al. 2014) and continues with the Bird 10,000 Genomes Project (https://b10k.genomics.cn/) has created new possibilities to use publicly available whole genome sequencing data to resolve many questions regarding avian biology, evolution, and adaptations (Zhang, Li, et al. 2014). Among these, the evolution of avian immune function is of particular interest. Being challenged by similar classes of pathogens as found in mammals, highly diversified immune genes in birds are rich in mixed signatures of symplesiomorphy, convergence, and apomorphic adaptations to novel functions (Burri et al. 2010; Cheng et al. 2015).

The coevolution between host and pathogen likely involves mainly molecules that form direct interface between the host and the pathogen structures, for example, the pattern recognition receptors (PRRs) and their ligands, the pathogen-associated molecular patterns (PAMPs; Janeway and Medzhitov 2002). Forming the first line of the host immune defence, PRRs appear to constantly evolve toward specific and appropriate recognition of certain PAMPs (Wang, Zhang, Liu, et al. 2016). The Toll-like receptor (TLR) protein family belongs to one of the most essential and functionally most characterized PRRs (Palm and Medzhitov 2009; Coscia et al. 2011). After specific PAMP binding, TLRs trigger signaling pathways activating transcription factors, such as NF-κB to induce expression of target genes that are key to triggering an inflammatory immune response and subsequent activation of acquired immunity (Iwasaki and Medzhitov 2015).

TLRs are transmembrane proteins each with the characteristic horse-shoe-shaped ectodomain (ECD), where the direct contact between the receptor surface and specific microbe molecules occurs, the transmembrane domain and an intracellular Toll–interleukin 1 receptor (TIR) domain that enables the downstream signal transmission (Kawai and Akira 2010). The capacity to bind a particular type of ligand in TLRs is dependent on the structure of the ECD. This has been used to divide the vertebrate TLRs into three different groups based on the pattern of hydrogen bonds formed among the asparagine residues in the neighboring leucine-rich repeats, which stabilize the shape of the ECD (Wang, Zhang, Liu, et al. 2016). The ten avian TLR family members belong either to the single-domain TLRs possessing a complete asparagine ladder (TLR3, TLR5, TLR7, TLR15, and TLR21) or to the three-domain TLRs with the ladder interrupted in the central part of the ECD (TLR1A, TLR1B, TLR2A, TLR2B, and TLR4; Wang, Zhang, Liu, et al. 2016). While some avian TLRs conservatively recognize the same ligands as mammalian TLRs (such as TLR4 detecting bacterial lipopolysaccharides, LPS; TLR5 flagellin; TLR3 dsRNA or TLR7 ssRNA; reviewed by Keestra et al. 2013), others were reported to form distinct paralogues with related ligand specificity (e.g., heterodimer-forming TLR1A/TLR1B together with TLR2A/TLR2B, which recognize di/triacylated lipopeptides; Keestra et al. 2007; Higuchi et al. 2008) or achieve recognition of similar ligands as their mammalian analogues through convergence (e.g., avian TLR21 that recognizes CpG DNA similarly as mammalian TLR9; Brownlie et al. 2009; Keestra et al. 2010). Finally, TLR15 that is unique to birds evolved to gain a novel function in recognition of extracellular proteases (de Zoete et al. 2011).

Gene duplication is an important mechanism allowing genes to evolve novel functions (Zhang 2003; Ellegren 2008). Gene duplications are common in TLRs, as with other immune genes. Duplicated TLR1 family members are known in mammals (TLR1, TLR6, TLR10, and TLR2) as well as in birds (TLR1A, TLR1B, TLR2A, and TLR2B; Coscia et al. 2011; Huang et al. 2011; Wang, Zhang, Liu, et al. 2016). Recently, duplication of TLR7 has been described in several avian taxa (Cormican et al. 2009; Grueber et al. 2012; Raven et al. 2017), while other birds possibly lack the duplicated TLR7. This intriguing finding contrasts with that found in mammals, which contain two other TLR7-like family members, TLR8 and TLR9, both missing in birds (Philbin et al. 2005; Temperley et al. 2008), suggesting that gene loss may also play a significant role in avian TLR evolution. TLR5 pseudogenization occurred independently several times within the passerines (Bainova et al. 2014) and possibly in parrots (Alcaide and Edwards 2011). The general pattern of TLR pseudogenization and gene duplication is largely unknown and may be far more common across various TLR genes and avian taxa.

Although generally conserved in their structure, TLRs were reported to exhibit very high levels of inter- and intraspecific genetic variation in birds (Alcaide and Edwards 2011; Huang et al. 2011; Grueber et al. 2014). Several authors have documented the functional significance of TLR genetic variation (Leveque et al. 2003; Walsh et al. 2008), characterized associations between TLR variation and disease susceptibility (Netea et al. 2012; Medvedev 2013) and suggested the maintenance of intraspecific polymorphism through balancing selection mediated by pathogens (Ferrer-Admetlla et al. 2008). Despite several attempts to understand the evolutionary significance of TLR genetic variation in birds (Alcaide and Edwards 2011; Grueber et al. 2014; Vinkler et al. 2014), the question of its adaptive value for avian specific PAMP recognition remains unresolved. Here, we use various approaches to analyze natural selection in the most comprehensive attempt to investigate adaptive pathogen-driven evolution in avian TLRs.

Facilitated by the conservative structures of TLR molecules, they present a suitable model for studying pathogen-driven microevolutionary processes at the DNA level. The direct contact between TLRs and PAMPs is restricted to certain amino acid sites (Gay and Gangloff 2007) where even single amino acid substitutions may have profound effects on receptor binding properties (Keestra et al. 2008; Walsh et al. 2008; Resman et al. 2009; Meng et al. 2010). Thus, positively selected sites can be predicted to emerge from the generally negatively selected background and can be used for prediction of functionally relevant positions in nonmodel animals lacking precise protein crystallographic data. Several studies have investigated the nature of positively selected sites (PSS) at the interspecific level in various TLRs in selected bird taxa (Alcaide and Edwards 2011; Huang et al. 2011; Grueber et al. 2014; Vinkler et al. 2014; Wang, Zhang, Liu, et al. 2016). Yet, colocation study of PSS on the receptor surface with the predicted functional residues already known in the model species (human and mouse), has never been done across the avian phylogeny (but see Vinkler et al. 2014, for a similar approach in Galloanserae birds).

Given the current lack of comprehensive evidence on TLR family molecular evolution across the avian phylogeny, our comparative study reported here performed evolutionary analysis on whole-genome sequences for 48 species representing 34 avian orders (supplementary material S2: table S23, Supplementary Material online) mainly gained through the Avian Phylogenomics Project (Zhang, Jarvis, et al. 2014). The TLR sequences extracted from these genomes together with those added from other public resources (in total 63 species) allowed us to infer the history of avian TLR gene gain and gene loss in the context of TLR evolution in other vertebrates. By critical assessment of positive selection acting on all members of the TLR family in birds, we describe the adaptive microevolutionary changes in these immune receptors in their molecular context. Thus, this study is the first to comprehensively predict functionally relevant genetic variation in avian TLRs providing insights into the coadaptation between host and pathogen through evolution of ligand recognition.

Results and Discussion

Evolution of TLR1 Family Gene Duplications

The TLR1 gene family is widely duplicated in vertebrates including birds (TLR1A, TLR1B, TLR2A, and TLR2B; Temperley et al. 2008; Cormican et al. 2009; Huang et al. 2011; Wang, Zhang, Liu, et al. 2016). Presently, however, the timing of these gene duplication events is unclear: the paralogues within the TLR1 and TLR2 subfamily were suggested to have duplicated before the mammalian divergence from sauropsids (Huang et al. 2011) or after this divergence, independently in birds and mammals (Temperley et al. 2008; Cormican et al. 2009; Mikami et al. 2012; Wang, Zhang, Liu, et al. 2016). This disagreement results from the lack of accounting for gene conversion as a mechanism of paralogue sequence homogenization in some studies (Temperley et al. 2008; Cormican et al. 2009; Wang, Zhang, Liu, et al. 2016). We found significant support for gene conversion between the two paralogues both in the TLR1 subfamily (TLR1A and TLR1B) and TLR2 subfamily (TLR2A and TLR2B; supplementary material S1: table S1, Supplementary Material online; for detailed info see supplementary material S3, Supplementary Material online) in birds. Therefore, we constructed both TLR1 and TLR2 phylogenetic trees based on the protein sequences of the nonconverted regions only (fig. 1; supplementary material S1: figs. S1 and S2, Supplementary Material online).

. 1.

Phylogenetic tree based on nonconverted regions of TLR1 subfamily members. The bootstrap values of maximum likelihood analysis obtained using PhyML and the posterior probability of Bayesian analysis obtained using MrBayes (in percentage per each node) are provided. Birds are represented by zebra finch (TaGu, Taeniopygia guttata) and chicken (GaGa, Gallus gallus), crocodiles by alligator (AlMi, Alligator mississippiensis), turtles by painted turtle (ChPi, Chrysemys picta), mammals by human (HoSa, Homo sapiens) and horse (EqCa, Equus caballus), amphibians by clawed frog (XeTr, Xenopus tropicalis), bony fish by zebrafish (DaRe, Danio rerio) and cartilaginous fish by shark (ChGr, Chiloscyllium griseum). The analysis was performed using a single amino acid sequence per TLR and species. Based on the results, we suggest renaming TLR1A to TLR10 and TLR1B to TLR1 in birds.

Our results based on the analysis of 179 sequences of 83 vertebrate species (for the list see supplementary material S2: table S24, Supplementary Material online) show that in the TLR1 family the sequences of avian species cluster together based on paralogue identity (supplementary material S1: fig. S1, Supplementary Material online). Avian TLR1A and reptilian TLR10 cluster together with mammalian TLR10, while avian TLR1B and reptilian TLR1/TLR6 (inconsistent nomenclature in reptilian TLR1) cluster together with mammalian TLR1 and TLR6 paralogues (shown in fig. 1; detailed in supplementary material S1: fig. S1, Supplementary Material online). This is consistent with the results previously published by Huang et al. (2011) confirming that the first duplication event within the TLR1 subfamily occurred before mammal–reptile divergence. On the basis of these findings, we suggest renaming avian TLR1A to TLR10 (from now on marked as TLR10 [TLR1A] in the text), while modifying the name of the avian TLR1B gene to TLR1 (marked as TLR1 [TLR1B]). Our analysis also indicates that the duplication of TLR1 in amphibians (namely in Xenopus) apparently arose independently of the TLR1 duplication in amniotes (supplementary material S1: fig. S1, Supplementary Material online).

In addition, in the TLR2 subfamily there are two copies present in birds (TLR2A and TLR2B). In mammalian genomes, on the other hand, only one functional TLR2 is maintained and for some species TLR2-like pseudogenes (a second copy) were described (Roach et al. 2005). The duplication of TLR2 in sauropsids could have arisen independently of mammals, or alternatively, a duplication event predating divergence between birds and mammals could have occurred (similar to the TLR1 subfamily), followed with pseudogenization leaving only a single TLR2 copy in mammals (Huang et al. 2011). Our analysis of 137 homologous sequences in 78 vertebrate species (for the list see supplementary material S2: table S24, Supplementary Material online) contradicts the second scenario (supplementary material S1: fig. S2, Supplementary Material online) and supports the independent duplication of TLR2 in sauropsids. This is also supported by the phylogenetic analysis of a ca. 80 amino-acid-long alignment of the conserved sequence in human TLR2/TLR2-like pseudogene and chicken TLR2A/TLR2B (supplementary material S1: fig. S3, Supplementary Material online); the resulting phylogenetic tree suggests that TLR2 duplication in birds was independent of a parallel mammalian TLR2 duplication that was followed by a pseudogenization of the second TLR2 copy (supplementary material S1: fig. S4, Supplementary Material online). However, it is difficult to draw any final conclusions on the validity of the two scenarios proposed for TLR2 evolution in amniotes since the support for the branching of the avian TLR2 paralogues is limited (bootstrap support for the mammal-sauropsid split = 47, supplementary material S1: fig. S2, Supplementary Material online).

Avian TLR7 Gene Duplications

Birds lack two of the three vertebrate TLR7 family members (i.e., TLR8 and TLR9; Philbin et al. 2005; Temperley et al. 2008). Instead, the avian TLR7 locus appears to have recently duplicated in some passerines (Cormican et al. 2009; Grueber et al. 2012) and waders (Raven et al. 2017). Within the genome sequences investigated in this study, we also found two TLR7 copies in other avian taxa: Cuculiformes and Mesiornithiformes (fig. 2). The number of amino acid substitutions distinguishing the two TLR7 loci ranges between species from 21 in zebra finch (Taeniopygia guttata) to only 3 in rifleman (Acanthisitta chloris; see fig. 2; supplementary material S1: table S2 and supplementary material S4, Supplementary Material online). Furthermore, the two TLR7 copies likely differ slightly in their functions, because 11 of the sites identified as distinct between the two copies in any avian species match the known ligand-binding positions (Wei et al. 2009; Gupta et al. 2016; Zhang et al. 2016), five are situated in the Z-loop region responsible for ligand binding and dimerization and the other 10 variable sites are identical with PSS detected in birds (supplementary material S1: table S3, Supplementary Material online). To bring independent support to the evidence of TLR7 gene duplication, we performed a PCR-based quantitative copy number variation (CNV) analysis comparing TLR7 with two single-copy TLRs (TLR3 and TLR4) in selected species (collared flycatcher, Ficedula albicollis; ground tit, Pseudopodoces humilis; white-throated sparrow, Zonotrichia albicollis and zebra finch). Our results confirm that the zebra finch (a species with two TLR7 loci based on the genome sequence analysis) possesses two TLR7 copies (supplementary material S1: tables S4–S6 and fig. S5, Supplementary Material online), while the other investigated species (with a predicted single TLR7 locus) possess only one TLR7 gene copy. In contrast to the TLR1 family, we did not detect any gene conversion between the two copies of the duplicated TLR7 genes (see supplementary material S3, Supplementary Material online). Despite this, the constructed tree of avian TLR7 does not show any evidence of two separate TLR7 phylogenetic clusters, but surprisingly both copies of TLR7 always clustered separately for each species (supplementary material S1: fig. S6 and table S7, Supplementary Material online). Therefore, we infer that the duplication in TLR7 occurred several times independently in recent avian history and possibly even multiple times in some cases (as, e.g., in ruddy turnstone, Arenaria interpres, where TLR7 was reported to be triplicated; Raven et al. 2017).

. 2.

Avian TLR7 duplication. In the schematic avian phylogenetic tree the species with duplicated TLR7 are highlighted in teal rectangles. The numbers of amino acid substitutions distinguishing the two copies of the duplicated TLR7 are shown behind the species name. The analysis was performed using a single sequence per gene and species.

TLR5 Pseudogenization in Birds

Recently, it has been shown that there is no functional gene for TLR5 in some passerine species (Bainova et al. 2014). In this study, we have also found several avian taxa possessing only nonfunctional TLR5 pseudogenes with stop-codons in their sequences. Apart from passerines, we revealed TLR5 pseudogenes in other clades: that is, Psittaciformes, Cariamiformes, Trogoniformes, Phaethontiformes, Eurypygiformes, and Apodiformes (fig. 3). The sequence data show that the stop-codon positions differ among the investigated species with differences also in the pseudogenization mechanism (single nucleotide substitution or frame-shifting indels, see supplementary material S1: table S8, Supplementary Material online). The stop-codons in TLR5, thus, arose independently in the evolutionary history of the different avian taxa. This, however, does not appear to be the result of relaxed selection acting in TLR5. Firstly, the number of single nucleotide variants (SNVs) in the alignment of the TLR5 functional sequences is comparable to other TLRs (indicating comparable mutation rate; supplementary material S1: table S9, Supplementary Material online). Secondly, the total number of PSS is higher than in other TLRs and ω (dN/dS ratio) is again similar to other TLRs (see below and supplementary material S1: table S10, Supplementary Material online) and higher than the genome average (Zhang, Li, et al. 2014). However, since most of the PSS contain only conservative substitutions without dramatic effects on the amino acid site physicochemical properties, it appears that negative selection may prevent the loss of function in those species where the TLR5 functional gene is maintained. Altogether this suggests that there may be a selection for TLR5 loss of function in certain avian evolutionary lineages. The loss of a TLR5 functional gene is not limited only to birds. The TLR5 pseudogene has also been described in one human allele, where its presence is associated with an increased risk of pneumonia infection caused by flagellated bacteria Legionella pneumophila (Hawn et al. 2003; Zhang et al. 2013), while at the same time possessing a TLR5 nonfunctional allele might be advantageous for decreasing the probability of autoimmune disease development (Hawn et al. 2005). TLR5 pseudogenization may be allowed by the high redundancy of pathogen detection, where apart from TLR5, flagellin is also recognized by other PRRs, for example, NLRC4 inflammasome (Zhao et al. 2011; Yang et al. 2014). In birds other flagellin-recognizing PRRs have not yet sufficiently been studied to support this hypothesis.

. 3.

TLR5 pseudogenization in birds. Species possessing only TLR5 pseudogene are highlighted in red rectangles within the schematic representation of avian phylogenetic tree. The position of the first stop codon is indicated by the number provided behind the species name (position numbering according to the chicken reference). The analysis was performed using a single sequence per gene and species.

Diversifying Selection in Avian TLRs

We tested pervasive, positive, and diversifying selection acting in all avian TLRs. Numbers of PSS differ between TLRs (from 0.3% to 5.8% PSS per TLR; table 1; supplementary material S1: fig. S7 and supplementary material S5: table S26, Supplementary Material online). The results suggest that positive selection is acting more on TLRs exposed toward the cell surface (mainly TLR1 [TLR1B], TLR2A, TLR5, TLR4, and TLR2B), than in the TLRs situated in endosomes (TLR21, TLR3, and TLR7). This may be because the endosomal TLRs are specialized for detection of less complex ligands, which show low structural variation (such as ssRNA in TLR7, dsRNA in TLR3 or CpG DNA regions in TLR21; Brownlie and Allan 2011) interacting with the genetic variability in TLRs. This supports the previous findings by Mikami et al. (2012) in vertebrates. The only exception is TLR15, which is unique to birds and reptiles (Boyd et al. 2012). Although situated on the cell surface, TLR15 is activated (unlike other TLRs) by the ECD proteolytic cleavage with pathogen-derived proteases (de Zoete et al. 2011). Since TLR15 harbors only low numbers of PSS (table 1), we hypothesize that this pathogen-recognition mechanism may be linked to reduced positive selection acting on the ECD. Nonetheless, the cleavage site is variable across bird taxa and contains two PSS (supplementary material S1: table S11, Supplementary Material online; see also supplementary material S1: text S1 and fig. S8, Supplementary Material online).

Table 1.

The Number of Positively Selected Sites in Avian TLRs.

	Speciesa	aa lengthb	PSSc	PSS/TLRd(%)
TLR10 [TLR1A]	45	794	31	3.9
TLR1 [TLR1B]	44	638	37	5.8
TLR2A	40	793	43	5.4
TLR2B	42	781	35	4.5
TLR3	51	895	19	2.1
TLR4	54	843	38	4.5
TLR5	46	861	45	5.2
TLR7	51	1041	31	3.0
TLR15	53	868	24	2.8
TLR21	14	907	3	0.3

Number of species (one sequence per species).

The protein sequence length in reference chicken TLRs (NCBI accession numbers are listed in supplementary material S1: table S14, Supplementary Material online).

Number of positively selected sites detected in investigated species per each gene.

The percentage of PSS per whole receptor amino acid sequence.

Our results are also mostly consistent with the findings of Wang, Zhang, Liu, et al. (2016), who suggested grouping of TLRs based on their ECD architecture and showed that in mammals the single-domain TLRs (TLR3, TLR5, TLR7, TLR15, and TLR21) are under stronger purifying selection than the three-domain TLRs (TLR10 [TLR1A], TLR1 [TLR1B], TLR2A, TLR2B, and TLR4). Here, we found that also in birds the positive selection is acting more on the three-domain TLRs, with the exception of TLR5, which is under strong positive selection in birds (table 1 and supplementary material S1: fig. S7, Supplementary Material online). Since some avian species have lost a functional TLR5 through pseudogenization (Bainova et al. 2014; and this study, fig. 3), we hypothesize that TLR5 has a specific role in avian immunity when compared, for example, to mammals. The variation in selection acting at TLR5 in different taxa might, for example, reflect the differences in selection against overactivation of gut immunity with flagellated symbiotic microbiota (Iqbal et al. 2005).

Being responsible for direct and specific recognition of structurally heterogeneous PAMPs (Reddick and Alto 2014), the pathogen-mediated selective pressures are particularly diversified in the TLR ECDs. As previously shown in general for vertebrates (Mikami et al. 2012), we also find the majority of TLR PSS in birds to be situated in ligand-binding ECDs (supplementary material S1: fig. S7, Supplementary Material online; the exact positions of PSS are visualized in supplementary material S1: fig. S9, Supplementary Material online and listed in supplementary material S6: table S27, Supplementary Material online). We may predict stronger functional effects of PSS with nonconservative amino acid substitutions that change the physicochemical properties of the particular residues. The highest numbers of these nonconservative PSS are in TLR4 and TLR2A (fig. 4; see also supplementary material S1: table S12 and supplementary materials S7 and S8, Supplementary Material online). Especially in TLR4, most PSS with dramatic changes in charge and hydrophobicity fall into the known functionally important sites of the receptor, with high probability of changing the TLR4 ligand-binding properties. On the other hand, most substitutions in PSS in TLR15, TLR3, and TLR10 [TLR1A] do not significantly change the amino acid physicochemical properties. Despite avian TLR1 [TLR1B] being significantly shorter than TLR10 [TLR1A], we found more PSS in TLR1 [TLR1B] and these were less conservative in their charge than in TLR10 [TLR1A], especially in the predicted ligand-binding region that is spanning the region avoiding gene conversion (supplementary material S1: text S2 and supplementary material S7, Supplementary Material online). This, together with high TLR2A PSS variation, can be possibly explained by distinct binding capacities of the heterodimers formed by the duplicated TLR1/TLR2 subfamily members, where the combination TLR1 [TLR1B]/TLR2A (unlike any other combination) is able to recognize peptidoglycans and efficiently recognizes diacylated bacterial lipoproteins (BLP; Higuchi et al. 2008). Surprisingly, despite its high number of PSS, TLR5 is the most PSS-conservative TLR across all avian taxa (supplementary material S1: table S12, Supplementary Material online). Given that flagellin is the only known TLR5 ligand (Hayashi et al. 2001), this may indicate functional constraints in those species with preserved functional TLR5, limiting any adaptations to larger numbers of relatively minor changes. The PSS identified in this study are consistent with PSS identified by other studies and other taxa especially in TLR15, TLR4, and TLR5 (PSS listed in supplementary material S6: table S27, Supplementary Material online and those located in ECD visualized in fig. 5; visualization of ECDs along with intracellular domains, ICDs, and transmembrane, TM, regions is given in supplementary material S9, Supplementary Material online). Interestingly, from the five amino acid substitutions related to Salmonella enterica resistance variation in chickens (Leveque et al. 2003) only position E301 is also positively selected within birds in general, where negatively charged E/D are substituted with positively charged K/R (fig. 4). The relevance of this position was also supported in other studies of birds (Grueber et al. 2014) and mammals (Vinkler et al. 2009; Wlasiuk and Nachman 2010) and partially with the results of numerous human studies (Arbour et al. 2000; Zareparsi et al. 2005; Shen et al. 2010; Cario 2013; Rupasree et al. 2015; but see also Ohto, Yamakawa, et al. 2012) showing potential importance of variation at a neighboring position D299G (D303 in chicken sequence).

. 4.

Physicochemical properties of the positively selected sites (PSS). All PSS are shown in all avian TLRs—amino acid substitutions are colored according to their physicochemical properties: acidic in red, basic in blue, neutral in purple, polar in green, and hydrophobic in black. The size of an letter corresponds to the procentual proportion of that particular amino acid within the sequence alignment. The numbering is adopted from reference chicken TLRs (for NCBI IDs see supplementary material S1: table S14, Supplementary Material online). PSS which correspond to functionaly important residues (black dot—ligand binding; gray dot—dimerization) are highlighted in bold and orange (identical site) or yellow (topological proximity closer than 5 Å from a functionaly important residue). Ectodomain (ECD), intracellular domain (ICD), and transmembrane (TM) region are visualized; nonconverted region in TLR1/2 is highlighted in pink; cleavage site in TLR3, TLR7, and TLR15 is indicated by a black line tipped with arrows (supplementary material S1: table S11 and text S1, Supplementary Material online). Nonconservative PSS are marked by stars (supplementary material S1: table S12, Supplementary Material online).

. 5.

Positively selected sites and functionaly important sites visualized on 3D extracellular domain structures of avian TLRs. PSS detected in this study are shown in blue. By orange coloration are highlighted the PSS identified at sites with previously described function. Other previously reported functionaly important residues are highlighted in black (ligand binding residues) or in gray (dimerization residues). The total numbers of PSS for each TLR are shown in blue rectangles, where in parentheses are numbers of PSS detected also in other avian/mammalian/both studies (for references see supplementary material S6: table S27, Supplementary Material online). The numbers of PSS with previously described function are provided in orange rectangles, where the number of PSS in topological proximity <5 Å is shown in parentheses. For more detailed information (including ICD and TM), see supplementary material S9, Supplementary Material online.

We found the highest agreement between the PSS and known functional site distribution in TLR5, TLR4, TLR2A, and TLR1 [TLR1B] (23, 18, 18, and 16 positions, respectively; for details see supplementary material S1: table S13 and fig. S10, Supplementary Material online). We show that the PSS frequently evolve toward an amino acid of similar physicochemical properties being gained in distantly related avian taxa, while distinct properties can be found in closely related taxa (supplementary material S1: text S2 and supplementary material S8, Supplementary Material online). This convergent evolution may result from analogous selective pressures induced by partially similar microbial communities in different avian taxa (Waite and Taylor 2014). For example, there are several PSS in avian TLR4 (D345, E375, and G379) that are good candidates for functionally important sites evolving under convergence either toward negative or toward positive charge (supplementary material S8, Supplementary Material online). Being either lipid IVa-binding sites (Meng et al. 2010; Ohto, Fukase, et al. 2012; Scior et al. 2013) or LPS-binding sites (Park et al. 2009; Garate and Oostenbrink 2013; Paramo et al. 2013; Scior et al. 2013), these positions have been identified as positively selected also in several other studies in birds (Grueber et al. 2014) and mammals (Nakajima et al. 2008; Vinkler et al. 2009; Wlasiuk and Nachman 2010; Areal et al. 2011; Shen et al. 2012; Fornůsková et al. 2013).

The lack of agreement between the identified PSS and previously described TLR functional sites in other TLRs may be caused either by the lack of functional studies in some TLRs (e.g., TLR15 and TLR21) or by interspecific variation in TLR-ligand binding (e.g., TLR1, TLR3, and TLR7, where the structural and experimental studies were performed in mammals and not in birds; see supplementary material S1: text S2, Supplementary Material online). In both cases the analysis of PSS may indicate the target sites that are candidates for a functional role and may be the focus for further research. Despite the low number of PSS among the previously reported functional sites of TLR7, which may result from the constraints of conservative ligand binding, two nonconservative PSS located on the TLR7 Z-loop at the proximity of the TLR7 cleavage site (supplementary material S1: text S2, Supplementary Material online; Ewald et al. 2008; Kanno et al. 2013) might affect formation of secondary ligand-binding sites (Zhang et al. 2016).

Although most of the PSS functional sites are located in the ECD and involved in direct ligand binding, in TLR3 we found two PSS that are lying in the transmembrane region and are located in close proximity to the known receptor dimerization residues (fig. 4). Since TLR3 dimerization depends on the dsRNA length (Pirher et al. 2008) and the number of TLR3 molecules involved in the interaction (TLR3 dimer or trimer; Mineev et al. 2014), selection on these residues may play a crucial role in ligand recognition in TLR3. In general, however, most PSS in avian TLR3 are very conservative (fig. 4; for details see also supplementary material S1: table S12 and supplementary material S7, Supplementary Material online), which is consistent with the results reported by Wang, Zhang, Liu, et al. (2016) who identified the TLR3 family as the most conservative TLR family within vertebrates. The avian PSS which can be suggested as candidates for relevant functional effect based on the nonconservativeness of the changes in their physicochemical properties are described in supplementary material S1: text S2, Supplementary Material online.

Conclusions

To our knowledge, this is the first study investigating adaptive evolution in all members of the avian TLR protein family in high number of species and most avian orders. Although we were unable to extract data to all TLR genes from all avian whole-genome assemblies investigated, our study covers 87% of the genes in these species, providing the most comprehensive analysis published so far. Our results provide important novel insights into the history of TLR family evolution including the processes of gene duplication, pseudogenization and diversifying selection. Most importantly, we show that the main duplication event for TLR1 arose before avian and mammalian lineages split into separate clades (giving rise to TLR1 [avian TLR1B] and TLR10 [avian TLR1A]), while avian gene duplication of TLR2 to TLR2A and TLR2B may have occurred in parallel to TLR2 duplication in some mammals that resulted in pseudogenization of the second copy of TLR2. The latter conclusion however remains uncertain due to the power limitation imposed by the short sequences (not affected by gene conversion) aligned in our analysis. Furthermore, we confirmed the unique avian gene duplication in TLR7 based on qPCR copy-number variation analysis. Surprisingly, this recent TLR7 gene duplication in birds emerged independently in many species and representing several clades (Passeriformes, Charadriiformes, Cuculiformes, and Mesiornithiformes). Similarly, the TLR5 pseudogenization, previously reported in several passerines and parrots, is seemingly more widespread in birds. Independent loss of functional TLR5 occurred also in seriemas (Cariamiformes), trogons (Trogoniformes), tropicbirds (Phaethontiformes), sunbitterns (Eurypygiformes), and swifts (Apodiformes).

The results of our analysis of positive selection in avian TLRs allowed us to predict functionally important, interspecifically variable positions. Consistent with some previous findings, these sites were located mainly at the ligand binding extracellular domain, mostly in TLRs where the extracellular domain is exposed to the cell surface and binds structurally diverse ligands (mainly in TLR1 [TLR1B], TLR2A, TLR2B, TLR4, and TLR5). The level of positive selection was low in endosomal TLRs which bind structurally more conservative oligonucleotide ligands (TLR3, TLR7, and TLR21). Positive selection also appears to act more on three-domain than on single-domain TLRs (except for TLR5 that is under very strong positive selection in birds). We also compared the PSS detected in this study with those previously described on an interspecific level in birds and mammals and with all sites previously identified as functionally relevant in TLRs to show strongest evidence for functional effect of the predicted PSS in TLR1 [TLR1B], TLR2A, TLR4, and TLR5. In agreement with the hypothesis that the variation in the TLR ligands is a driving force for selection in TLRs, nonconservative substitutions with the potential to dramatically change surface physicochemical properties were found mainly in two TLRs (TLR2A and TLR4) with highly variable ligands known in mammals. Altogether, our results indicate strong positive selection driving TLR evolution in birds. To better understand the significance of this adaptive evolution in avian TLR genes, further research would benefit from both in silico structural modeling allowing closer prediction of the positively selected physicochemical changes and functional testing of the variant effects in model in vitro systems.

Materials and Methods

Data Set

To create our data set, we used sequence data publicly available for up to 63 avian species per gene representing all orders of Neoaves (supplementary material S2: table S23, Supplementary Material online). Most of the TLR sequences were extracted from the whole-genome data generated by the Avian Phylogenomics Consortium (http://avian.genomics.cn/en/; Zhang, Jarvis, et al. 2014), genomes included in the B10K Project (Avianbase; http://b10k.genomics.cn/; Eöry et al. 2015). The nucleotide sequences of all avian TLRs were obtained by blasting (BLAST v 2.2.25+, NCBI, Zhang et al. 2000; blastn settings: E = 0.1, hits with the highest score and lowest E were accepted) the reference TLR sequences of chicken and zebra finch (for sequence ID number see supplementary material S1: table S14, Supplementary Material online) against the CDS database of all avian whole-genome sequences of the species included in the first B10K study (Zhang, Jarvis, et al. 2014). In cases where the blast search in the CDS databases was not successful, we ran blat search against the genomic scaffolds (e.g., the case of TLR5 pseudogenes; see supplementary material S1: text S3 and table S15, Supplementary Material online; settings: tile size = 11, step size = 5, allowing one mismatch in tile, the number of repetitions of a tile = 1,000,000). In several avian species, where TLRs were previously sequenced, the sequences were found using on-line Web BLAST search (blastn searching cds database with default settings; https://blast.ncbi.nlm.nih.gov/Blast.cgi; accessed in March 2013) and downloaded directly from the GenBank nucleotide database (Benson et al. 2012). The list of NCBI accession numbers of the sequences used is given in supplementary material S2: table S23, Supplementary Material online. Basic description of the TLR CDS analyzed is provided in supplementary material S1: text S4, tables S16 and S17 and more detailed in supplementary material S2: table S25, Supplementary Material online.

Multiple Sequence Alignment

First, all CDSs were translated to extract the protein sequences of all representatives of each TLR family member. Then, the multiple amino acid sequence alignment was done for each TLR separately using the ClustalW alignment tool implemented in Geneious v. 9.0.4 (Biomatters Limited; Kearse et al. 2012), or using MAFFT v. 6.850 (Katoh and Standley 2013) in the case of the duplicated genes from the vertebrate TLR1 family. For the construction of nucleotide codon alignments corresponding to the protein sequence alignments, we used the PAL2NAL webtool v. 14 (http://www.bork.embl.de/pal2nal/; accessed in July 2014; Suyama et al. 2006). All created multiple alignments are shown in supplementary material S10, Supplementary Material online.

Gene Conversion Analysis

The occurrence of gene conversion was statistically tested on nucleotide sequence alignment of all duplicated genes (TLR1, TLR2, and TLR7) using GENECONV v. 1.81 (S. A. Sawyer, Washington University in St. Louis, https://www.math.wustl.edu/∼sawyer/geneconv/). Full program setting is provided for all investigated genes in supplementary material S3, Supplementary Material online.

Phylogenetic Tree Analyses

Phylogenetic analysis was done only in the duplicated avian TLRs (i.e., in TLR1, TLR2, and TLR7), also including other vertebrate taxa, that is, mammals, reptiles, amphibians, and fishes (involved species and their sequences IDs are listed in supplementary material S2: table S24, Supplementary Material online). For the TLR1 family, the phylogenetic analysis was based on the sequence alignment of nonconverted amino acid regions only. Amino acid sequences were used instead of nucleotide sequences to avoid the biasing effect of rapid evolution on the third codon positions on broad evolutionary scales (i.e., higher vertebrate taxa). For the phylogenetic analysis, we used two approaches: maximum likelihood calculated in PhyML v. 3.0 software package (Bootstrap: 1,000; Tree type: SPR&NNI; Substitution model: LG; Guindon et al. 2010) and Bayesian estimation of phylogeny calculated using MrBayes v. 3.2.1 (number of generations: 1,000; burn-in fraction was set to default of 25%; the critical value for topological convergence diagnostic was 0.01; the best suited aa model for TLR1: Fixed(Jones) and for TLR2: Fixed(Wag); Ronquist et al. 2012). Schematic phylogenetic trees for visualization of TLR7 duplication (fig. 2) and TLR5 pseudogenization (fig. 3) were constructed using previously published avian phylogeny (Jarvis et al. 2014). All trees were graphically adjusted in FigTree v1.3.1 (A. Rambaut, University of Edinburgh, United Kingdom; http://tree.bio.ed.ac.uk/software/figtree/).

Copy Number Variation Analysis

To verify the increased number of TLR7 copies in selected passerine species, we performed a qPCR copy number variation (CNV) analysis. Tissue samples of four model species (Ficedula albicollis, Pseudopodoces humilis, Taeniopygia guttata, and Zonotrichia albicollis) used for this analysis were obtained from various genetic banks (for details see supplementary material S1: table S18, Supplementary Material online). DNA was extracted using the DNeasy Blood and Tissue Kit (QIAGEN) and stored at −20°C. As reference single copy genes, we used TLR3 and TLR4. Primers were designed to amplify a conserved region of a similar length for all studied genes (for details see supplementary material S1: table S19, Supplementary Material online). The specificity of the primers was previously verified by Sanger sequencing (Applied Biosystems 3130xl Genetic Analyzer) of a broader surrounding region of each gene in all investigated species (for PCR conditions see supplementary material S1: table S20, Supplementary Material online; NCBI accession numbers are listed in supplementary material S1: table S21, Supplementary Material online). Each sample was run in triplicate in LightCycler 480 Instrument II (Roche) using EvaGreen Dye (Biotum; Mao et al. 2007). The qPCR efficiency was calculated based on a dilution series (5-times dilution) for each gene and sample in LightCycler 480 software v1.5.1 using both 2nd Derivate function and automatic Fit Point method (the values measured are shown in supplementary material S1: table S4, Supplementary Material online). The TLR7 copy numbers were then calculated based on a modified version of the formula proposed by Pfaffl (2001; eq. 1, where Eff stands for the PCR efficiency and Cp stands for the crossing point).

Positive Selection Analysis and Estimates of Conservativeness of Amino Acid Substitutions

For the detection of selective pressures acting on each avian TLR, the codon alignment generated by PAL2NAL tool (http://www.bork.embl.de/pal2nal/; last accessed July 2014; Suyama et al. 2006) was used to ensure correct alignment of codons. All regions involving gaps were removed before the analysis and tested separately. The problem of missing sequence data did not impact our estimates since the number of sequences sampled was not too small in any of the TLRs and thus the dN/dS ratios could have been estimated. The position numbering of PSS followed the chicken reference sequence (for NCBI IDs see supplementary material S1: table S14, Supplementary Material online). To test for positive selection acting on individual residues at the interspecific level in avian TLRs, we used two methods based on the hierarchical Bayes (Bayes Empirical Bayes, BEB) approach implementing the Markov chain Monte Carlo routine—PAML (Phylogenetic Analysis by Maximum Likelihood; Yang 2007) and FUBAR (A Fast, Unconstrained Bayesian AppRoximation for Inferring Selection; Murrell et al. 2013). Being based on unrelated preconditions, these two tests provide independent estimates of positive selection. In PAML (Version 4.7), the codon-based substitution models (codeml) using comparison of neutral M8a (beta&ω = 1) with alternative M8 (beta&ω) model were adopted. The likelihood ratio test (LRT) for comparison of two nested models was calculated using the chi-square approximation: χ2 = 2×(lnLM8−lnLM8a), where LM8 and LM8a are the likelihood values. The degrees of freedom (df) were established as the difference in the numbers of parameters in used models (for details see supplementary material S1: table S10, Supplementary Material online). If the LRT is significant (≤0.05), positive selection is considered to be detected. The BEB approach (Yang 2007) was then used to determine site specific posterior probabilities indicating positive selection (≥0.9) at specific codons. The phylogenetic tree needed for PAML analysis was constructed based on the previously published avian phylogeny (Jarvis et al. 2014). FUBAR (Fast Unconstrained Bayesian AppRoximation) analysis was performed at the Datamonkey server (http://www.datamonkey.org/; accessed in July 2014; Pond and Frost 2005) with the significance level of posterior probability established by default to 0.9. The FUBAR algorithm was used because it is more robust and much faster than other available selection tests, which are based on the random effect likelihood (REL; Murrell et al. 2013) methods. The degree of dissimilarity in biochemical properties of amino acid substitutions was tested using the PRIME tool (PRoperty Informed Model of Evolution; accessed in July 2014) available at the Datamonkey server (Pond and Frost 2005). This tool builds on the same conceptual frameworks as MEME (Murrell et al. 2012), but allows the nonsynonymous substitution rate β to depend not only on the site in question but also on what type of residues are being exchanged. Both predefined sets of five amino-acid composite properties were used for PRIME analysis: 1) Polarity index, Secondary structure factor, Volume, Refractivity/Heat Capacity, and Charge/Iso-electric point (Atchley et al. 2005) and 2) Chemical Composition, Polarity, Volume, Iso-electric point, and Hydropathy (Conant et al. 2007) on the significant level of posterior probabilities ≥0.9. Amino acid physicochemical properties (chemistry, charge, and hydrophobicity) at all PSS were graphically visualized using a web-based application Weblogo v. 3.5 (http://weblogo.threeplusone.com/create.cgi; last accessed October 2017; Crooks et al. 2004).

We compared the PSS identified in this study with the results of other studies that focused on avian and mammalian TLR evolution identifying positive selection at the interspecific level (Nakajima et al. 2008; Vinkler et al. 2009; Wlasiuk et al. 2009; Wlasiuk and Nachman 2010; Alcaide and Edwards 2011; Areal et al. 2011; Huang et al. 2011; Shen et al. 2012; Fornůsková et al. 2013; Grueber et al. 2014; Vinkler et al. 2014; Wang, Zhang, Chang, et al. 2016) and studies that described the functionally relevant residues for TLR1 (Jin et al. 2007; Omueti et al. 2007), TLR2 (Underhill et al. 1999; Lorenz et al. 2000; Xu et al. 2000; Kang and Chae 2001; Tao et al. 2002; Gautam et al. 2006; Kang et al. 2009), TLR3 (Sarkar et al. 2003; Choe et al. 2005; de Bouteiller et al. 2005; Bell et al. 2006; Sun et al. 2006; Ranjith-Kumar et al. 2007; Liu et al. 2008; Pirher et al. 2008; Luo et al. 2012; Mineev et al. 2014), TLR4 (Poltorak et al. 1998; Ronni et al. 2003; Nishitani et al. 2006; Kim et al. 2007; Walsh et al. 2008; Park et al. 2009; Resman et al. 2009; Meng et al. 2010; Ohto, Fukase, et al. 2012; Ohto, Yamakawa, et al. 2012; Garate and Oostenbrink 2013; Paramo et al. 2013; Scior et al. 2013; Wang, Su, et al. 2016), TLR5 (Jacchieri et al. 2003; Andersen-Nissen et al. 2007; Keestra et al. 2008; Yoon et al. 2012; Ivicak-Kocjan et al. 2013; Song et al. 2017), TLR7 (Wei et al. 2009; Yu et al. 2013; Tseng et al. 2014; Gentile et al. 2015; Gupta et al. 2016; Zhang et al. 2016), TLR10 (Hasan et al. 2005; Nyman et al. 2008; Guan et al. 2010; Jang and Park 2014), TLR15 (Wang, Zhang, Chang, et al. 2016), and TLR21 (Keestra et al. 2010).

Protein Structure Modeling

To predict 3D structures of all avian TLRs, the I-TASSER v. 5.0 server (https://zhanglab.ccmb.med.umich.edu/I-TASSER/; last accessed November 2016; Roy et al. 2010) was used. Since all TLRs are transmembrane proteins, the ECD and ICD (and in the case of TLR3 also TM region) were modeled separately, always based on the chicken reference sequence (for NCBI IDs see the supplementary material S1: table S14, Supplementary Material online). The domains were identified by SMART v. 7.0 web tool (http://smart.embl-heidelberg.de/; last accessed November 2016; Letunic and Bork 2018), the amino acid ranges of ECD and ICD for each TLR are provided in supplementary material S1: table S22, Supplementary Material online. The I-TASSER model with the highest C value reflecting the confidence score for estimating the quality of predicted models was always downloaded and used for further analysis. The graphical visualization of important residues was then done using the PyMOL Molecular Graphics System (Version 1.7.6, Schrödinger, LLC).

Assessing Function of Positively Selected Sites

The list of the previously reported functionally important positions as well as PSS detected in other interspecific studies for all TLRs were obtained by a detailed review of the published literature (for the complete list of references see supplementary material S6: table S27, Supplementary Material online). The distance of any PSS detected in our study from these previously described functionally important residues was measured on the 3D structural models of TLRs (obtained in previous step) in the PyMOL Molecular Graphics System (Version 1.7.6, Schrödinger, LLC) using python command iterate and plugin distancetoatom, where PSS lying in distance closer to 5 Å were considered as closely connected to the functionally important residues, that is, having potential influence on the receptor function.

Click here for additional data file.