Genetic drift and bottleneck do not influence diversity in Toll‐like receptor genes at a small spatial scale in a Himalayan passerine

Abstract Genetic diversity is important for long‐term viability of a population. Low genetic diversity reduces persistence and survival of populations and increases susceptibility to diseases. Comparisons of the neutral markers with functional loci such as immune genes [Toll‐like receptors; TLR] can provide useful insights into evolutionary potential of a species and how the diversity of pathogens and selection pressures on their hosts are directly linked to their environment. In this study, we compare genetic diversity in neutral (eleven microsatellite loci) and adaptive (seven TLR loci) loci to determine genetic variation in a nonmigratory western Himalayan passerine, the black‐throated tit (Aegithalos concinnus), distributed across an elevation gradient with varying degree of pathogen‐mediated selection pressure. We further compare the diversity in TLR loci with a high‐elevation sister species, the white‐throated tit (Aegithalos niveogularis). Our results indicate a lack of population genetic structure in the black‐throated tit and signatures of a past bottleneck. In contrast, we found high diversity in TLR loci and locus‐specific (TLR7) signatures of pathogen‐mediated selection, which was comparable to diversity in the white‐throated tit. Levels of diversity at TLR5 locus corresponded very closely with neutral microsatellite variation. We found evidence of positive selection in TLR1LA, TLR5, and TLR7 loci highlighting the importance in pathogen recognition. Our finding demonstrates that reduction in neutral variation does not necessarily lead to reduction in functional genetic diversity and probably helps in revival of population in a widespread species.

Mismatch in genetic diversity patterns between neutral and adaptive loci can be influenced by selective forces. For example, the selection pressure driven by pathogens is strong on functional loci and significantly influences genetic diversity at specific loci compared to effects on the neutral genetic diversity. Furthermore, both neutral and functional genetic diversity are affected by random fluctuations in the number of alleles due to genetic drift, population bottleneck, demographic history, and gene flow in a population (Takahata, 1990;Takahata & Nei, 1990;Frankham, 1996;Weber et al., 2004;Grueber, Wallis, & Jamieson, 2013;Gonzalez-Quevedo et al., 2015;Cortázar-Chinarro et al., 2017). Therefore, quantification of genetic diversity provides valuable insights on how ecological and evolutionary forces shape the genetic diversity in widespread species as well as species of conservation concern.
In vertebrates, maintenance of immunity is a complex process driven by innate and adaptive immune responses. The toll-like receptors (TLRs) are an important component of the innate immune system that plays an important role in recognizing pathogens and lead the initial fight against infections (Akira et al., 2001). The TLRs are a family of evolutionarily conserved genes, belonging to the pattern-recognition receptor (PRRs) family, that encode immune cell receptors known to recognize specific molecular ligands present in pathogens, termed pathogen-associated molecular patterns (PAMPS) (Akira et al., 2001;Keestra et al., 2013). By means of cellular cascades upon TLR activation, the adaptive immune system is subsequently triggered through antigen presentation to T cells by the Major Histocompatibilty Complex (MHC) genes (Iwasaki & Medzhitov, 2015). Species that experience a change in pathogen pressure evolve a high MHC diversity that helps to identify and eliminate a broad spectrum of pathogens (Alcaide et al., 2014;Westerdahl et al., 2004). Nonetheless, the studies involving the MHC are limited by technical challenges, especially in the case of passerines, where a high rate of gene duplication has been recorded (Minias et al., 2019). More importantly, the variations in the MHC genes alone are incapable of explaining differences in immunity, due to the collective nature of immune response, whereby different molecular players interact with the pathogen over the course of an infection (Acevedo-Whitehouse & Cunningham, 2006).
The relatively simple genetic structure, few pseudogenized and duplicated genes, and importance in being first responders to infections have led to increased emphasis on TLR genes in a wide range of taxa, including ungulates, birds, cetaceans, and wild rodents Ishengoma & Agaba, 2017;Tschirren et al., 2013). In birds, the TLR repertoire is made up of ten genes, with the gene family observed to have evolved by ancient duplication events (e.g., TLR1LA and TLR1LB, TLR2A and TLR2B;Temperley et al., 2008;Yilmaz et al., 2005). However, more recent duplication and pseudogenization events in TLR7 and TLR5, respectively, and gene conversion in two paralogues-TLR1LA and TLR1LB, TLR2A and TLR2B-have occurred independently in several passerine species (Velová et al., 2018). Distinct TLR variants have shown to be associated with diseases (Netea et al., 2012), with certain variations serving as good indicators of infection status (e.g., Borrelia infection in bank voles; Tschirren et al., 2013), disease severity (in humans, Leoratti et al., 2008) as well as survival (e.g., New Zealand robin; Grueber et al., 2012).
Parasite-mediated selection pressure can alter the diversity and allele frequencies of immune genes in many ways. For example, the TLR genes interacting directly with pathogens are often under positive selection (e.g., Alcaide et al., 2014;Darfour-Oduro et al., 2015;Tschirren et al., 2011). While purifying selection is experienced across all loci, episodic positive selection has played an important role in the evolution of most avian TLRs (Grueber, Wallis, & Jamieson, 2014;Velvová et al., 2018). Therefore, positive selection mediated by pathogens on immune genes plays an important role in shaping co-evolutionary dynamics between the host and pathogen (Antonides et al., 2019;Babik et al., 2014;Cornetti et al., 2018;Shultz & Sackton, 2019). However, in some human populations, balancing selection is the main force in shaping the evolution of immunity genes (Ferrer-Admetlla et al., 2008).
The diversity of pathogens and selection pressures on their hosts are directly linked to environment and can vary even at a small spatial scale (e.g., Cornetti et al., 2018;Gonzalez-Quevedo et al., 2016).
The western Himalayan montane system is a species-rich highly seasonal ecosystem (Srinivasan et al., 2018). The risk of Plasmodium and Haemoproteus infections in the western Himalayan birds is highest at the low elevations and decreases with increasing elevation whereas Leucocytozoon infection risk increases with elevation (Ishtiaq & Barve, 2018). Therefore, the western Himalayan birds provide a fantastic opportunity to explore the extent of functional gene diversity distributed across an elevational gradient.
In this study, we compare patterns in neutral and functional diversity in a western Himalayan passerine-the black-throated tit [Aegithalos concinnus iredalei, a 5-7 gm passerine, common and widely distributed in Pakistan, Nepal, Bhutan, India, Myanmar, Thailand, Laos, Vietnam, and mainland China, but see recent split as Aegithalos iredalei; del Hoyo et al., 2020]. The black-throated tit is a nonmigratory (sedentary) species distributed in lower to midelevation (1,800-2,600 m a.s.l.) woodland forest. The black-throated tit is a co-operative breeder favoring male-biased dispersal (Dai et al., 2013;Li et al., 2012). However, there are no studies conducted on population dispersal or gene flow in India. Using eleven neutral and seven adaptive loci in the black-throated tit, we aim to understand evolution of genes at three levels: First, we quantify and compare TLR variation between the black-throated tit and a close congeneric, ecologically similar species in the western Himalaya, the white-throated tit [Aegithalos niveogularis; a 7-8 gm passerine with restricted range in the western Himalayas-mainly Pakistan and India, Harrap, 2020]. The white-throated tit is a sedentary high-elevation species (2,800-3,200 m) but not found in sympatry with the black-throated tit. Owing to spatial variation in pathogen-mediated selection pressure, we predict that black-throated tit populations will possess a diverse set of TLR genes compared to their high-elevation counterpart, the white-throated tit. We compare signatures of selection at the haplotype and protein levels for each species. Second, we explore elevational variation in TLR haplotypes and malaria infection risk in the black-throated tit. Third, we quantify neutral variation in the black-throated tits across an elevational gradient and explore pathogen-mediated selection by comparing neutral versus TLR variation in the black-throated tit.
We explore the role of drift and signatures of genetic bottleneck in shaping the neutral genetic diversity. In addition, we compared the patterns in neutral variation in black-throated tit in India with the genotype data (using similar neutral loci) of a closely related congeneric species, Aegithalos concinnus talifuensis that is widespread across mainland China and recently split as separate species A. concinnus (Dai et al., 2017). These data allowed us to compare signatures of bottlenecks, effect of geographical distance on genetic differentiation, and population structure in two regional populations on a similar spatial scale.

| Bird blood sampling
Blood samples were collected between 2009 and 2016 during the breeding season (April-May) using mist nets in the Garhwal region, Uttarakhand, India, in the western Himalaya. The black-throated tits were sampled from five sites spanning an elevational gradient from 1,800-2,600 m a. s. l. The distance between sampled locations ranged from 5 to 150 km. The white-throated tits were sampled from one location at 3,200 m a.s.l. (see details in Table 1 and Figure 1). At each sampling site, 6-12 mist nets were set up in high bird activity locations. Each bird caught was identified and ringed with a numbered metal ring. Blood was drawn from the sub-brachial vein (blood collected never exceeded 10% of body weight) and stored in SET buffer (20-40 μL blood in 500 μL buffer 0.15 M NaCl, 0.05 M Tris, 0.001 M EDTA, pH 8.0) for subsequent analysis. Captured individuals were released at the site immediately after processing.
Therefore, we designed species-specific primers for TLR2A (ae-giTLR2; Table S1) using Primer3web v 4.1.0 (Untergasser et al., 2012) and checked for specificity using NCBI's Primer-BLAST. These primers amplified TLR2A in the white-throated tit as well. We could not amplify TLR15 and TLR21 using alternate primers as well . PCRs were set up in 25 μL reaction volumes containing a final concentration of 1× Taq Buffer, premixed with  Table S1 for details), extension at 72°C for 80 s and completed with a final extension step of 72°C for 10 min. PCR amplifications were checked on a 1.5% agarose gel and visualized with SYBR Safe DNA Gel stain (Invitrogen, Waltham, MA, USA). PCR products were cleaned by using Qiagen's PCR Purification kit (Qiagen, Hilden, Germany) and were sequenced in both directions. Sequences were analyzed using Geneious R9.

| Comparison of TLR allelic richness at nucleotide and protein levels
We compared allelic richness at both nucleotide and protein levels.
To understand selection pressure at the protein level we translated the nucleotide sequences for each TLR locus into protein sequences using the "Biostrings" package in R (Pagès et al., 2019).
Allelic diversity metrices are sensitive to sample size, so we used the Chao richness index and rarefaction curves to estimate allelic richness at nucleotide level by each locus. We calculated allelic richness after correcting for differences in sample sizes using the lowest sample size for any loci (from rarefaction curves) by species and elevation, thereby avoiding data extrapolation using the software HP-RARE v1.0 (Kalinowski, 2005). We compared estimated allelic richness in TLRs in the black-throated tit at three levels: (a) neutral loci versus TLR loci (adjusted for sample size) across elevation (low, medium or high); (b) pairwise comparisons in TLR allelic richness for each elevation category; (c) pairwise comparisons of allelic richness in the black-throated tit by elevation with sister species-the whitethroated tit. We used Welch t test to test significance of unequal variances between two means (Welch, 1951).

| Haplotype network construction
We constructed haplotype networks of protein sequences for each locus to define genetic clusters using minimum spanning trees built in Arlequin v3.5 (Excoffier & Lischer, 2010). The networks were visualized and their modularity function in Gephi (Bastian et al., 2009) was used to identify the clusters.

| Gene conversion in paralogues
Among avian TLRs, gene conversion has been observed between paralogues TLR1LA-TLR1LB and TLR2A-TLR2B (Velvová et al., 2018). We statistically tested for gene conversion events between TLR1LA and TLR1LB in both tit species using the GENECONV v1.81 (Sawyer, 1989). The program was run using the default settings of 10,000 permutations and mismatch set to zero.

| Tests for positive selection
Positive selection acting on coding sequences of TLRs was examined by estimating ω (dN/dS, the ratio of the rate of nonsynonymous (dN) to the synonymous (dS) substitution per site) using the maximum likelihood-based method in CODEML program implemented in the Phylogenetic Analysis by Maximum Likelihood (PAML v4.9; Yang, 2007). We built phylogenetic trees for each TLR loci sequenced in both Aegithalos species as input tree for using the Jukes-Cantor substitution model in MEGA X . The site models allow the ω ratio to vary among sites (Nielsen & Yang, 1998;Yang et al., 2000). CODEML implements various evolutionary models to determine selection. We used Model M0, the null model which Likelihood ratio tests (LRT) were then used to conclude any significant differences between the neutral and selection models, based on a chi-square distribution of data. We used the Bayes Empirical Bayes (BEB) (Yang et al., 2005) analysis to detect the positively se-

| Molecular sexing and parasite screening in the black-throated tit
Aegithalos species do not exhibit sexual dimorphism. Birds were sexed by amplifying the intronic region of the CHD1Z/W gene using the primers 0057F/002R following PCR conditions as described in Round et al. (2007). We screened samples for infection with three blood parasite genera (Plasmodium, Haemoproteus, and Leucocytozoon) by amplifying a partial fragment 476 bp of the cytochrome b gene as described in Hellgren et al. (2004). The infection status was confirmed by sequencing the PCR product bidirectionally.

| Spatial variation in TLR haplotypes and infection risk in the black-throated tit
We used generalized linear model (GLM) to understand the association between TLR haplotypes, TLR haplotype clusters, and TLR cluster zygosity with infection status, sex, and elevation. First, we tested if the presence of TLR haplotypes by loci varied across elevation and sex which were modeled as a fixed effect. Second, we tested if the distribution of the two codon forms in positively selected sites (as identified by CODEML and MEME selection methods) were influenced by elevation and sex, modeled as fixed effects.
Finally, we tested if infection status, as a binary response variable, was influenced by TLR haplotypes, elevation, and sex modeled as fixed effects. We assigned TLR haplotypes into clusters using haplotype network analysis ( Figures S1-S7). Thus, for each TLR locus, all haplotypes were grouped in 2-5 clusters and association of infection status was then carried out in two ways: i) with specific TLR haplotype clusters or ii) with TLR cluster zygosity, encoded as 0 (homozygous) or 1 (heterozygous) if the same or different haplotype cluster was present, respectively. We fitted ten GLMs based on the number of clusters present for the TLR locus and six models with cluster heterozygosity. We used the MuMIN package (Bartoń, 2013) to identify the best model based on quasi small-sample corrected Akaike information criterion (qAIC c ). Models with a difference (ΔAICc) of ≤2 are as parsimonious as the best-fit model (lowest AICc) and are considered strong evidence that the quality of the candidate models differs from one another (Burnham & Anderson, 2002) and the confidence interval function was used to determine the 95% confidence intervals of estimated parameters. Models were considered significant if the upper or lower confidence interval did not span zero. All GLM analyses were carried in the lme4 package (Bates et al., 2015) using a binomial error distribution and logit link function. In models that showed overdispersion of data, quasibinomial distribution was used to correct this effect. The significance of the models was tested using Wald's chi-square test (Bolker et al., 2009).

| Neutral variation in the black-throated tit
We used eleven putatively neutral microsatellites-TG01040, TG03031, TG04004, TG04041, TG01147 (Dawson et al., 2010), and Ase18, Ase37, Escµ6, Man13, Pca3, and PmaD22 (Simeoni et al., 2007), to characterize the population structure in blackthroated tit populations. These markers have been used in determining neutral variation in the black-throated tit populations in China (Dai et al., 2017). Furthermore, using similar markers allowed us to compare the genetic structure in two populations (see below in Population genetic structure). We multiplexed loci in four replicates (to account for allele dropout) with each locus tagged with fluorescently labeled dye (6-Fam, VIC or NED) at the 5' end of the forward primer (Table S3 for

| Population genetic structure in the blackthroated tit
We used five methods to assess population structure in the blackthroated tit. The first two clustering methods used a Markov chain Monte Carlo (MCMC) algorithm, built-in STRUCTURE v2.3.4 (Pritchard et al., 2000), and the discriminant analysis of principal components (DAPC) (Jombart et al., 2010). In addition, three distance-based methods were used: F ST estimations, the analysis of molecular variance (AMOVA), and isolation by distance (IBD).
Using the Bayesian clustering method STRUCTURE, we explored the possibility of each individual sample belonging to one of several clusters (K) based on allele frequencies. We ran K = 1 to K = 10, for twenty independent iterations with 500,000 MCMC steps, burn-in of 50,000, using an admixture model having correlated allele frequencies. We used Structure Harvester (Earl & vonHoldt, 2012) to determine best K following Evanno et al. (2005). CLUMPAK (Kopelman et al., 2015) was used to visualize STRUCTURE results.
For populations with high levels of gene flow or low population differentiation (F ST < 0.03), STRUCTURE often fails to accurately assign individuals to discrete populations (Latch et al., 2006). Therefore, DAPC was performed using "adgenet" package in R (Jombart, 2008).
We used PCA value of five, which successfully categorized the majority of the data with a confident classification of individuals into separate clusters by imposing a filter of membership probability q > 0.7.
We used AMOVA to check for genetic structuring within the population, by specifying individuals from genetic clusters as predicted by DAPC in Arlequin.

| Comparison of genetic variation in the blackthroated tit across two regions: India and China
To accurately compare the spatial structure and population genetic variation in the present distributional range of the species, 78 samples of Aegithalos concinnus talifuensis populations in mainland China (Dai et al., 2017) were randomly sampled from the genotype data.
We retained monomorphic loci (TG01040, TG04041, Man13) in our dataset to account for heterozygosity loss in our data. We compared these data at two levels:

| Spatial autocorrelation
We examined the effect of geographical distance on the genetic differentiation between individuals in each population by region.
Regression was performed for developing an average kinship coefficient (Fij; Loiselle et al., 1995) between pairs of individuals that were closer together (the first distance interval, ~5 km), and the logarithm of geographic distance was conducted in SPAGeDi v 1.2 (Hardy & Vekemans, 2002). Estimates close to 1 and −1 indicate high and low relatedness, respectively, and values close to zero indicate no relationship between relatedness and distance class.

| Population bottleneck
There was no evidence of population structure, so we combined all populations by regions to test for genetic bottlenecks. We used two methods to test for signatures of population bottleneck at neutral loci. First, we calculated M, a ratio based on the number of alleles to the allele size range (Garza & Williamson, 2001) in Arlequin to detect signatures of any past genetic bottleneck. The second method was used to detect recent bottlenecks using stepwise mutation model (SMM) and two-phase model (TPM), both run at variance = 12, probability = 90% and for 10,000 replicates in BOTTLENECK v1.2 (Piry et al., 1999). Significance of the models was assessed using the sign test and Wilcoxon signed-rank test for heterozygosity excess and the shape of the distribution curve indicating bottleneck was noted.

| TLR variation in the black-throated and whitethroated tits
For functional genetic diversity, we amplified 70-78 samples of the black-throated tit for TLR loci and 6-8 samples of the white-throated tit (Table 2; Tables S4 and S5). For the black-throated tit populations, an average of 13 polymorphisms resulted in a mean of 38 haplotypes across seven loci (Table S6). Haplotype and nucleotide diversity varied between 0.56-0.99 and 0.0007-0.004, respectively. However, in high-elevation sister species, the white-throated tit, an average of four polymorphisms resulted in five haplotypes averaged across seven loci (  Table 1 for population code). In the remaining sampling sites, 4-6 loci showed significant deviation from HWE after Benjamini-Hochberg FDR correction (MAN: TLR1LA, TLR2, TLR4, TLR7; ANU: TLR1LA,   TLR2, TLR3, TLR4, TLR5, TLR7; KAN: TLR1LA, TLR2, TLR3, TLR4, TLR7; p < .05 in all cases). For the black-throated tit, TLR1LB haplotypes encoded one unique protein whereas the rest of the loci coded between 7-30 unique protein sequences (Table 2 and Table S6).
Stop codons were detected in two haplotypes of TLR2 (haplotypes 2 and 3) which were prevalent in 19.2% of the total genotyped black-throated tit populations. These two haplotypes were present as homozygotes or in combination with each other in nearly half of the samples, indicating a possible overlap of TLR2 function with another locus. In the white-throated tit, haplotypes from each TLR locus translated between 2-13 unique protein sequences (Table 2).
Rarefaction analysis of TLR loci revealed that saturation of the curve was closest for TLR2 and TLR 4 in the white-throated tit. However, we had not successfully reached an asymptote for the haplotypes recorded in the remaining loci (Table S7). Abbreviations: H, number of haplotypes; H d , haplotype diversity; k, average number of nucleotide differences; NA, ω could not be determined; P, number of unique protein sequences; S, number of segregating sites; *, significance at p < .05; π, nucleotide diversity; ω, ratio of non-synonymous substitution rate (dN) to synonymous substitution rate.
for individual TLR locus in the black-throated tit, with one mutational step from each other ( Figures S1-S7). We had too few samples and haplotypes to make a reliable network for the white-throated tit.
Gene conversion analysis showed that multiple gene conversion events were detected between TLR1LA and TLR1LB in the blackthroated tit (FDR corrected p < .05), spanning between a length of 425-609 bp and one event detected in the white-throated tit, spanning 627 bp.

| Pathogen-mediated selection in TLR genes
ω values using M0 model indicated that TLRs were under purifying selection in both Aegithalos species (Table 2) Comparisons of TLR allelic richness (corrected to n = 6) across elevation in the black-throated tit and white-throated tit. Here, 1,800, 2,000 m, and 2,600 m refer to low-, medium-, and high-elevation categories in the black-throated tit and 3,200 m location refer to the whitethroated tit (p < .001 in all cases). In the black-throated tit, individual codons in TLR1LA, TLR5, and TLR7 were identified to be under positive selection. Some of the codons identified by CODEML were also validated by MEME to be under positive selection at a significance of p < .1 (Table 3). The codon identified to be under selection for TLR7 in both species was found to be at a synonymous position using protein alignment.

Analysis with GARD identified two recombination breakpoints in
all TLR genes in the black-throated tit and 2-3 breakpoints in all TLR genes of the white-throated tit (except TLR1LB, TLR2, and TLR3, which could not be tested due to less than three haplotypes).

| Spatial variation in TLR haplotypes and infection risk in the black-throated tit
Based on molecular sexing, male (n = 40) and female (n = 38) birds were nearly equally represented in the black-throated tit popula-  Figure 5) but no significant effect of sex. However, positively selected codons did not show any significant effect of elevation or sex (Table S8). For clusters and heterozygosity models, infection status was not influenced by any parameter except a negative association with elevation in TLR1LA, TLR3, TLR4 as top model (Table S9).

| Neutral variation in the black-throated tit
We successfully genotyped 78 black-throated tit samples, and allele size was determined in the range of 101-343 bp (Table 5 for   allele size and Table S10 for genotype information). Of the eleven loci used in this study, three (TG01040, TG04041, and Man13) were monomorphic and were dropped from subsequent analyses (but see Bottleneck analysis). We found no evidence of null alleles or allelic dropout in any of the loci. We found no significant difference  We retained all loci as the subsequent analysis was done across elevation. Pairwise F ST across locations was found to be very low and varied between 0 and 0.021 (Table S11). Inbreeding coefficient (F IS ) across all loci was low (mean = 0.008; Table 5

F I G U R E 3
Comparisons of allelic richness (size corrected) using eight neutral loci and seven TLR loci (*p < .05; **p < .01; and ***p < .001) in the black-throated tit

| Population genetic structure and gene flow in the black-throated tit
STRUCTURE analyses within species indicated seven clusters (K = 7) with weak genetic structure and high gene flow between populations ( Figure S8). DAPC analysis indicated six genetic clusters with considerable overlap among four clusters (Figure 6), as well as a lack of segregation by sex ( Figure S9). Nearly 88% of the individuals could be categorized with >70% membership probability to one of the six different clusters. We considered four groups for AMOVA: (a) DAPC clusters 2 and 5; (b) clusters 1 and 3; (c) cluster 4; (d) cluster 6 ( Figure 6 for DAPC cluster details). We found 13.84% of variance observed was attributable to significant differences in populations within groups (p < .001) whereas 4.25% with no significant difference at the group level (p > .05; Table S12). The Mantel test showed no evidence of any decrease in genetic differentiation with an increase in geographic distance across sampled populations (r = 0.004, p = .46; Figure S10).

| Spatial autocorrelation
We found no relationship between relatedness and spatial distance between pairs of individuals in India (b = −0.0005, SE = 0.0004;

F I G U R E 4
Relationship between pairwise population structure at Toll-like receptor (TLR) and microsatellite loci among black-throated tit populations. See Table 1 Figure S11a), and a significant negative relationship (b = −0.04. SE = 0.001; p < .001; Figure S11b) between pairs of individuals in China.

| Signatures of bottleneck
The M-ratio of allelic richness to allelic size range was low for both

| D ISCUSS I ON
We characterized genetic variation in neutral and adaptive loci in a widespread Himalayan passerine, the black-throated tits, which showed very low genetic diversity, absence of genetic structure, and evidence of a past genetic bottleneck. By contrast, diversity in TLR loci was high and exhibited locus-specific (TLR1A, TLR5, and TLR7) signatures of selection; however, there was no evidence of selection driven by malaria although this could be influenced by some other pathogens. Furthermore, we quantified variation in adaptive loci to understand selection pressure at three levels: (a) there was no change in TLR diversity within species across an elevational gradient, (b) there was no association in spatial variation in TLR diversity with infection risk, and (c) comparisons with a high-elevation sister species exhibited patterns of selection pressure shaping diversity at a specific TLR locus.

| TLR variation in two Aegithalos species
Across the seven TLR loci, we found 12-95 alleles in the black-  Xu et al., 2020) and is more comparable to those reported in widespread species (e.g., lesser kestrel, house finch, Alcaide & Edwards, 2011). We did not find contrast in TLR variation between two species. While our sample size for each elevation category was probably low to detect the selection, the levels of variation in TLR loci were relatively similar in the two sister species.
Furthermore, we found no variation between the allelic richness F I G U R E 5 Proportion of TLR7 haplotypes present in the black-throated tit showed significant variation across elevation using generalized linear model. Significant values are indicated by asterisk (*p < .05; ***p < .001) and in comparison to low-elevation category TA B L E 5 Summary of genetic diversity measures for eight polymorphic putatively neutral loci in the black-throated tit

| TLR diversity and infection status
While elevation-dependent decrease of certain pathogens has been documented (e.g., Borrelllia, Cornetti et al., 2018;Plasmodium and Haemoproteus, Ishtiaq & Barve, 2018), it is possible that high-eleva- We were not able to establish a similar relationship in our system using malaria infection data, possibly due to low numbers of infected individuals recorded and low sample size. One possible explanation for balancing selection is that of heterozygote advantage, where heterozygous states are better suited to combat a wider variety of pathogens (Doherty & Zinkernagel, 1975 but see Ejsmond & Radwan, 2015;Esjmond et al., 2014;Spurgin & Richardson, 2010).
However, we were unable to find associations between infection status and heterozygosity as well.

| Neutral versus TLR variation in the blackthroated tit
The variation in TLR loci was significantly higher than neutral loci.
We detected high diversity in TLR loci in the black-throated tits.
The significant differences in genetic diversity between neutral and adaptive markers showed the influence of pathogen-mediated selection in TLRs. Studies on bottlenecked populations demonstrated that TLR diversity was nearly the same (Gilroy et al., 2017;Gonzalez-Quevedo et al., 2015), or even less diverse than neutral loci (Hartmann et al., 2014). Both neutral and functional genetic diversity are affected by random fluctuations in the number of alleles  (Lande, 1988(Lande, , 1994. The presence of these stop codon haplotypes could have long-reaching negative biological consequences, by disruption of protein function resulting in increased susceptibility to certain infections (Hawn et al., 2003(Hawn et al., , 2005Merx et al., 2007).
Determining the role of selection on TLR diversity is however complicated by the presence of recombination sites in all TLR genes.

Recombination events can significantly influence genetic diversity
in immune genes such as MHC (Wutzler et al., 2012), but this is less common in TLRs (Alcaide & Edwards, 2011;Fisher et al., 2011;Gilroy et al., 2017;Gonzalez-Quevedo et al., 2015) and therefore warrants further investigation. Furthermore, in line with previous observa-  (Dai et al., 2017). The relatively high genetic variation in China could be attributed to hybridization between two lineages with male-biased gene flow.
While we found no relationship between inbreeding coefficient and geographic distance in India, the neighboring spatially close (5 km) individuals were more genetically related than any random individual in the 20 km distance class of China populations.
Absence of a significant pattern in India could be attributed to fragmented habitat recently colonized by a few individuals with reduced genetic diversity.
The drift-selection equilibrium can be tipped in the favor of drift influencing diversity and is prominently felt in small or reduced populations, such as those experiencing a population bottleneck (Gonzalez-Quevedo et al., 2015). Our bottleneck analyses on two regional populations showed contrasting results-the Indian population exhibited signatures of a past bottleneck whereas the China population showed a heterozygote excess as measured by two mutation models, indicating a population expansion (Cornuet & Luikart, 1996), which is in line with observations from previous studies (Dai et al., 2011). In the absence of sufficient ecological and population historical data, we cannot determine the exact timing or estimate the extent of the bottleneck event. A bottleneck was detected in the Indian population using M-ratio but the lack of heterozygosity excess shows a signature of old but severe bottleneck with a subsequent recovery. It is possible that the black-throated tit population declined in the distant past due to habitat disturbance such as forest fires or logging. Uttarakhand forests have been exposed to major forest fires every year (Chitale & Behera, 2019;Gupta et al., 2018) which could lead to fragmented habitat and influence genetic patterns and structure, local extinctions, and subsequent colonization.
In conclusion, our study shows that even though a past bottleneck event altered variations in neutral loci, diversity in most TLR genes was relatively high despite the effects of genetic drift.
Measures of functional gene diversity provide critical information in identifying vulnerable species, as reduced TLR diversity would hint at a lowered ability to recognize new or diverse pathogens, ultimately affecting individual fitness (Dalton et al., 2016). TLR diversity in the black-throated tit is high, in contrast to observations from several endangered species (Dalton et al., 2016;Gilroy et al., 2017), indicating that this species is possibly not subject to any severe threat from disease. However, the white-throated tit warrants further investigation with robust sample size. In blackthroated tits, evolution of TLRs seemed to be slightly more complex as compared to similar studies in other birds, with evidence of recombination, gene conversion, purifying, and positive selection observed (Alcaide & Edwards, 2011). Improved sampling size and comparative analysis with other adaptive loci such as MHC can provide deeper insights into the evolutionary pressure acting on immune genes. This is especially important in an environment of altering pathogenic burden caused by climate change, leading to the range expansion of parasites. Such studies can eventually prove valuable for understanding immunogenetics in migratory birds, where drastic changes in the climate, geography, and environment can affect ability to combat infections.

ACK N OWLED G M ENTS
We thank Wellcome Trust/DBT India Alliance Fellowship (IA/ I(S)/12/2/500629) to FI for financial support for this research. Our sincere thanks to Uttarakhand Forest Department for providing permission to collect avian blood samples and Dr. Y V Jhala, Wildlife Institute of India, for facilitating this research. We thank Shruti Sengupta for help with lab work and Dr. Sabuj Bhattacharyya for help during microsatellite analysis.

CO N FLI C T O F I NTE R E S T S
No competing Interests. Both authors approved the final version of the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
All data are available in the Table S2. DNA sequences with accession number MN629487-MN646077 are submitted in GenBank database.