Selective evolution of Toll-like receptors 3, 7, 8, and 9 in bats.

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.

Introduction

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 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) (Pandey et al. 2015; 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 C-terminal 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 double-stranded 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.

Materials and methods

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 schreibersii, and Pipistrellus abramus), two 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).

Table 1

Sample information used in this study

Speciesa	Accession number
	TLR3	TLR7	TLR8	TLR9
Didelphimorphia
Monodelphis domestica	XM_007496035.1	XM_007500926.1	XM_007500927.1	XM_007500513.1
Dasyuromorphia
Sarcophilus harrisii	XM_012552074.1	XM_003765580.1	XM_003765579.2	XM_003762588.1
Dermoptera
Galeopterus variegatus	XM_008563573.1	XM_008571332.1	XM_008571342.1	NA
Lagomorpha
Oryctolagus cuniculus	NM_001082219.1	NA	NA	HM448400.1
Primates
Cercocebus torquatus atys	EU204935.1	EU204942.1	EU204945.1	EU204946.1
Gorilla gorilla	NM_001279752.1	KF321040.1	AB445669.1	AB445676.1
Homo sapiens	NM_003265.2	AK313858.1	BC101075.2	NM_017442.3
Macaca mulatta	NM_001036685.1	AB445665.1	AB445672.1	AB445679.1
Pan troglodytes	NM_001130470.1	NM_001130133.1	NM_001130472.1	NM_001144866.1
Pongo pygmaeus	AB445635.1	AB445663.1	AB445670.1	AB445677.1
Rodentia
Chinchilla lanigera	XM_005373438.2	XM_005411374.2	XM_005411373.2	XM_005410346.2
Cavia porcellus	NM_001173029.1	XM_003462893.3	XM_003462892.2	NA
Mus musculus	NM_126166.4	NM_001290755.1	AY035890.1	NM_031178.2
Octodon degus	XM_004639797.1	XM_004635044.1	XM_004635045.1	XM_012517625.1
Rattus norvegicus	NM_198791.1	NM_001097582.1	NM_001101009.1	NM_198131.1
Carnivora
Canis lupus familiaris	XM_540020.4	NM_001048124.2	XM_003435448.3	NM_001002998.1
Felis catus	NM_001079829.1	NM_001080133.1	EF484949.1	NM_001009285.1
Leptonychotes weddellii	XM_006741822.1	XM_006733443.1	NA	NA
Panthera tigris altaica	XM_007098990.1	XM_007090078.1	NA	XM_007081936.1
Ursus maritimus	XM_008710465.1	XM_008696372.1	XM_008696373.1	NA
Cetartiodactyla
Bos taurus	NM_001008664.1	NM_001033761.1	EF583902.1	NM_183081.1
Lipotes vexillifer	XM_007466549.1	XM_007455054.1	XM_007455053.1	NA
Ovis aries	NM_001135928.1	AM981305.1	NM_001135929.1	NM_001011555.1
Sus scrofa	NM_001097444.1	NM_001097434.1	NM_214187.1	NM_213958.1
Tursiops truncatus	XM_004321214.1	XM_004317713.1	XM_004317714.1	XM_004313247.1
Balaenoptera acutorostrata scammoni	XM_007173978.1	XM_007180334.1	XM_007180333.1	NA
Vicugna pacos	XM_006198002.1	XM_006212620.1	XM_006212621.1	XM_006220180.1
Physeter catodon	XM_007126456.1	XM_007121207.1	XM_007121208.1	NA
Eulipotyphla
Condylura cristata	XM_004682500.2	XM_004689851.2	XM_004689852.2	NA
Erinaceus europaeus	XM_007517692.1	XM_007517472.1	XM_007517473.1	XM_007519193.1
Sorex araneus	XM_004610405.1	XM_004612185.1	XM_004612186.1	XM_012934195.1
Perissodactyla
Ceratotherium simum simum	XM_004428766.1	XM_004435114.1	XM_004435115.1	NA
Equus caballus	NM_001081798.1	XM_005613953.1	NM_001111301.1	NM_001081790.2
Chiroptera
Carollia perspicillata	NA	NA	KU163598.1	NA
Desmodus rotundus	KR349157	KR349160	KR349163	KR349164
Eidolon helvum	NA	NA	KU163606.1	NA
Eptesicus fuscus	XM_008151907.1	XM_008156577.1	XM_008156576.1	NA
Miniopterus natalensis	XM_016211320.1	XM_016215635.1	XM_016202926.1	XM_016202449.1
Myotis brandtii	XM_005863096.1	XM_005880946.1	XM_005880947.1	XM_014529721.1
Myotis davidii	XM_006772707.1	XM_006763796.1	XM_006763795.1	XM_006770629.1
Myotis lucifugus	XM_006092654.1	XM_006088607.1	XM_006088606.1	XM_014448429.1
Myotis myotis	NA	NA	KU163595.1	NA
Noctilio leporinus	NA	NA	KU163601.1	NA
Noctilio albiventris	NA	NA	KU163600.1	NA
Nyctalus leisleri	NA	NA	KU163592.1	NA
Nyctalus noctula	NA	NA	KU163591.1	NA
Pteropus alecto	NM_001290169.1	NM_001290164.1	NM_001290163.1	GU045608.1
Pteropus vampyrus	XM_011363986.1	XM_011362910.1	XM_011362911.1	XM_011364642.1
Pteronotus parnellii	NA	NA	KU163599.1	NA
Pipistrellus nathusii	NA	NA	KU163593.1	NA
Rhinolophus hipposideros	NA	NA	KU163605.1	NA
Rhinolophus ferrumequinum	NA	NA	KU163604.1	NA
Rhinolophus euryale	NA	NA	KU163603.1	NA
Rousettus leschenaultii	AB472355.1	AB472356.1	NA	AB472357.1
Rousettus aegyptiacus	XM_016143919.1	XM_016119579.1	XM_016119482.1	XM_016133206.1
Saccopteryx bilineata	NA	NA	KU163602.1	NA
Vampyressa pusilla	NA	NA	KU163597.1	NA
Vampyressa bidens	NA	NA	KU163596.1	NA
Vespertilio murinus	NA	NA	KU163594.1	NA
Myotis ricketti	KU302512	KU577291	KU577298	KU577307
Rhinolophus affinis	KU302513	KU577292	KU577299	KU577308
Rhinolophus sinicus	KU302514	KU577293	KU577300	KU577309
Eonycteris spelaea	KU302515	KU577294	KU577301	NA
Ia io	KU302516	KU577295	KU577302	KU577310
Hipposideros armiger	KU302517	KU577296	KU577303	KU577311
Miniopterus schreibersii	KU302518	KU577297	NA	KU577312
Cynopterus sphinx	NA	NA	KU577304	NA
Pipistrellus abramus	NA	NA	KU577305	NA
Myotis altarium	NA	NA	KU577306	NA

NA species that lack a complete coding sequence for TLR in the NCBI database or were not identified in this study

aSpecies in bold were sequenced in this study

Complete coding sequence acquisition

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).

Phylogenetic analysis

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.

Selection detection

The selection pressure was usually measured by the nonsynonymous/ synonymous substitution rate ratio at the 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.

3-D structure of bat TLRs

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).

Results

TLR3

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).

Fig. 1

ML tree of TLR3 conducted from complete coding sequences. The bootstrap values are shown above or below the nodes. The branches with ω > 1 are labeled in red rhombus. The ω values are shown in red near the positively selected branches (conducted from model M2″)

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 “two-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.

Table 2

Likelihood ratio tests of selection pressure on TLRs under branch-specific models

Gene (no. of species, no. of sites)	Models	ω (d N/d S)	lnL	np	Models compared	2Δ(lnL)	P values
TLR3 (50, 2688)	M0	ω = 0.26	−32,563.11	99
	M0′	ω = 1	−33,756.95	98	M0′ vs. M0	2387.68	0
	M2	ω 1 = 0.26, ω 0 = 0.26	−32,563.11	100	M2 vs. M0	0	1.00
	M2′	ω 1 = 1, ω 0 = 0.26	−32,829.78	99	M2′ vs. M2	533.34	2.81 × 10−13
	M2″	Variable ω	−32,541.60	132	M2″ vs. M2	43.02	0.092
TLR7 (49, 3117)	M0	ω = 0.26	−33,696.11	97
	M0′	ω = 1	−34,886.75	96	M0′ vs. M0	2381.28	0
	M2	ω 1 = 0.32, ω 0 = 0.25	−33,688.30	98	M2 vs. M0	15.62	7.74 × 10−5
	M2′	ω 1 = 1, ω 0 = 0.25	−33,875.27	97	M2′ vs. M2	373.94	7.30 × 10−10
	M2″	Variable ω	−33,633.78	130	M2″ vs. M2	109.04	2.45 × 10−10
TLR8 (64, 1722)	M0	ω = 0.28	−27,166.03	128
	M0′	ω = 1	−28,051.22	127	M0′ vs. M0	1770.38	0
	M2	ω 1 = 0.39, ω 0 = 0.23	−27,126.50	129	M2 vs. M0	79.06	0
	M2′	ω 1 = 1, ω 0 = 0.23	−27,300.11	128	M2′ vs. M2	347.22	4.92 × 10−9
	M2″	Variable ω	−27,018.90	195	M2″ vs. M2	215.20	0
TLR9 (39, 1581)	M0	ω = 0.47	−21,463.97	77
	M0′	ω = 1	−21,697.49	76	M0′ vs. M0	467.04	0
	M2	ω 1 = 0.64, ω 0 = 0.42	−21,447.19	78	M2 vs. M0	33.56	6.91 × 10−9
	M2′	ω 1 = 1, ω 0 = 0.42	−21,469.96	77	M2′ vs. M2	45.54	1.50 × 10−11
	M2″	Variable ω	−21,407.26	106	M2″ vs. M2	79.86	7.00 × 10−7

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

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.

Table 3

Likelihood ratio tests of selection pressure on bat TLRs under site-specific models

Gene (no. of species, no. of sites)	Models	lnL	Proportion of PSS	Models compared	2Δ(lnL)	P values
TLR3 (17, 2703)	M1a	−10,187.03
	M2a	−10,187.03	0	M1a vs. M2a	0	1
	M7	−10,184.63
	M8	−10,184.35	0.79%	M7 vs. M8	0.56	0.76
TLR7 (17, 3132)	M1a	−10,730.02
	M2a	−10,723.93	1.75%	M1a vs. M2a	12.18	2.27 × 10−3
	M7	−10,732.24
	M8	−10,723.34	3.41%	M7 vs. M8	17.80	1.36 × 10−4
TLR8 (34, 3144)	M1a	−11,691.66
	M2a	−11,668.97	4.86%	M1a vs. M2a	45.38	1.40 × 10−10
	M7	−11,701.94
	M8	−11,667.04	11.84%	M7 vs. M8	69.80	6.66 × 10−16
TLR9 (15, 1602)	M1a	−7662.31
	M2a	−7625.24	13.06%	M1a vs. M2a	74.14	1.11 × 10−16
	M7	−7663.90
	M8	−7625.51	13.32%	M7 vs. M8	76.78	0

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

TLR7

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).

Fig. 2

ML tree of TLR7 conducted from complete coding sequences. The bootstrap values are shown above or below the nodes. The branches with ω > 1 are labeled in red rhombus. The ω values are shown in red near the positively selected branches (conducted from model M2″)

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 “two-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).

Table 4

Positive selection for bat TLR3, TLR7, TLR8, and TLR9 under site-specific models

Gene	PSSa	Site in alignmentb	P (ω > 1)		Domain
			M2a	M8
TLR7	238K	234	0.662	0.962*	LRR7
	282V	278	0.794	0.978*	LRR8
	303N c	299	0.949	0.994**	LRR9
	385Q	381	0.563	0.951*	LRR12
	387Y	383	0.799	0.978*	LRR12
	394N	390	0.620	0.957*	LRR12
	453A	449	0.710	0.966*	LRR14
	461Q	457	0.645	0.960*	LRR14
	493L	488	0.797	0.978*	LRR14
	737C	732	0.797	0.978*	LRR24
TLR8	419K d	15	0.989*	0.996**	LRR14
	479T d	72	0.983*	0.994**	LRR14
	545S	137	0.899	0.962*	LRR17
	606V d	198	0.996**	0.998**	LRR19
	675R c	267	0.997**	0.999**	LRR22
	697Y d	289	0.984*	0.994**	LRR23
	721R	312	0.904	0.961*	LRR24
TLR9	270R d	60	0.966**	0.992**	LRR8
	277R	67	0.871	0.961*	LRR8
	298H	88	0.894	0.971*	LRR9
	301I	91	0.862	0.958*	LRR9
	325T c	115	0.951*	0.987*	LRR10
	326R c	116	0.969*	0.992**	LRR10
	331K	121	0.915	0.977*	LRR10
	355S d	145	0.999**	1.000**	LRR11
	358P	148	0.898	0.971*	LRR11
	364H d	154	0.968*	0.992**	LRR11
	382D d	172	0.984*	0.997**	LRR12
	392E d	182	0.991**	0.998**	LRR12
	405Q c	195	0.973*	0.994**	LRR13
	407Q	197	0.913	0.977*	LRR13
	416R c	206	0.971*	0.993**	LRR13
	418C	208	0.866	0.960*	LRR14
	431P	220	0.866	0.960*	LRR14
	434R c	223	0.960*	0.990**	LRR14
	447G	236	0.838	0.950	LRR14
	451S	240	0.934	0.983*	LRR14
	461A	250	0.928	0.981*	LRR14
	547Q	333	0.921	0.978*	LRR17
	602R d	388	0.973*	0.993**	LRR20
	651K c	437	0.956*	0.989*	LRR21
	658V d	444	0.976*	0.994**	LRR22
	727K d	513	0.972*	0.993**	LRR24
	735A	521	0.936	0.983*	LRR25

*P (ω > 1) > 0.95; **P (ω > 1) > 0.99

aPositively selected sites detected by models M2a or M8 with P (ω > 1) ≥ 0.95 (scored by BEB). Sites both selected by models M2a and M8 are in bold. The amino acids refer to Rhinolophus affinis

bThe alignments of bat TLR3, TLR7, TLR8, and TLR9 can be found in the top of mlc files (result files from PAML)

cThe sites only selected in bats

dThe sites both selected in bats and mammals

Fig. 5

3-D structure and positively selected sites of bat TLR ectodomain. a TLR3, b TLR7, c TLR8, and d TLR9. The positively selected sites are shown in balls. The sites only selected in bats are shown in red. The sites both selected in bats and mammals are colored in purple. The structure and amino acids refer to Rhinolophus affinis

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

TLR8

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 order but did not support the monophyly of 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.

Fig. 3

ML tree of TLR8 conducted from partial coding sequences. The bootstrap values are shown above or below the nodes. The branches with ω > 1 are labeled in red rhombus. The ω values are shown in red near the positively selected branches (conducted from model M2″)

Model M0 of TLR8 gave an estimate of the average ω over all branches and sites of mammals as 0.28. The ω of mammals 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.

TLR9

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).

Fig. 4

ML tree of TLR9 conducted from partial coding sequences. The bootstrap values are shown above or below the nodes. The branches with ω > 1 are labeled in red rhombus. The ω values are shown in red near the positively selected branches (conducted from model M2″)

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 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).

Discussion

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 different with model M0, it could be concluded that the TLR3 evolved under purifying selection based on the results of branch specific models. In the site-specific models, Escalera-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 Escalera-Zamudio et al. (2015) (302S and 332W). 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 ligand-binding 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.

Conclusions

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.

(DOCX 17 kb)

(DOCX 17 kb)

(DOCX 14 kb)

(DOCX 15 kb)

(DOCX 14 kb)

(DOCX 15 kb)

(DOCX 32 kb)

(ZIP 102 kb)

(ZIP 1.62 MB)