Adaptive evolution of virus-sensing toll-like receptor 8 in bats.

Recently, bats have gained attention as potential reservoir hosts for emerging zoonotic single-stranded (ssRNA) viruses that may prove fatal for humans and other mammals. It has been hypothesized that some features of their innate immune system may enable bats to trigger an efficient early immune response. Toll-like receptors (TLRs) represent a first line defense within the innate immune system and lie directly at the host-pathogen interface in targeting specific microbe-molecular patterns. However, the direction and strength of selection acting on TLRs are largely unknown for bats. Here, we studied the selection on viral ssRNA sensing TLR8 based on sequence data of 21 bat species. The major part (63 %) of the TLR8 gene evolved under purifying selection, likely due to functional constraints. We also found evidence for persistent positive selection acting on specific amino acid sites (7 %), especially when compared to viral TLR evolution of other mammals. All of these putatively positively selected codons were located in the ligand-binding ectodomain, some coincidenced or were in close proximity to functional sites, as suggested by the crystallographic structure of the human TLR8. This might contribute to the inter-species variation in the ability to recognize molecular patterns of viruses. TLR8 evolution within bats revealed that branches leading to ancestral and recent lineages evolved under episodic positive selection, indicating selective selection pressures in restricted bat lineages. Altogether, we found that the TLR8 displays extensive sequence variation within bats and that unique features separate them from humans and other mammals.

Introduction

Zoonotic viruses that emerge from wildlife pose a serious threat for humans and domestic livestock. Over the past years, bats have attracted growing attention as a potential reservoir of emerging zoonotic viruses that may have a high impact on human health such as the severe acute respiratory syndrome-like (SARS) and Middle East respiratory syndrome-like (MERS) coronaviruses, Ebola and Marburg filoviruses, Hendra, Nipah, and other newly recognized paramyxoviruses, hepaciviruses, pegiviruses, and hantavirus, as well as new lineage of influenza A virus (reviewed in Calisher 2015). In addition, bats have long been known as vectors for lyssaviruses causing rabies. Apparently, all aforementioned viruses are single-stranded RNA (ssRNA) viruses and, except for lyssaviruses, cause little or no pathology in bats (e.g., Wynne and Wang 2013; Brook and Dobson 2015). This raises the question why certain viruses can infect and persist in apparently healthy bats whereas they may cause severe symptoms and even death in other vertebrate species (Wang et al. 2011; Bean et al. 2013; Plowright et al. 2015).

Typically, mechanisms of infection and resistance involve direct interactions of the host and the pathogen at a molecular level. In general, viruses are first recognized by the innate immune system via pattern recognition receptors (PRR), including the Toll-like receptor (TLR) family, which sense microbe-associated molecular patterns (MAMPs; Takeda et al. 2003; Kawai and Akira 2010). Resistance to virus infection can evolve by selection of mutations that alter amino acids in the binding region of receptor molecules. This sort of positively selected residues point to specific molecular recognition interfaces between host and viral components that have adapted and counter adapted in a long series of classical Red Queen conflicts and can be used to infer the history of the host–virus “arms race” during their co-evolution (Daugherty and Malik 2012; Hancks et al. 2015). Thus, PRRs may rapidly evolve to evade pathogen antagonism and to ensure activation of a robust immune response. Although bats appear to share the immune cell and gene repertoire with other mammals, it has been repeatedly hypothesized that the long-term co-existence of bats and viruses must have imposed strong selective pressures on the bat genome and that those genes related to the first line of antiviral defense within the innate immune system, e.g., TLRs, may reflect this co-evolution (Wibbelt et al. 2010; Cowled et al. 2011; Baker et al. 2013; Wynne and Wang 2013; Zhang et al. 2013).

Binding of MAMPs to TLRs triggers a series of events leading to increased expression of proinflammatory genes (Kawai and Akira 2006). Within TLRs, the extracellular part, consisting of multiple leucine rich repeats (LRR), is involved in pathogen recognition and evolves faster than the highly conserved intercellular Toll/interleukin-1 receptor (TIR) domain which is responsible for intracellular signalling (Mikami et al. 2012). Most mammals have 10–12 different TLRs, each recognizing different ligands, hereby nucleic acids of viruses are recognized by endosomal TLRs (3, 7, 8, and 9; Alexopoulou et al. 2001; Lund et al. 2004; Xagorari and Chlichlia 2008). Single nucleotide polymorphisms (SNPs) in all TLRs have been correlated with disease susceptibility demonstrating that single amino acid alterations may have an effect on immune response mechanisms (humans: Schroder and Schumann 2005; Wong et al. 2010; Taylor et al. 2012; Uciechowski et al. 2011; Wang et al. 2014; rodents: Tschirren et al. 2013).

Studies of the evolutionary dynamics of TLRs have revealed that selective pressures vary between TLRs and taxa. In general, a high degree of purifying selection is found for all TLRs and is particularly reported for nucleic acid-sensing TLRs due to functional constraints (Roach et al. 2005; Barreiro et al. 2009; Mukherjee et al. 2009). However, an increasing number of studies detected specific residues within the LRR binding region to be under positive selection in almost all TLRs (mammals: Areal et al. 2011; apes: Wlasiuk and Nachman 2010; rodents: Tschirren et al. 2011; bats: Zhang et al. 2013, Escalera-Zamudio et al. 2015; birds: Alcaide and Edward 2011; Grueber et al. 2012; amphibians: Babik et al. 2014) indicating that the ligand-binding domain is subject to classical Red Queen conflicts as mentioned above (Daugherty and Malik 2012; Hancks et al. 2015). Indeed, a recent study highlighted that nucleic acid sensing TLRs in bats exhibited unique mutations fixed in ligand-binding sites (Escalera-Zamudio et al. 2015). They found further evidence for episodic positive selection acting specifically upon the bat lineage (based on eight bat species belonging to three families) suggesting potential functional differences between bat and other mammalian TLRs. Further, several ancestor genes involved in innate immunity (e.g., TLR7, NLRP3, and MAP3K7) were found to be under stronger positive selection in the genomes of the Australian black flying fox (Pteropus alecto) and David’s Myotis (Myotis davidii) than their orthologues in seven other mammalian species (Zhang et al. 2013). Interestingly, TLR7, which is besides TLR8 and RIG-1 (an intracellular PRR of the retinoic acid-inducible gene I family, Jensen and Thomsen 2012) responsible for the detection of ssRNA viruses, was one of those being under positive selection. This is remarkable because TLR7 has been suggested to be under strong purifying selection in other mammalian taxa (Barreiro et al. 2009; Georgel et al. 2009; but see Areal et al. 2011). Here, we investigated traces of selection acting upon the viral ssRNA sensing TLR8 gene to shed further light on evolutionary dynamics of the viral sensing TLRs within different bat orders and their potential implications in ligand-binding affinities. In addition to the aforementioned studies, we expanded the number of investigated bat taxa to 21 bat species, belonging to seven families, to increase resolution. We included three species of the Rhinolophoidae, which are known as a main reservoir for beta-coronaviruses. Specifically, we studied the molecular variation, the strength and direction of selection, focussing on the presence of putatively positively selected amino acids in the LRR-binding region. Finally, we discuss the altered amino acid residues and their potential implications in ligand-binding affinities in detail and in comparison to new data on two ssRNA ligand-binding sites of the human TLR8 molecule (Tanji et al. 2015).

Materials and methods

Tissue collection, DNA and RNA extraction

Tissue samples from 15 species where obtained both from free-ranging bats and also from museum specimen (Table S1, Supplementary Material). Wild bats were captured at the following locations: Bulgaria (Rhinolophus euryale, R. hipposiderus, R. ferrumequinum, and Myotis myotis, permission number permit of the Ministerstvo na Okolnata Sreda i Vodita, Sofia and RIOSV Ruse; No. 57/18.04.2006, No. 100/04.07.2007), Costa Rica (Saccopteryx bilineata, Pteronotus parnelli, Carollia perspicilatta, Vampyressa pusilla, Vampyressa bidens, Ministerio del Ambiente y Energia; # 036-2009-SINAC: No. 137–2010–SINAC and 168–2011–SINAC) and from Panama (Noctilio albiventris, Noctilio leporinus, Ministerio de Agricultura, Ganaderia, Acuacultura y Pesca #04508). We captured wild bats using mist-nets or harptraps, depending on the specific efficacy, when bats emerged from their roosts. A 4-mm wing-tissue sample was collected using a sterile biopsy punch and the sample immediately stored in 96 % ethanol until DNA isolation. Liver samples (stored in RNAlater at −20 °C) of M. myotis and N. albiventris were used for RNA analyses. Animals were euthanized in accordance with appropriate guidelines (Gannon and Sikes 2007) and under the license of national authorities. In addition, we received frozen tissue material from a museum specimen from Eidolon helvum (Natural History Museum Berlin). In Germany, we collected tissue material from four bat species that were obtained from central depositories where bats killed at wind turbines are deposited. Samples were stored at −80 °C after collection (Nyctalus noctula, Nyctalus leisleri, Vespertilio murinus, and Pipistrellus nathusii, Landesumweltamt Brandenburg, Vogelschutzwarte Buckow). Total DNA was extracted from the biopsies using standard protocols (DNA extraction Kit, Genial). RNA of M. myotis and N. albiventris was extracted using RNA extraction-kit (RNeasyKit, Quiagen). Possibly remaining genomic DNA was removed by digestion and the cDNA-library established by reverse transcriptase (RevertAid H Minus First cDNA Synthesis Kit, Fermentas).

Available TLR8 sequences of other bat species were included in our analyses (Pterpous alecto GenBank accession number XM_006906028, Eptesicus fuscus XM_008156576, Myotis davidii XM_006763795, M. brandtii XM_005880947, M. lucifugus XM_006088606) and compared to the equine (XM_014728737) and human (XM_005274543) orthologous.

Primer design

We first aligned the published bat TLR8 sequences from the Australian black flying fox (P. alecto) and three species of mouse-eared bats (M. davidii, M. brandtii, and M. lucifugus) to identify highly conserved regions. We then designed a diverse array of primer sets to amplify the whole TLR8 coding region (3,056 bp without primers out of 3,123 bp of the human TLR8-coding sequence) in all selected bat species (Supplementary Table S2). Established primers spaced overlapping fragments between 700 and 800 bp in length.

TLR8 amplification and sequencing

We amplified the TLR8 coding region in two to three overlapping amplicons using a proofreading, long-range polymerase (Myfi-Polymerase, Thermoscintific). PCR reactions were conducted in a total volume of 25 μL consisting of 10–50 ng genomic DNA, 10 μmol/L of primers and 5 × MyFi reaction buffer with 1 mmol/L dNTPs, 3 mmol/L MgCl2, and 2 U/μL MyFi DNA Polymerase. PCR conditions were as follows: initial denaturation at 95 °C for 3 min followed by 36 cycles of denaturation at 95 °C for 15 s, annealing between 50 and 58 °C for 15 s and elongation at 72 °C for 4–5 min depending on fragment length and a final elongation step by 72 °C for another 10 min. PCR products were electrophoresed in a 1 % agarose gel. Positive amplicons were purified (Exonuclease I, FastAP thermosensitive alkaline phosphatase, Thermoscientific) and bi-directly sequenced using Big-Dye Terminator technology and an ABI PRISM 310 Automated Genetic Analyser (Applied Biosystems, Foster City CA, USA). Bat TLR8 sequences have been deposited in GenBank under the following accession numbers KU163591-KU163606.

Statistical analysis

Nucleotide sequences were edited and aligned manually using Chromas Lite 2.1.1 (Technelysium Pty Ltd 1998–2012) and MEGA 6.0 (Tamura et al. 2013). Sequence percent identities compared to other mammalian species (human and horse) were obtained by homology Blast search (http://blast.ncbi.nlm.nih.gov) and calculated by MEGA. To estimate evolutionary divergence of the TLR8 gene sequences within bats, we calculated besides sequence percent identities the nucleotide diversity π by estimating the mean number of base differences per site using a bootstrap procedure with 500 replicates with Mega 6.0.

Under neutrality, coding sequences are expected to present a ratio (ω) of non-synonymous (dN) over synonymous (dS) base pair substitutions that does not significantly deviate from 1 (ω = dN/dS = 1), while significant deviations may be interpreted as either the result of positive (ω > 1) or purifying (ω < 1) selection (Ellegren 2008). We estimated the relative rates of dN and dS according to the Nei and Gojobori (1986) method, applying the Jukes-Cantor correction for multiple hits in 500 replicates (Jukes and Cantor 1969) implemented in MEGA. Evidence for persistent positive selection in the past was assessed through the program CODEML (included in the software package Paml 3.15, Yang 2007), which estimates ω among sites applying different models of codon evolution (neutral, purifying, or positive). First we tested whether ω differs among sites by comparing model M0, which assumes a constant ω across all sites to model M3 (discrete), which allows ω to vary among sites, estimating the proportion of conserved, neutral and positive selected codons from the data. To test for the presence of sites evolving under positive selection we used the nested model pair M7, allowing codons only to evolve neutrally or under purifying selection and M8, which adds a class of sites under positive selection assuming a beta distribution for ω. The two nested models pairs were compared using the likelihood ratio test (LRT, α = 95 %) with two degrees of freedom (Yang et al. 2005). For a better comparison to number and position of positively selected sites detected in TLRs of other taxa we applied a set of methods to detect candidate-codons under persistent positive selection. Based on the methodology adopted by other authors and in previous studies (Wlasiuk and Nachman 2010; Areal et al. 2011; Matos et al. 2013; Castel et al. 2014; Darfour-Oduro et al. 2015), only codons identified by three or more methods were considered to be putatively under positive selection. It has been suggested that the application of several selection assessment methods, followed by the classification of sites detected in consensus as positively selected can minimize false positives (e.g., Kosakovsky Pond and Frost 2005, Maldonado et al. 2014). We applied Model M8 implemented in the PAML package and the SLAC (single likelihood ancestor counting), FEL (fixed-effect likelihood), REL (random effect likelihood), and FUBAR (Fast, Unconstrained AppRoximation analyses) methods implemented in the Hyphypackage (www.datamomkey.org; Kosakovsky Pond et al. 2005; Delport et al. 2010; Murell et al. 2013). Codons with p values <0.1 SLAC and FEL, with Bayes Factor >50 for REL, and with posterior probabilities >0.9 for M8 and FUBAR were considered as candidates under positive selection. We further tested the degree of dissimilarity of amino acid substitutions according to their physiochemical properties using PRIME (Property Informed Model of Evolution), with the amino acid property set according to Atchley et al. (2005), also available in the Datamonkey web server. A change on this properties was considered significant if the posterior probabilities exceeded >0.9. The evolutionary conservation of amino acid positions was then predicted using the ConSurf tool (Ashkenazy et al. 2010), with the assumption that positively selected residues were least conserved.

Finally, to test for episodic positive selection Branch-Site REL and MEME (mixed effects model of evolution) implemented in the Hyphypackage were used. The Branch-Site REL model estimates proportion of sites under selection along tree branches and allows evolutionary rate to simultaneously vary along phylogenetic branches and sites (Kosakovsky Pond et al. 2011). The MEME method identifies lineage-specific events of positive selection at specific sites, even though the same site is under purifying or neutral selection in other lineages (Murell et al. 2012). Codon sites were considered to be under episodic positive selection at a significant level p < 0.05. A Neighbour Joining tree was used as working topology in all analyses.

The domain architecture for TLR8 was determined by using a Simple Modular Architecture Research Tool (SMART, Schultz et al. 1998), which combines homology searches with predictive algorithms to identify putative protein domains. TM domains were additionally identified by TMPRED (Hofmann and Stoffel 1993), which estimates transmembrane regions on the basis of hydrophobicity of predicted TM-helices. The LRRfinder tool (www.lrrfinder.com, Offord et al. 2009) was used to calculate, in addition to the domain architecture, the surface accessibility (buried or exposed) and the predicted secondary structure (helix, coil or strand) for each amino acid, which we used to discuss the potential functionality of amino acids that have been predicted to have evolved under persistent positive selection. Tertiary protein structure predictions were generated using Swiss Model (Biasini et al. 2014) based on the known crystallised human (hTLR8) structure. Positions of sites of persistent positive selection affecting particular amino acid positions in the protein were visualized by PyMol 1.2.8 (De Lano 2002). Amino acid positions and LRR numbering were assigned according to the human h TLR8 molecule throughout the article (Tanji et al. 2013, 2015).

Results

Domain architecture and diversity

Similar to other mammalian TLRs, TLR8 of all investigated bat species showed the typical domain architecture with a signal peptide followed by a LRR ectodomain for ligand binding, a transmembrane domain (TM), and a cytoplasmatic TIR domain for signal transduction. Number and locations of LRRs within the ectodomain varied substantially among species (15–19 LRR, Supplementary Fig. S1). We could not detect TM domains in the Noctilio and Rhinolophus species by SMART analyses. However, based on TMPRED analyses, we identified the TM domain in the appropriate region for species of these two genera (Supplementary Fig. S2). Expression was confirmed by RT-PCR in N. albiventris and M. myotis, in spleen and liver.

The open-reading frame of the amplified single large TLR8 exon displayed no stop codons. We detected several indels (1–3 amino acids), which always concerned three contiguous nucleotides or a multiple of it, thus frame shifts never occurred (Supplementary Fig. S3). We obtained 1,018 amino acids (aa) out of the 1,024 aa that were previously described for hTLR8. However, the obtained lengths of amino acid sequences varied between 1,015 aa (Vespertilionidae) and 1,020 aa (Phyllostomidae), as a result of the presence of indels. Sequence length variations were all located in the LRR ectodomain.

The nucleotide sequences showed 79–85 and 76–81 % sequence identity to their equine and human orthologues. The adopted TLR8 protein sequences showed less similarity with 74–76 % (equine) and 70–72 % (human) amino acid sequence identity, respectively, with higher similarities in the TM and TIR-domain compared to the LRR ectodomain (Supplementary Table S3a, b, c). Mean nucleotide (nt) and amino acid (aa) similarities within bats reached 84 % (nt) and 80 % (aa), and ranged from 78 % (nt) and 72 % (aa) between M. myotis and R. hipposideros to 99 % (nt and aa) identity between E. helvum and P. alecto, displaying a remarkable intra-order diversity. In all comparisons of the LRR ectodomain, amino acid sequences showed a lower level of similarity between taxa than the corresponding nucleotide sequences. The opposite was true for the TM and TIR domain, were amino acid sequences showed higher percent sequence identities between all taxas than the nucleotide sequences (Supplementary Table S3a, b, c).

Evidence for persistent positive selection

Non-synonymous substitutions did not exceed those of synonymous substitutions when the whole sequence was considered or when only the LRR-region was analyzed with highest values for ω in the LRR ectodomain (Table 1). However, we detected significant ω heterogeneity within the TLR8 sequences (M0 vs M3, p < 0.001). According to model M3 the large majority of codons evolved under purifying selection (63 %, ω = 0.06), 30 % of sites evolved neutrally (ω = 0.83) and 7 % of all sites were detected to evolve under positive selection (ω = 2.34; Table 2). When comparing models that estimate ω with a beta distribution also model M8 allowing for positive selection fitted the data significantly better than model M7 that considered only neutral sites (p < 0.001). To identify specific sites that might have evolved under persistent positive selection, we used different analytical approaches (see method section, Supplementary Table S4). Estimated PSS occurred clustered in certain regions within the amino acid sequence alignment and were all located in the ligand-binding LRR ectodomain (Fig. 1, Supplementary Fig. S3). This clustering was even pronounced in the 3D structure of the molecule (see below). 24 specific codon sites were estimated as positive selected sites (PSS) by at least three model analyses and were thus elected as putatively PSS in bats (Table 3, Supplementary Table S4). Amino acid character alterations were present in all 24 PSS including non-polar and polar as well as basic polar and acidic polar side chain polarities (Table 3). This was underlined by the amino acid conservation scores of the ConSurf analyses, which indicated all PSS to be least conserved (scores >1.5, Table 3).

Table 1

Estimates of average evolutionary divergence within the different domains of the chiropteran TLR8

TLR domain	N	Length bp (aa)	π	dS (±SE)	dN (±SE)	ω
LRR ectodomain	21	2,460 (820)	0.177 ± 0.004	0.451 ± 0.023	0.150 ± 0.007	0.517
TM/TIR	21	612 (204)	0.086 ± 0.006	0.409 ± 0.045	0.025 ± 0.004	0.114
All	21	3,072 (1,024)	0.159 ± 0.004	0.442 ± 0.019	0.112 ± 0.006	0.437

aa amino acid, bp base pair, N number of inferred haplotypes, π number of base differences per site, d number of synonymous substitutions per synonymous site, d number of nonsynonymous substitutions per nonsynonymous site, SE standard error estimates by bootstrap procedure with 500 replicates estimated by MEGA, ω d /d ratio computed by SLAC

Table 2

Results of selection analyses on TLR8 based on codon evolution models by PAML

Model	Parameter estimates	Log likelihood	−2lnΔLa
M0	ω = 0.387	−1,7356.9	1,032.2***
M3b	p 0 = 0.635, ω0 = 0.061 p 1 = 0.299, ω1 = 0.825 p 2 = 0.067, ω2 = 2.347	−1,6840.8
M7	p = 0.178, q = 0.332	−1,6887.6	86.6**
M8	p 0 = 0.916, p = 0.253, q = 0.643 p 2 = 0.084, ω = 2.134	−1,6844.3

***p < 0.0001, **p < 0.001

-2lnΔL is distributed approximately as χ 2 with two degrees of freedom

p proportion of sites were ω < 1, p proportion of sites were ω = 1, p proportion of sites were ω > 1

Fig. 1

Predicted protein domain architecture of the TLR8 gene of Rhinolophus ferrumequinum and putatively positive selected sites in bats. TM domains are shown in dark blue; other domains are labeled accordingly. Protruding loops are marked in orange. Predicted regions involved in ligand binding (red rectangles) and dimerization (circles: orange inactivated, red activated dimer) according to the human TLR8 molecule (Tanji et al. 2013, 2015) are marked below the gen cartoon. Amino acid sites with statistically significant ω values in bats are indicated by triangles (M8: posterior probabilities shaded triangles >95 %, blanc triangles >90 %; FEL, PRIME: p values shaded triangles <0.05, blanc triangles <0.1; REL: Bayes factor shaded triangles >50; SLAC: p values shaded triangles <0.1; FUBAR: posterior probabilities shaded triangles >90 %, blanc triangles >85 %). Putatively positive selected sites (PSS) in bats estimated by more than three models are shown as red x. The four clustered PSS regions are marked by red rectangles. PSS in mammals according to Areal et al. (2011) in Suidae (Darfour-Oduro et al. 2015) and in primates (Wlasiuk and Nachman 2010) are mapped as green triangles

Table 3

Putatively positively selected amino acid sites across the Chiropteran TLR8 gene

Positiona	Domain locationa	Amino acid alterations	Amino acid properties		Surface accessibilityb	Conservation scoresc
			Polarity	Side chain charge
43	LLR-NT	E-S-G-A-V-F	Non-/polar/neutral	Acidic/non-acidic	Exposed	3.269
104	Loop in LRR 2	F-C-W-I-K-S	Non-/polar/neutral	Non-acidic/basic	Buried	3.258
109		R-I-T-K	Non-/polar/neutral	Non-acidic/basic	Exposed	3.293
236	LRR 7	T-I-A-N-S-V	Non-/polar/neutral	Non-acidic	Exposed	2.704
238		E-K-G-R-S	Non-/polar/neutral	Acidic/non-acidic/basic	Exposed	3.302
277	LRR 8	S-Y-Q-V-F-K	Non-/polar/neutral	Non-acidic/basic	Exposed	2.776
312	LRR 10	H-C-Y-D-N-S	Non-/polar/neutral	Acidic/non-acidic	Exposed	3.123
338		S-N-H-V-R	Non-/polar/neutral	Non-acidic/basic	Exposed	2.538
352	Loop in LRR 11	Q-E-D-A-T	Non-/polar/neutral	Acidic/non-acidic	Exposed	2.170
354		P-L-A-V	Non-polar/neutral	Non-acidic	Exposed	2.812
356	LRR 11	Y-H-N-F-R	Non-polar/polar	Non-acidic/basic	Exposed	2.120
368		S-L-R-C-H-F-Y	Non-/polar/neutral	Non-acidic/basic	Exposed	2.447
412	LRR 13	T-S-G-N	Non-/polar/neutral	Non-acidic	Exposed	1.765
420	LRR 14	L-G-V-A	Non-polar/neutral	Non-acidic	Buried	1.544
421		N-K-S-E	Polar/neutral	Acidic/non-acidic/basic	Exposed	2.971
447	Z-loop	V-G-L-I-P	Non-polar/neutral	Non-acidic	Buried	2.496
481		H-K-R-T-Q	Polar/neutral	Non-acidic/basic	Exposed	2.703
608	LRR 19	E-S-R-V-G-K-A-T	Non-/polar/neutral	Acidic/non-acidic/basic	Exposed	3.542
677	LRR 22	K-T-R-C-S-Q-H	Non-/polar/neutral	Non-acidic/basic	exposed	3.348
685		Q-L-Y-R-E	Non-/polar/neutral	Acidic/non-acidic/basic	Exposed	2.411
691	LRR 23	Q-N-T-K-Y	Non-/polar/neutral	Non-acidic/basic	Exposed	1.961
699		H-C-R-Q-S	Non-/polar/neutral	Non-acidic/basic	Exposed	3.195
712		K-S-H-P	Non-/polar/neutral	Non-acidic/basic	Exposed	2.193
723	LRR 24	R-K-I-L-H-S	Non-/polar/neutral	Non-acidic/basic	Exposed	2.853

Relative to the human TLR8 protein

According to the LRR-finder-tool

According to ConSurf analyses (a high value indicates less conservation)

The eight PSS detected by Escalera-Zamudio et al. (according to their M8 analyses) matched with the estimated PSS of our M8 model analyses, except position 506, which was identified to have undergone episodic positive selection by our analyses (see below and Supplementary Fig. S3). Four PSS (position 236, 338, 481, and 712) matched with predicted positive selected codons in mammals (Areal et al. 2011) and others were situated in close proximity (43, 238, 351, 354, 420, 421, 420, 447, and 608). Of the three PSS described in TLR8 of other taxa (491 in primates (Wlasiuk and Nachman 2010), 178, 388 in suidae (Darfour-Oduro et al. 2015)) one (388) showed amino acid alterations and the other two were conserved in bats (Fig. 1, Supplementary Table S4, Fig. S3).

Evidence for episodic positive selection

We performed branch-site REL and MEME analyses to detect signatures of episodic positive selection in specific bat lineages within TLR8 evolution. We found evidence for episodic positive selection in the ancestral lineages of different higher taxonomic clades and families, including the ancestors of Yinchiroptera, Rhinolophoidae, Noctilionoidae, and Pteropodidae with strongest support for the Rhinolophoid lineage (Table 4, Fig. 2). Analyses also indicated that terminal branches of single species have been under positive selection pressure, including M. brandtii, R. hipposideros, S. bilineata, and V. murinus. The branch leading to S. bilineata might concern the whole family of Emballonuridae, because only a single species of this family was included in our analyses. The MEME analyses revealed 34 sites of episodic positive selection along branches (Fig. 2, Table 4, Supplementary Table S4). Besides lineage-specific codons, some of the identified sites coincided with codons detected to evolve as well under persistent positive selection. Position 351 involved in ligand binding and dimerization was positively selected in the ancestral branches leading to Yinchiroptera, Rhinolophoidae and within the Myotis branch. Position 539 which has been described as a ligand-binding-site mutation in bats by Escalera-Zamudio et al. (2015) was identified as a site underlying episodic positive selection in our analyses (Table 4, Supplementary Table S4, Fig. S3).

Table 4

Detection of putatively episodic positive selection within the Chiropteran TLR8 evolution

Branch-site REL analysis					MEME analysis
Branch	ω+a	Pr (ω = ω+)b	LRTc	p valued	Codons evolving under episodic positive selectione
Ancestral lineage of Rhinolophoidae	8.9	0.13	38.12	<0.0001	77-109-236-338- 351 -356- 416-418- 447-677-691- 508
Ancestral lineage of Pteropodidae	7.9	0.12	35.61	<0.0001	39-61-279-338-397-407-418-421-447-539-677-685-712-723 -730
Ancestral lineage of Noctilionoidae	90.2	0.02	25.59	<0.0001	279- 338-420-439-539
Ancestral lineage within Myotis	174.1	0.01	20.54	0.00	77 -338 -351- 356- 416-418-419-421-508 -677-685-723
Ancestral lineage of Yinpterochiroptera	14.7	0.05	14.58	0.00	32 -277-338 -351-508-677-712
S. bilineata	2.9	0.24	12.64	0.01	39-88-236-277-344-354-412-421- 685-691-723
R. hipposideros	30.2	0.02	12.54	0.01	236-386
V. murinus	7.4	0.06	12.39	0.01	88- 109- 167- 354 -397-421- 447-677-723
M. brandtii	31.8	0.02	9.67	0.03	157-419

Maximum likelihood estimation of the rate class with unconstrained ω

Maximum likelihood estimation of the proportion of sites evolving at ω+

Likelihood ratio test for ω+ = 1 (null) versus ω+ unrestricted (alternative)

p value corrected for multiple testing using the Holm–Bonferroni method

Position relative to human protein sequence; codons in bolt were detected exclusively by the MEME analysis, whereas codons in italic were identified also to evolve under persistent positive selection (see Table 3)

Fig. 2

Phylogenetic tree for Chiropteran TLR8. Neighbor-joining tree was estimated based on a distance matrix of Kimuras two-parameter evolutional model using the human TLR8 sequence (XM_005274543) as outgroup to root the tree. Bootstrap percent probabilities were based on 500 replicates. The scale bar represents a genetic distance of 0.02. Estimated lineages evolving under episodic positive selection by the branch-site REL analyses are shown in red and mean ω values are represented. Amino acid sites under episodic positive selection detected singularly by MEME analyses are mapped above branches (arrows); see also Table 4

Predicted molecule structure and location of PSS according to human TLR8

Analyses of surface accessibility confirmed that all PSS, except three (aa position 104, 420, and 447), should prepossess exposed positions in the three-dimensional molecule structure (Table 3). However, position 104 was located in a protruding loop within LRR2 and 447 was located within the Z-loop and might thus be incorrectly defined by the LRR-finder-tool.

Projection of PSS onto the 3D protein model of R. euryale revealed that the sites under selection localized to approximately four distinct regions within the LRR ectodomain (Fig. 3). The main region concerned several PSS located in a section between LRR10-18 including the Z-loop, which has been found to be involved in dimerization and in binding uridine-rich degraded ssRNA at two distinct binding sites (Tanji et al. 2013, 2015). Five PSS (312, 338, 368, 420, 421, and 481) were situated in a concavity surrounded region of LRR 10, 11, and 14 and parts of the Z-loop, which are considered to bind ssRNA oligonucleotides (Tanji et al. 2015). Here, position 338 had the strongest support by convergent findings in all our analyses and was also indicated as PSS in mammals (Areal et al. 2011, Supplementary Table S4). Out of the 21 amino acids predicted to be directly involved in RNA-binding at that site, 12 exhibited amino acid mutations and nine were conserved in bats (Supplementary Table S4). Three PSS (352, 354, and 356) were located within or next to another binding site predicted to bind the base uridine and small chemical ligands, also known to be responsible for TLR8 dimerization (Figs. 3 and 4). Here, PSS 352 has been found to interact directly with the base uridine, as does position 351, which was found to be under episodic positive selection in bats. Some of the predicted amino acids representing this uridine binding site in hTLR8 were also conserved in bats (11 out of 16 aa) and five positions showed variable amino acid substitutions (Supplementary Table S4). Interestingly, PSS 447 was located in the Z-loop within a part (position 442–457) predicted to be cleaved in hTLR8 (Tanji et al. 2013).

Fig. 3

Predicted tertiary molecule structure of a chiropteran TLR ectodomain (R. euryale) and sites of persistent positive selection. Tertiary structure was generated by Swiss-Model based on known human-crystallized TLR8 structure, Z-loop is shown in orange, cleaved part in black, ligand-binding site predicted to be involved in binding ssRNA oligonucleotids are marked in green, and the site predicted to be involved in dimerization and binding the base uridine are indicated in yellow (Tanji et al. 2015). Positions of sites under persistent positive selection in bats are shown in pink; purple indicate positions that coincide with amino acids of human binding sites. a Front, b side, and c top view

Fig. 4

Predicted tertiary molecule structure of a chiropteran TLR8-dimer (R. euryale) and sites of persistent positive selection at the dimerization interface. Tertiary structure was generated by Swiss-Model based on known human-crystallized TLR8-dimer, ligand-binding site predicted to be involved in dimerization and binding the base uridine (Tanji et al. 2015) are indicated in light (α-chain) and dark green (ß-chain). Positions of sites under persistent positive selection in bats are shown in light (α-chain) and dark (ß-chain) blue; light (α-chain) and dark (ß-chain) purple indicate positions that coincide with amino acids of the human binding sites. a Front, b side, and c top view

Besides the PSS found within the ligand-binding region, we found several PSS clustering in three additional regions of the protein (Fig. 3). One region concerned LLR1 (43) and a protruding loop within LRR2 (104, 109). A second region concerned three amino acids located next to each other within LRR7 and 8 (236, 238, 277) at the protein surface. Of those, position 236 was also described as PSS in mammals (Areal et al. 2011). A last region showed five PSS (691, 677, 685, 712, and 732) in close proximity within LRR23 and 24, also one (677) with strong support (Supplementary Table S4). Position 677, 699, and 723 were in close neighborhood to each other located at the surface towards the dimerization interface. Position 712 was also identified as PSS in mammals (Areal et al. 2011).

Discussion

Similar to TLR8 of other mammal species, the TLR8 gene of bats is highly conserved concerning gene and domain architecture, suggesting an analogous function. In all investigated bat species, the coding sequence of TLR8 is encoded within a single large exon, which is in agreement with the gene structure of other mammals. However, indels lead to structural length differences in the LRR ectodomain between different bat families. Also, the number and location of LRRs within the ectodomain varied substantially between species, suggesting alterations in molecule structures which might have an effect on ligand-binding affinities. The observed percent sequence identities at the amino acid level were highest in the TM and TIR domain when compared to that of other mammals and within bat species, indicating the importance of a conserved function in these domains for triggering the intracellular immune cascade. These contrasts to finding in the LRR ectodomain were less similarity was found at the amino acid level indicating selection mechanisms towards amino acid mutations. Altogether intraorder variability between bat taxa was high especially in the ligand-binding LRR ectodomain, and reached percent sequence identity values as different as to that to other mammal taxa, like human and horse. This finding might point to a long history of species-specific adaptation processes in the co-evolution with viral pathogens within the species rich order of Chiroptera.

As a result of the constraints imposed by the receptor protein function, purifying selection has been argued as the dominant selective force in the evolution of all TLRs, particularly for nucleic acid sensing TLRs (e.g., Georgel et al. 2009; Wlasiuk and Nachman 2010; Fornuskova et al. 2013; but see Areal et al. 2011). In bats, we found an excess of synonymous over non-synonymous substitutions when considering the whole TLR8 gene and the LRR ectodomain exclusively, which is consistent with processes of purifying selection. However, the ratio of non-synonymous over synonymous base pair substitutions reached a relatively high value of 0.517 for omega, a value that is much higher than those observed in most other proteins in mammals (range of ω = 0.1–0.25, Ellegren 2008). This fact was underlined in our analytical approach of codon evolution, where maximum likelihood models that allowed for positive selection fitted the data significantly better than models that considered only neutral sites indicating besides purifying selection also the occurrence of a positive selection pressure in the history of the TLR8 gene evolution. Overall, the large majority of codons evolved under neutrality or purifying selection, but 7 % of all sites were detected to evolve under persistent positive selection, a high value compared, e.g., to the TLR8 evolution in primates, where only 3 % evolved under positive selection (Wlasiuk and Nachman 2010) and might reflect the longer evolutionary history of bats compared to that of primates (first radiation of bats in the early Eocene compared to the Oligocene in primates).

All positively selected sites were located in the LRR ectodomain several within or next to predicted sites involved in ligand binding, which is consistent with findings in TLRs of other taxa (e.g., Areal et al. 2011; Wlasiuk and Nachman 2010; Tschirren et al. 2011, Escalera-Zamudio et al. 2015; Alcaide and Edward 2011; Grueber et al. 2012; Babik et al. 2014, Vinkler et al. 2014). As the recognition of viral RNA is essential for host defense, only such mutations that are advantageous in the recognition of co-evolving viruses should have become fixed and reflected as differences between species (Morozumi and Uenishi 2009). This was supported by the observation that the amino acid alterations within PSS led to physio-chemical changes. Physio-chemical alterations are thought to be important in defining the molecular structure and are therefore expected when natural selection promotes a functional change in the structure of proteins (Smith 2003, Vinkler et al. 2014; Darfour-Oduro et al. 2015). Our result indicate a similar or even higher level of persistent positive selection acting in TLR8 as described for non-viral TLRs in other mammalian orders (e.g., Tschirren et al. 2011; Areal et al. 2011; Smith et al. 2012, Darfour-Oduro et al. 2015) and contrast the picture of a more constrained evolution for viral TLRs as reported for other mammalian orders (humans: Barreiro et al. 2009; primates: Wlasiuk and Nachman 2010; rodents: Fornuskova et al. 2013; Suidae: Morozumi and Uenishi 2009; Darfour-Oduro et al. 2015, but see mammals: Areal et al. 2011). However, we lack information about the diversity of non-viral TLRs in bats, thus an intra-order comparison of viral and non-viral TLR evolution is still needed for the order of Chiroptea.

To get a closer insight in the strength of selection acting on TLR8 evolution within bats, we examined whether episodic positive selection has contributed to the diversity pattern. We found evidence that the branches leading to higher taxonomic clades have undergone episodic positive selection at several amino acids sites, indicating that the ancestral species within these lineages, especially those of Rhinolophoidae and Pteropodidae, had to undergo adaptive changes at these sites in response to their pathogen environment. Additionally, we detected episodic positive selection in terminal branches leading to single species (M. brandtii, R. hipposiderus, M. brandtii, and V. murinus). Our results indicate that repeated episodes of positive selection, interspersed with purifying selection shaped the observed TLR8 diversity pattern throughout the evolution of bats, and we assume that this adaptive process is currently still in operation.

The arrangement of sites under positive selection can predict locations of binding interactions between host and pathogen proteins (Werling et al. 2009; Daugherty and Malik 2012). Projection of persistent PSS onto the predicted 3D protein model revealed that the sites under selection localized to approximately four distinct regions within the LRR ectodomain (Fig. 3). Three of them concerned parts of the protein that have not been described before to be involved in ligand binding in humans (Tanji et al. 2015). However, some of them coincidence with PSS detected also in mammals (236 and 712), or were in close proximity (43 and 608) and might thus be of functional relevance. Some others were located in close neighbourhood in a protruding loop (104 and 109) or in the dimerization interface (677, 699, and 723) and might have functional relevance for structural stability, because protruding loops and the dimerization interface are thought to be involved in dimerization or stabilisation of the dimeric structure (Gautam et al. 2006; Tanji et al. 2013).

The fourth and main region with 12 PSS coincided with the predicted ligand-binding region of LLR10-18 and the Z-loop in humans (Fig. 3). Crystallographic analyses of the human TLR8 revealed that this region binds degradation products of ssRNA at two distinct sites (Tanji et al. 2015). One site within the concave surface of LRR10-14 interacts with ssRNA oligomers of different length and different sequence motifs, but seems to prefer sequences containing purine-pyrimidine dinucleotides. The other site is known to interact with small chemical ligands (e.g., imiquimod, R484) and to bind non-phosphorylated uridine mononucleoside and is furthermore important for the correct dimerization of the TLR8 dimer. These two sites have a synergistic immunostimulatory effect, because uridine-binding has to be stimulated by ligand-binding of RNA-oligonucleotides, but only uridine binding provokes reorganization of the dimer into the active conformation (Geyer et al. 2015). Thereby, it is conceivable that species-specific pyrimidine rich RNAs determine whether TLR activation is initiated or not. Some PSS (312, 338, 368, 420, 421, and 481) where located in next proximity of this binding site and three of them (236, 338, and 481), matched with PSS in mammals which did not include any chiropteran sequence (Areal et al. 2011). Here, we found strongest support for PSS at position 338. Whether these PSS might be involved in ligand binding of RNA oligomers and thus responsible for shaping a more effective immune response is highly speculative. It might also be possible that they provoke conformational differences in affecting the relation among predicted amino acids directly involved in ligand binding with the result of species specific differences in ligand-binding affinities.

Even more interesting was the location of five PSS in the site predicted to bind uridine and to be involved in dimerization and conformational activation of TLR8. Here, PSS 352 is known to directly interact with uridine. Others (354 and 356) were again in close proximity to amino acid sites predicted to interact with uridine and/or involved in dimerization. Interestingly, position 351, also predicted to interact with uridine was detected to be positive selected within the branch of Rhinolophoidae and Pteropodidae, both known as main reservoirs of ssRNA viruses. Another PSS (447) was situated within a part of the Z-loop (442–457) known to be cleaved in humans (Tanji et al. 2013), which was also defined as episodically selected in Rhinolophoidae and Pteropodidae. According to the 3D localization, a likely involvement in dimerization seems to be possible (Fig. 4). The amino acids at the cleavage site in the Z-loop (441/442) showed variable amino acid alterations in bats including a deletion of three amino acids in Vespertilioninae (corresponding to human SYA) and coincided with a deletion of five amino acids in mouse and rat (human: RQSYA). The deletion of this motif has been stated to be responsible for the dysfunction of TLR8 in rodents and has been suggested as an important determinant for distinct ligand recognition between rodent and non-rodent TLR8s (Liu et al. 2010; Govindaraj et al. 2011).

Whether and how the interspecific variation of TLR8 detected in bats play a role in TLR8 activation after binding to the RNAs from different viral pathogens is a complex question. However, it is well known that non-synonymous amino acid alterations at sites which are not directly involved in pathogen binding might indirectly influence the pathogen-binding ability of a receptor by, e.g., changing the conformation of the dimer structure in both the inactivated and activated form (Gautam et al. 2006; Tschirren et al. 2011; Geyer et al. 2015). This could also account for the three other detected PSS rich regions located outside the defined ligand-binding region of the hTLR8. There is also accumulating evidence that in distinct TLRs different species recognizes either a different ligand or that a particular TLR ligand interaction does not induce the same response within different hosts, as a result of structural differences between the proteins (e.g., TLR4: Lien et al. 2000, Lizundia et al. 2008; TLR5: Anderson-Nissen et al. 2007; TLR8: Govindaraj et al. 2011). Concerning TLR8, for instance imiquimod failed to activate hTLR8, but activated bovine and porcine TLR8 (Zhu et al. 2008, 2009); mouse TLR8 failed to respond to ligand stimulation in the absence of polyT-oligodeoxynucleotides, most probably because of the missing RQSYA motif (Gordon et al. 2006, Liu et al. 2010); and the synthetical ligand R848 induced human and sheep TLR8 activation 10 to 15-fold, whereas bovine and cat TLR8 were induced less than 5-fold (Liu et al. 2010). These variable results indicate complex mechanisms involved in TLR8 dimerization and activation which vary between different species and might concern also TLR8 function within bats. Differences in either ligand–recognition and/or receptor–ligand interactions as well as modifications in the dimer conformation are thought to reflect co-evolution of a host TLR with species-specific microbes, allowing for the detection of species-specific pathogens or discrimination between pathogenic agents and tolerated commensals (Werling et al. 2009; Govindaraj et al. 2011). Within the co-evolution of bats and their viruses, it could thus be speculated that viruses may inactivate TLR8 simply by inhibiting correct dimerization or conformation of ligand-binding sites.

Further studies are urgently needed to investigate whether TLR8 polymorphism of bats provoke any functional differences in immune response compared to humans and other mammals, which would improve our understanding of bats being potential reservoirs of certain viruses. Our findings do not necessarily elucidate how bats can remain asymptotic to viral infections, but they indicate that the TLR8 gene is subject to diversifying as well as purifying selection pressures. In addition, our results provide important insights into which codons might be candidates for influencing species-specific ligand-binding affinities of TLR8 in bats, which in turn may influence vulnerability or tolerance to specific viruses.

Below is the link to the electronic supplementary material.

Table S1 Investigated bat species, sample accession numbers due to the IZW biobank and geographic origin. Table S2 Designed primers to amplify the chiropteran TLR8 gene. Table S3 Nucleotide/amino acid sequence percent identities of the whole (a), the LRR ectodomain (b) and the TM and TIR domain (c) of the TLR8 sequences between investigated bat taxa and between the equine and human orthologous. Table S4 Positive selected sites identified by different model analyses for the Chiropteran TLR8 gene in relation to PSS of other taxa and to relevant binding an dimerization sites of the human TLR8 (hTLR8) molecule according to Tanji et al. (2013, 2015). (PDF 3309 kb)

Figure S1 Predicted protein domain architecture of bat TLR8 genes determined by SMART analyses (Schultz et al. 1998). Figure S2 Estimated transmembrane domains within the TLR8 gene. Noctilio leporinus (A) and Rhinolophus euryale (B) identified by TMPRED. Figure S3 Amino acid sequence variability of the TLR8 gene between 21 bat species and predicted positive selected sites. (PDF 111 kb)