key: cord-0961955-kon3ooo3
authors: Jiang, Haiying; Li, Juan; Li, Linmiao; Zhang, Xiujuan; Yuan, Lihong; Chen, Jinping
title: Selective evolution of Toll-like receptors 3, 7, 8, and 9 in bats
date: 2016-12-24
journal: Immunogenetics
DOI: 10.1007/s00251-016-0966-2
sha: fb0bce3654938b8c93e5b324d2a23e90f95bb15b
doc_id: 961955
cord_uid: kon3ooo3

Previous studies have shown that bats are reservoirs of a large number of viruses, many of which cause illness and mortality in humans and other animals. However, these bat-associated pathogens cause little, if any, clinicopathology in bats. This long-term adaptation should be reflected somewhat in the immune system. Toll-like receptors (TLRs) are the first line of immune defense against pathogens in vertebrates. Therefore, this study focuses on the selection of TLRs involved in virus recognition. The coding sequences of TLR3, TLR7, TLR8, and TLR9 were sequenced in ten bats. The selection pressure acting on each gene was also detected using branch- and site-specific methods. The results showed that the ancestor of bats and certain other bat sublineages evolved under positive selection for TLR7, TLR8, and TLR9. The highest proportion of positive selection occurred in TLR9, followed by TLR8 and TLR7. All of the positively selected sites were located in the leucine-rich repeat (LRR) domain, which implied their important roles in pathogen recognition. However, TLR3 evolved under negative selection. Our results are not in line with previous studies which identified more positively selected sites in TLR8 in mammalian species. In this study, the most positively selected sites were found in TLR9. This study encompassed more species that were considered natural reservoirs of viruses. The positive selection for TLR7, TLR8, and TLR9 might contribute to the adaptation of pathogen-host interaction in bats, especially in bat TLR9. ELECTRONIC SUPPLEMENTARY MATERIAL: The online version of this article (doi:10.1007/s00251-016-0966-2) contains supplementary material, which is available to authorized users.

The bats (order Chiroptera) are the second largest order of mammals, with 1150 species (www.iucnredlist.org). The order Chiroptera is divided into two suborders, Yinpterochiroptera and Yangochiroptera, according to morphological, behavioral, and molecular evidence. The Yinpterochiroptera includes rhinolophoid microbats and the megabats (pteropodids), while the other bats are seen as part of Yangochiroptera (Lei and Dong 2016; Teeling et al. 2005) . Previous studies have shown that bats are natural reservoirs of a large number of viruses (Han et al. 2015) . Approximately 150 viruses have been found in bats (Luis et al. 2013, Supplementary data) . Many of these viruses can cause illness and mortality in humans and other animals. Many infectious diseases outbroke in recent years were considered to originate from bats, such as the Severe Acute Respiratory Syndromes coronavirus (SARS-CoV; which caused a pandemic in [2002] [2003] (Xing-Yi et al. 2013) , the Middle East respiratory syndrome coronavirus (MERS-CoV; which was reported in 2013) (Wang et al. 2014) , and the Ebola virus (the culprit of the outbreak of Ebola virus disease in West Africa in 2014) (Leroy et al. 2005) . Remarkably, bats that carry viruses usually do not have biological lesions (Halpin et al. 2000; Leroy et al. 2005; Towner et al. 2007 ). Moreover, these bats did not develop overt clinical signs after experimental infection with Hendra and Nipah viruses (Middleton et al. 2007; Williamson Electronic supplementary material The online version of this article (doi:10.1007/s00251-016-0966-2) contains supplementary material, which is available to authorized users. et al. 2000) . To date, only the lyssavirus and Tacaribe virus (not zoonotic disease) were found to cause clinical signs of disease in bats (Ann et al. 2012; Johnson et al. 2008) . It seems that these bat-associated viruses cause little damage in bats themselves. We speculate that the bats have evolved certain mechanisms to prevent viral infection and control viral reproduction.

There are four types of host pattern recognition receptors (PRRs) involved in virus recognition, including Toll-like receptors (TLRs), retinoic acid-inducible gene I-like receptors (RLRs), nucleotide oligomerization domain-like receptors (NLRs), and AIM2-like receptors (ALRs) ; Thompson et al. 2011 ). Among these PRRs, the TLRs have been most extensively studied. The TLRs represent the first line of defense against pathogens and survey the extracellular and endosomal compartments. Structure of TLRs consists of a signal peptide, an N-terminal ligand recognition ectodomain (ECD), a transmembrane domain and a Cterminal intracellular Toll/interleukin-1 receptor (TIR) domain that binds signaling molecules and initiates innate cellular immune responses. The ectodomain is composed of 19-26 leucine-rich repeats (LRRs), which are ligand-binding domains, capped with the N-and C-terminal ends (LRR-NT and LRR-CT motifs) (Botos et al. 2011; Tanji et al. 2013) .

To date, 13 members of mammalian TLRs have been described. Among the TLRs, TLR2, TLR3, TLR4, TLR7, TLR8, and TLR9 are involved in virus recognition. Generally, TLR2 is activated by gram-positive bacterial lipoteichoic acid and GPI anchors of parasitic protozoans. TLR4 is well known for its critical role in the control of gram-negative bacterial infection by recognizing lipopolysaccharide (LPS). Recently, TLR2 and TLR4 are reported to mediate the cellular responses to virus infection by recognizing glycoproteins (Blanco et al. 2010; Boehme et al. 2006; Jude et al. 2003; Karen et al. 2002; Murawski et al. 2009; Rassa et al. 2002; Thompson et al. 2011) . TLR3 recognizes doublestranded RNA (dsRNA), the intermediate produced during the replication of most viruses. TLR7 and TLR8 recognize viral ssRNA. TLR9 detects unmethylated CpG motifs that are found in certain dsDNA viral, bacterial and protozoan genomes (Bowie 2007; Gupta et al. 2015; Thompson and Iwasaki 2008; Thompson et al. 2011) .

We hypothesize that host TLRs response to virus may be one of the genetic components underlying the excellent virus tolerance of bats. Signatures of an accelerated evolution of these TLRs would support this hypothesis because positive selection indicates a modification in gene function that resulted in elevated fitness. Because the predominant viral activators are nucleic acids, this study focuses on TLR3, TLR7, TLR8, and TLR9 that are responsible for detecting RNA and DNA pathogens. We sequenced the complete coding sequence (CDS) of these TLRs in ten bats and detected the selection pressure acting upon each gene.

Sample collection and RNA extraction Ten healthy bats were trapped in Guangdong or Yunnan provinces of China in 2013-2014 (Table 1) . These samples belong to three families, five Vespertilionidae species (Myotis ricketti, Myotis altarium, Ia io, Miniopterus s c h re i b e r s i i , and P i p i s t re l l u s a b r a m u s ) , t w o Rhinolophidae species (Rhinolophus affinis and Rhinolophus sinicus) two Pteropodidae species (Eonycteris spelaea and Cynopterus sphinx), and one Hipposideridae species (Hipposideros armiger). All the experimental animal programs were approved by the committee on the Ethics of Animal Experiments of the Guangdong Institute of Applied Biological Resources following basic principles. Animals were dissected after they were anesthetized with diethyl ether. The liver, kidney, lung, spleen, brain, and intestine tissues were removed quickly, preserved in RNA-Be-Lock-A (Sangon Biotech (Shanghai) Co., Ltd., China), stored at −80°C. Total RNA was extracted from frozen liver tissue using RNAiso Plus (Takara Biotechnology (Dalian) Co., Ltd., Japan).

Total RNA was reverse-transcribed into complementary DNA (cDNA) using a PrimeScript® 1st Strand cDNA Synthesis Kit (Takara Biotechnology (Dalian) Co., Ltd., Japan). The CDS of TLRs were obtained by generating overlapping amplicons. The primers were referred from previous study (Iha et al. 2010) or designed from the highly conserved regions of mammalian TLRs using Oligo 7.60 (Molecular Biology Insights, Inc., USA). The 5′ and 3′ RACE were performed using SMARTer™ RACE cDNA Amplification Kit (Clontech Laboratories, Inc., USA) and LA Taq (Takara Biotechnology (Dalian) Co., Ltd., China). The primer details are listed in the additional files (Tables S1, S2, S3, and S4). The PCR products were cloned into pMD18-T vectors and were Sanger sequenced in both directions. The sequences were assembled by SeqMan (DNASTAR, Inc., USA). All of the TLR gene sequences were submitted to GenBank (Benson et al. 2006) with the accession numbers KU302512-KU302518 and KU577291-KU577312. The complete CDS of TLRs were identified using ORF finder (http://www.ncbi.nlm.nih.gov/gorf/gorf.html) and BLAST analyses (https://blast.ncbi.nlm.nih.gov/Blast. cgi). The lengths of coding sequences for each species obtained in this study are available in additional files (Table S5 ). 

The TLR sequences of 26 other bat species, 31 representative Eutheria mammalian species, and 2 Metatheria outgroups were retrieved from GenBank (Benson et al. 2006 ). The GenBank accession numbers are listed in Table 1 . The codon sequences were aligned using MUSCLE implemented in Mega 6.0 (Tamura et al. 2013 ). All of the alignment gaps were removed. The best-fit substitution model for each dataset was selected by BIC criteria using Mega 6.0 (Tamura et al. 2013) . The maximum likelihood (ML) trees were generated based on the codon sequences by PhyML 3.0 (Guindon et al. 2010 ) with 1000 bootstraps.

The selection pressure was usually measured by the nonsynonymous/ synonymous substitution rate ratio at the NA species that lack a complete coding sequence for TLR in the NCBI database or were not identified in this study a Species in bold were sequenced in this study protein level (ω = d N /d S ). If ω = 1, the selection has no effect on fitness, the coding sequences are under neutral evolution. If ω < 1, the coding sequences are under purifying selection. If ω > 1, the coding sequences are under positive selection (Yang 2006) . To identify the selection pressure, the CODEML program was performed using PAMLX, a graphical user interface for the PAML package (Xu and Yang 2013) . The one-ratio model (model M0, assuming the same ω for all branches and sites, model = 0 and Nsites = 0) analysis was performed firstly to test the overall selection pressure in mammals. Model M0′ was performed with ω = 1 overall (model = 0, Nsites = 0, fix_omega = 1, and omega = 1). The likelihood ratio test (LRT) was performed between M0 and M0′ to test whether the TLRs of mammals are under neutral evolution.

Branch models Three two-ratio model (model = 2 and Nsites = 0) were performed. Model M2 allows the clade of Chiroptera have a uniform ω 1 , while the rest of mammals have the background ω 0 . Model M2′ allows the ω of bats to be fixed to 1 and a uniform ω 0 for the rest of mammals. In model M2″, the ω is allowed to vary among each branch of bats and set the rest of mammals as the background ω 0 . Models M2 and M0 were compared to examine whether the bat branches have different ω values compared with the rest of mammals using a LRT. Models M2′ and M2 were compared to test whether the bats are under neutral evolution. Models M2″ and M2 were compared to test whether the ω of bats branches are varied. Specific branches that evolve under positive selection were colored in red by mapping the ω value for each branch (ω > 1) onto the topologies.

Site models Models M1a, M2a, M7, and M8 (model = 0 and Nsites = 1, 2, 7, and 8) were performed to calculate the ω ratio among codons. Models M1a and M7 assume codons under neutral or purifying selection. Models M2a and M8 allow codons to evolve under positive selection. The models were compared using a LRT (M2a versus M1a and M8 versus M7) with two degrees of freedom (df = 2) (Yang 2006) . The positively selected sites were identified under Bayes empirical Bayes (BEB) with a posterior probability of P > 95%. The sites both selected by M2a and M8 were confirmed to be the positively selected sites. The datasets consisted of mammals and the dataset consisted of only bats were used to calculate the positively selected sites, respectively.

To identify the position and function of positively selected sites of TLRs, we annotated the amino acid sequences using LRRfinder (Offord et al. 2010 ) and simulated the 3-D structure of the bat TLRs using SWISS-MODEL Workspace (Arnold et al. 2006; Biasini et al. 2014 ). The PyMOL molecular graphics system (Schrodinger 2015) was used to edit the 3-D structure. The crystal structures of the human TLR3 (SWISS-MODEL Template Library (SMTL) id 2a0z.1), human TLR8 (SMTL id 3wn4.2), and horse TLR9 (SMTL id 3wpb.1) were used as templates (Bell et al. 2005; Kokatla et al. 2014; Ohto et al. 2015) . We took the horseshoe bat (R. affinis) as an example and labeled the positively selected sites on the 3-D structure models (The site numbers and the amino acids refer to complete coding sequences of R. affinis).

Seven bat TLR3 complete CDS were obtained, ranged from 2715 to 2733 bp (Table S5 ). The aligned dataset of TLR3 consisted of 50 species, including 17 bats; 2688 bp were included in the final dataset. The ML tree of TLR3 is shown in Fig. 1 . Each order of mammals was monophyletic, but the relationships among orders were weakly supported. All of bats were clustered in a clade and separated into Yangochiroptera and Yinpterochiroptera with high bootstrap values (BP ≥ 97), which was consistent with the species tree of bats (Lei and Dong 2016; Teeling et al. 2005) . The phylogenetic tree of TLR3 supported the monophyly of Laurasiatheria and Euarchontoglires. The ML tree of TLR3 placed Eulipotyphla at the root of Laurasiatheria and the Chiroptera as the sister group of Cetartiodactyla. In conclusion, the phylogenetic tree of TLR3 was consistent with species tree of mammals (Nery et al. 2012) .

The one-ratio model (M0) analysis of TLR3 gave an estimate ω = 0.26 (Table 2) . This was an average over all the mammalian branches and sites. The average ω across the mammalian tree was significantly smaller than one (model M0′), indicating strong negative selection pressure acting on mammalian TLR3 (P = 0 < 0.05, Table 2 ). To detect the selection pressure acting on bats, we performed a two-ratio model (model M2) and labeled the clade of Chiroptera as the foreground ω (ω 1 ) and the rest of the mammals as the background ω (ω 0 ). The foreground ω 1 represented the average ω over all branches and sites within bats. The results gave an estimate that both of the foreground and background ω were 0.26 (Table 2) . Model M2 did not fit the data better than model M0 (P = 1.00 > 0.05, Table 2 ). In addition, the average ω of bats was significantly smaller than one (model M2′, P = 2.81 × 10 −13 < 0.05). The small value of ω showed the domination of purifying selection in the evolution of mammalian and bat TLR3. To examine if there was any specific branch under positive selection within the bat lineage, we also performed another Btwo-ratio^model and labeled each branch of bats as the foreground ω and the rest of the mammals as the background ω (model M2″). However, none of the branches was positively selected in bats (Fig. 1) . Again, the M2″ was not a better fit than model M2 (P = 0.092 > 0.05, Table 2 ). Thus, the branch-specific models showed that the bat TLR3 evolved under strong purifying selection.

Four random-site models were also conducted to detect the positively selected sites for TLRs (Table 3) . However, the selection models (models M2a and M8) were not significantly different from the neutral models (models M1a and M7; P > 0.05). Therefore, the null hypothesis of neutral selection cannot be rejected. Again, the site-specific models provided no evidence for positive selection for bat TLR3.

For TLR7, seven bats were sequenced in this study. The complete CDS of TLR7 ranged from 3135 to 3180 bp (Table S5 ).

The alignment consisted of 3117 bp and 49 species, including 17 bats (Table 1 ). The ML tree of TLR7 is shown in Fig. 2 . Similar to the ML tree of TLR3, the monophyly of each order of mammals, Laurasiatheria and Euarchontoglires were well supported, but the relationships among orders were weakly supported. The relationships within Chiroptera were the same as those resulted from TLR3. However, the relationship among the orders was different compared with the results of TLR3. The ML tree of TLR7 showed a sister relationship of Chiroptera and Eulipotyphla with a low bootstrap value (BP = 24).

On average, the ω value of mammalian TLR7 was estimated as 0.26 which was significantly smaller than one (M0′ model, P = 0 < 0.05, Table 2 ). Therefore, the mammalian TLR7 evolved under purifying selection. In model M2, ω 1 = 0.32 and ω 0 = 0.25 (Table 2) . Model M2 fitted remarkably better than Model M0 (P = 7.74 × 10 −5 < 0.05), indicative of a divergence in the selection pressure between bats and mammals. The selection pressure for TLR7 was somewhat relaxed on bats. When the ω 1 of bats was fixed to 1, model M2′ fitted significantly less well than model M2 (P = 7.30 × 10 −10 < 0.05). This outcome suggested that the selection pressure on bat TLR7 was partially but not completely relaxed. Furthermore, we tested another Btwo-ratio^model that allowed ω vary across bat branches (model M2″) to investigate the specific positively selected branch. Model M2″ was significantly better than model M2 that ω was assumed to be uniformed among bat branches (P = 2.45 × 10 −10 < 0.05). In model M2″, the bat ancestor was detected to evolve under positive selection (Fig. 2) .

In the site-specific models, the selection models (models M2a and M8) fitted significantly better to the data than the neutral models (models M1a and M7). Model M2a provided 1.75% of sites for positive selection, while the proportion was 3.41% in model M8. Posterior probabilities of the inferred positively selected sites were estimated by the BEB procedure (Table 4 ). Model M8 identified 10-amino acid sites under positive selection with P > 95%. However, the posterior probabilities of these sites were low in model M2a except the 303N (P = 0.949). Therefore, only the 303N was considered as under positive selection for bat TLR7. But this site was not detected as positively selected site for mammalian TLR7 (Table S7) . Moreover, this positively selected site was located in LRR9 (Fig. 5) .

Briefly, the bat TLR7 evolved under negative selection with the bat ancestor and one site (303 N) were positively selected.

For TLR8, we sequenced the complete CDS of M. ricketti and R. affinis (3120 bp) and partial CDS of other seven bats (1776-3012 bp, Table S5 ). The aligned dataset consisted of 34 bats and 30 other mammals (Table 1) ; 1722 bp were included in the final dataset after removing gaps. The ML tree of TLR8 is shown in Fig. 3 . The topology of TLR8 also supported the monophyly of each o r d e r b u t d i d n o t s u p p o r t t h e m o n o p h y l y o f Euarchontoglires. The relationships within bats were the same as those of TLR3 and TLR7. However, Chiroptera was first diverged within Laurasiatheria in TLR8 gene tree.

Model M0 of TLR8 gave an estimate of the average ω over all branches and sites of mammals as 0.28. The ω of mammals M0 all branches have the same ω, M0′ all branches have the same ω = 1, M2 the Chiroptera clade has the ω 1 while other branches have ω 0 , M2′ the Chiroptera clade has the ω 1 = 1 while other branches have ω 0 , M2″ each branch of bats has ω 1 while other branches have ω 0 , lnL the natural logarithm of the likelihood value, np no. of parameters, 2Δ(lnL) twice the difference in 1nL between the two models compared was significantly smaller than one (model M0′, P = 0 < 0.05, Table 2 ), indicating strong negative selection in mammalian TLR8. In model M2, ω 1 = 0.39 and ω 0 = 0.23 (Table 2 ). In addition, model M2 was a significantly better fit than models M0 and M2′ (P < 0.05, Table 2 ). These outcomes suggested that although both of mammals and bats evolved under purifying selection on average, the selection pressure acting on bats was different from mammals. Correspondingly, eight specific branches within bats were detected as positive selection in model M2″ (Fig. 3) . As the same with the result of TLR7, the ancestor of bats was under positive selection with ω = 1.22. There were three branches under strong positive selection with ω = 999, the terminal branch of Noctilio leporinus, the ancestral lineage of Pipistrellus nathusii and Nyctalus spp., and the ancestor of Rhinolophus hipposideros, Rhinolophus ferrumequinum, and Rhinolophus euryale. In addition, the ancestor of Emballonuridae, Noctilionoidae, Pteropodidae, and Phyllostomidae, the ancestor of family Vespertilionidae, and two ancestral branches within genus of Myotis were also detected to be under positive selection (Fig. 3) .

Signatures of positive selection were also presented in the random-site models. There were 4.86% and 11.84% of sites identified as positively selected sites by models M2a and M8, respectively. Five positively selected sites were identified by the BEB procedure (P > 95%) with the models M2a and M8 (Table 4 ). Four of them were also positively selected in mammals. However, the site of 675R was specifically selected in bat TLR8. Moreover, all of the positively selected sites detected in TLR8 were located in LRRs (Fig. 5) .

In summary, bat TLR8 provided strong evidence of positive selection in both branch and site-specific models.

For TLR9, we sequenced the complete CDS of R. affinis, R. sinicus, and H. armiger (3099-3102 bp) and partial CDS of M. ricketti, I. io, and M. schreibersii (2271-2901 bp, Table S5 ). The final alignments of TLR9 consisted of 1581 bp. Fifteen bats and fourteen other mammals were included. The ML tree of TLR9 is shown in Fig. 4 . The gene trees of TLR9 showed much difference with mammalian species trees (Nery et al. 2012; Nishihara et al. 2006; Prasad et al. 2008) . The Laurasiatheria and Euarchontoglires overlapped although the support values were low. The order Lagomorpha (Euarchontoglires) was first diverged, followed by Eulipotyphla and Chiroptera (Laurasiatheria). The Primates and Rodentia were sister with Perissodactyla (Fig. 4) . However, in mammalian species trees, the Laurasiatheria and Euarchontoglires should be monophyletic (Nishihara et al. 2006; Prasad et al. 2008) .

Model M0 of TLR9 gave an estimate of the average ω across mammals as 0.47, which was significantly smaller than one (model M0′, P = 0 < 0.05, Table 2 ). Although the ω of the bat clade (ω 1 = 0.64) was estimated to be greater than that in other mammals (ω 0 = 0.42), the ω 1 was significantly smaller than one (P < 0.05, Table 2 ). This suggested that the bat TLR9 lnL the natural logarithm of the likelihood value, PSS proportion sites with ω > 1 detected by models M2a and M8, 2Δ(lnL) twice the difference in 1nL between the two models compared evolved under purifying selection on average. We then conducted model M2″ to test the selection pressure acting on specific lineages. There were seven branches evolved under positive selection, including the ancestor of Chiroptera, the ancestor of Yinpterochiroptera, the ancestral lineage of Yangochiroptera, the terminal branch of Pteropus vampyrus, the ancestral lineage of Rousettus, the ancestor of Myotis, and the ancestor of Myotis brandtii and Myotis lucifugus (Fig. 4) .

In the random-site models, a high proportion of sites were identified for positive selection, 13.06% in model M2a and 13.32% in model M8. In the BEB procedure, models M2a and M8 identified 14 and 26-amino acid sites potentially under positive selection, respectively (P > 95%). Among these selected sites, 14 sites overlapped in models M2a and M8. Therefore, these 14 sites were confirmed as positively selected sites (Table 4 ). The number of positively selected sites for TLR9 was distinctly more than that of TLR7 and TLR8. Moreover, six of these sites were specifically selected in bats while the others were common to mammals and bats. The high proportion of positively selected sites and branches inferred the strong positive selection in bat TLR9. Moreover, all of the sites under positive selection occurred in LRRs (Fig. 5 ).

Previous studies have suggested that genes involved in immune and defense systems are good candidate for positive selection because of the evolving pressure to resist pathogens (Barreiro and Quintana-Murci 2010) . However, on average, all TLRs are under negative selection with a few positively selected sites in vertebrates (Alcaide and Edwards 2011; Areal et al. 2011; Fornůsková et al. 2013; Grueber et al. 2014; Nakajima et al. 2008; Roach et al. 2005) . The selective evolution of bat TLRs have been discussed in some previous studies on single or multiple loci (Cowled et al. 2011; Escalera-Zamudio et al. 2015; Schad and Voigt 2016; Zhang et al. 2013) . In this study, we detected specific selective branches and sites in bats with half more taxa (up to 15-34 bats in different genes were included) and obtained some different results. Some of the bats belong to families where TLR sequences are not available before (e.g., Hipposideridae and Rhinolophidae). Thus, the evolutionary patterns within bats at species, family, and order levels could have been analyzed more in-depth.

According to the mammalian species tree, Laurasiatheria and Euarchontoglires should be monophyletic, but the relationships among orders within Laurasiatheria are still unresolved. In most molecular taxonomic studies, Eulipotyphla is consistently placed at the root of Laurasiatheria while the relationships of Cetartiodactyla, Perissodactyla, Carnivora, and Chiroptera are inconsistent and weakly supported (Nery et al. 2012; Nishihara et al. 2006; Prasad et al. 2008) . The gene tree of TLR3 corresponded to the mammalian species tree (Fig. 1) . As for TLR7, the phylogeny also similar to the species tree with a slight difference in the place of Eulipotyphla (Fig. 2) . The gene tree of TLR8 showed some difference compared with the species tree (Fig. 3) . However, the most different signals occurred in the phylogeny of TLR9 (Fig. 4) . This pattern was also observed in a previous study (Escalera-Zamudio et al. 2015) . The difference between the phylogeny of TLR9 and mammalian species tree might be a consequence of the accumulation of positively selected branches and sites. However, the relationships within bats were consistent among TLRs.

On average, both the mammalian and bat TLR3, TLR7, TLR8, and TLR9 were under purifying selection, consistent with previous studies (Areal et al. 2011; Cowled et al. 2011; Escalera-Zamudio et al. 2015; Schad and Voigt 2016; Zhang et al. 2013) . But there were some specific branches and sites in bat TLR7, TLR8, and TLR9 evolved under positive selection.

However, there was not any specific branch or site detected under positive selection in bat TLR3. This result is inconsistent with the description of a previous study (Escalera-Zamudio et al. 2015) . In the previous study, they gave an estimate of ω = 1 for the bat clade for the reason that the LRT between models M1 and M2 was significant (Escalera-Zamudio et al. 2015) . Actually, their model M1 is equivalent to model M2′ in the present study, while their model M2 is the same with our study. Because of the significance of LRT between models M2 and M2′ the null hypothesis (ω = 1 for the bat clade, M2′) should be rejected. Namely, our results based on branch specific models for TLR3 actually is consistent with the work of Escalera-Zamudio et al. (2015) . Moreover, plus the result of another LRT that model M2 did not significantly -Zamudio et al. (2015) showed that the LRT comparing models M8a and M8 was significant and detected six positively selected sites in bat TLR3. When we checked their result file of model M8, these sites were only selected by Naive Empirical Bayes (NEB) analysis but not selected by BEB analysis with P > 95% (Escalera-Zamudio et al. 2015) . As for mammals, five positively selected sites were detected (Table S7) . One positively selected site (87T) coincided with the dsRNA-TLR3 interaction sites (Liu et al. 2008) , which suggested that this site had an important role in recognizing viruses for mammals. The positively selected sites 4N and 507H were also detected by Areal et al. (2011) . However, the difference was that we did not detect any sites in the transmembrane and TIR domains in mammals. What is more, the five positively selected sites in mammals were detected to have a low value of ω in bat TLR3 (Table S7 ). These outcomes suggested the function conservation in bat TLR3. The crystal structure of complex between dsRNA and TLR3-ECD shows that the ectodomain of TLR3 interacts with the sugar-phosphate backbone of dsRNA, but not with individual bases, which explains the lack of specificity for nucleotide sequences in TLR3 (Liu et al. 2008) . Therefore, it is not difficult to understand the conservation of bat TLR3. In the other hand, the high conservation of bat TLR might be driven primarily by the pressure to maintain the signaling networks to binding with accessory receptor molecules and signaling intermediates (Kawasaki and Kawai 2014) rather than by the pressure of virus. With regard to TLR7, both mammals and bats were under strong purifying selection, although the ω for bats was significantly greater than the background ω for the rest of the mammals (Table 2 ). This result is consistent with previous study for the same reason analyzed for TLR3 (Escalera-Zamudio et al. 2015) . Moreover, the bat ancestor was detected to be under positive selection for TLR7. This was also proved by a study of comparative analysis of bat genomes (Zhang et al. 2013) . We identified only one positively selected site (303N) located in LRR, the ligand-binding domains for TLR7 (Botos et al. 2011 ). However, this site was not in compliance with the site descripted by Escalera-Zamudio et al. (2015) . The unique positively selected site for TLR7 detected by Escalera-Zamudio et al. (2015) (391Q) was also detected in model M8 (394N in the present study) with P (ω > 1) = 0.957 (Table 4 ). However, P (ω > 1) for 394N was low in model M2. This outcome might be a consequence of adding more bat species. However, the positively selected sites identified in mammalian and bat TLR7 were not overlapped (Table S7) . These results implied the important role of the positively selected site (303N) in virus adaptation of bats. In addition, the number of positively selected sites in mammals was similar with previous study of mammalian TLR7 with less species (Areal et al. 2011 ). However, the crystal structure of TLR7 and ligand-binding site in TLR7 has not been resolved yet.

For TLR8, the higher taxonomic branches have been positively selected, especially the two nodes with ω = 999. One of them was the ancestor of R. hipposideros, R. ferrumequinum, and R. euryale, which was also detected for episodic positive selection previously (Schad and Voigt 2016) . The other one, the ancestor of Nyctalus spp. and P. nathusii, was first detected in the present. Additionally, one single terminal species, N. leporinus was also detected for strong positive selection. The common element among TLR7, TLR8, and TLR9 was that the bat ancestor was positively selected in the long-term evolution, which was also detected by Escalera-Zamudio et al. (2015) . However, in a similar study for TLR8, the ancestor of bats was not detected (Schad and Voigt 2016) . In the site models, 675R was specifically selected in bats while the other four positively selected sites were consistent between bats and mammals (Table 4 ). This specific site for bat TLR8 was also detected by Escalera-Zamudio et al. (2015) (666C) . Moreover, all of these selected sites were located in LRRs, indicating their important role in recognizing pathogens (Botos et al. 2011) . Sironi et al. (2015) showed that some of the positively selected sites for mammalian TLR4 match to the structure differences between human and mice which is another kind of virus reservoir. The similar phenomenon was also found in bat TLR8 that 419K, 479T, and 606V were located in the LRRs involved in ligand recognition according to the crystal structure of human TLR8 (Tanji et al. 2013) .

For TLR9, seven branches within Chiroptera were under positive selection according to this study, while only Myotis showed an estimate of ω > 1 in the study by Escalera-Zamudio et al. (2015) . Furthermore, the number of positively selected sites identified in this study was significantly higher than that in previous studies. In the most related studies, four positively selected sites were identified, while fourteen sites were detected in this study (Areal et al. 2011; Escalera-Zamudio et al. 2015) . Of particular interest, six of these positively selected sites were specific in bats (Table 4) . Among them, 325T and 355S were also identified by Remarkably, all of the positively selected sites were located in LRR domains. According to the results from crystal structure of mammalian TLR9, ligands were recognized by LRRNT-LRR10 and LRR20-LRR22 (Ohto et al. 2015) . Half of the positively selected sites for TLR9 were located in these regions. The specific positively selected sites may contribute to the high tolerance of viruses in bats especially 325T, 326R, and 651K, which were located in the ligandbinding regions. The high proportion of positively selected sites and branches in bat TLR9 has never been identified. This result may be a consequence of the larger samples affected by different pathogens, for the fact that all the bats added in this study have been reported as hosts from which viruses were isolated (Luis et al. 2013) . The other reason lies on the recognition of TLR9. TLR9 not only recognizes virus DNA but also detects bacterial and protozoan genomes with unmethylated CpG motifs (Bowie 2007; Gupta et al. 2015; Thompson and Iwasaki 2008; Thompson et al. 2011) . However, the known ligand-TLR9 interaction sites are studied from a few specific ligands (Ohto et al. 2015; Tanji et al. 2013) . The identification of the exact interaction sites between the virus and TLR7, TLR8, and TLR9 is still in progress. The positively selected sites identified in this study could constitute a potential target for investigating bat viral resistance, especially the sites specifically detected in bat TLRs.

In conclusion, bat TLR3 evolved completely under purifying selection. Although the negative selection has been acting dominantly on TLRs, positively selected signals were present in bat TLR7, TLR8, and TLR9. Weak negative selection of TLRs is consistent with selection to conserve important binding specificities such as the conservation of interactive region in TLR3 (Liu et al. 2008) . The bat TLR9 showed the strongest positive selection both in branch and site models. The bat ancestor has experienced positive selection in the evolution of TLR7, TLR8, and TLR9. The positive selection for bat TLRs was a consequence of adaptation to pathogen-host interaction, especially in TLR9. Additionally, the nucleic viruses usually evolve rapidly (Domingo and Perales 2014) while the nucleic sensing TLRs show slow evolution and divergence (Escalera-Zamudio et al. 2015) . Therefore, we propose that the selection of TLRs was a combination of functional conservation and the pressure of viruses. However, we did not detect any positive selection branches in rodents, another virus reservoir (Luis et al. 2013) . Therefore, other mechanisms must exist to explain the virus tolerance of bats in addition to the TLRs. A further survey of the selective evolution of other types of pathogen recognized receptor genes in bats need to be elucidated to better understand the ability of viral resistance in bats.

Molecular evolution of the Toll-like receptor multigene family in birds

Tacaribe virus causes fatal infection of an ostensible reservoir host, the Jamaican fruit bat

Signatures of positive selection in Toll-like receptor (TLR) genes in mammals

The SWISS-MODEL workspace: a web-based environment for protein structure homology modelling

From evolutionary genetics to human immunology: how selection shapes host defence genes

The molecular structure of the Toll-like receptor 3 ligandbinding domain

SWISS-MODEL: modelling protein tertiary and quaternary structure using evolutionary information

New insights for development of a safe and protective RSV vaccine

Human cytomegalovirus envelope glycoproteins B and H are necessary for TLR2 activation in permissive cells

The structural biology of Toll-like receptors

Translational mini-review series on Toll-like receptors: recent advances in understanding the role of Toll-like receptors in anti-viral immunity

Molecular characterisation of Toll-like receptors in the black flying fox Pteropus alecto

The evolution of bat nucleic acid-sensing Toll-like receptors

Contrasted evolutionary histories of two Tolllike receptors (Tlr4 and Tlr7) in wild rodents (MURINAE)

Episodic positive selection in the evolution of avian Toll-like receptor innate immunity genes

New algorithms and methods to estimate maximumlikelihood phylogenies: assessing the performance of PhyML 3.0

Cross talk between Leishmania donovani CpG DNA and Toll-like receptor 9: an immunoinformatics approach

Isolation of Hendra virus from pteropid bats: a natural reservoir of Hendra virus

Bats as reservoirs of severe emerging infectious diseases

Molecular cloning and expression analysis of bat Tolllike receptors 3, 7 and 9

Experimental study of European bat lyssavirus type-2 infection in Daubenton's bats (Myotis daubentonii)

Subversion of the innate immune system by a retrovirus

Hemagglutinin protein of wild-type measles virus activates Toll-like receptor 2 signaling

Toll-like receptor signaling pathways

Structure-based design of novel human Toll-like receptor 8 agonists

Phylogenomic analyses of bat subordinal relationships based on transcriptome data

Fruit bats as reservoirs of Ebola virus

Structural basis of Toll-like receptor 3 signaling with doublestranded

A comparison of bats and rodents as reservoirs of zoonotic viruses: are bats special?

Experimental Nipah virus infection in pteropid bats (Pteropus poliocephalus)

Respiratory syncytial virus activates innate immunity through Toll-like receptor 2

Natural selection in the TLR-related genes in the course of primate evolution

Resolution of the Laurasiatherian phylogeny: evidence from genomic data

Pegasoferae, an unexpected mammalian clade revealed by tracking ancient retroposon insertions

LRRfinder: a web application for the identification of leucine-rich repeats and an integrative Toll-like receptor database

Structural basis of CpG and inhibitory DNA recognition by Toll-like receptor 9

Microbial sensing by Toll-like receptors and intracellular nucleic acid sensors

Confirming the phylogeny of mammals by use of large comparative sequence data sets

Murine retroviruses activate B cells via interaction with Toll-like receptor 4

The evolution of vertebrate Toll-like receptors

Adaptive evolution of virus-sensing Toll-like receptor 8 in bats

The PyMOL Molecular Graphics System, version 1

Evolutionary insights into host-pathogen interactions from mammalian sequence data

MEGA6: molecular evolutionary genetics analysis version 6.0

Structural reorganization of the Toll-like receptor 8 dimer induced by agonistic ligands

A molecular phylogeny for bats illuminates biogeography and the fossil record

Toll-like receptors regulation of viral infection and disease

Pattern recognition receptors and the innate immune response to viral infection

Marburg virus infection detected in a common African bat

Bat origins of MERS-CoV supported by bat coronavirus HKU4 usage of human receptor CD26

Experimental Hendra virus infectionin pregnant Guinea-pigs and fruit bats (Pteropus poliocephalus)

Isolation and characterization of a bat SARS-like coronavirus that uses the ACE2 receptor

PAMLX: a graphical user interface for PAML

Computational molecular evolution

Comparative analysis of bat genomes provides insight into the evolution of flight and immunity

Acknowledgements We thank our colleagues Dr. Libiao Zhang, Qi Liu, and Yi Chen for sample collection. We also thank Qianyin Mai from Guangzhou University for critical reading.

Funding This project was supported by the Planning Funds of Science and Technology of Guangdong Province (2013B031500006 and 2016B070701016) and the Funds for Environment Construction and Capacity Building of GDAS ' Research Platform (xm-2015081).

The authors have declared that no competing interests exist.