Pathogen‐mediated selection favours the maintenance of innate immunity gene polymorphism in a widespread wild ungulate

Abstract

Toll‐like receptors (TLR) play a central role in recognition and host frontline defence against a wide range of pathogens. A number of recent studies have shown that TLR genes (Tlrs) often exhibit large polymorphism in natural populations. Yet, there is little knowledge on how this polymorphism is maintained and how it influences disease susceptibility in the wild. In previous work, we showed that some Tlrs exhibit similarly high levels of genetic diversity as genes of the Major Histocompatibility Complex (MHC), and signatures of contemporary balancing selection in roe deer (Capreolus capreolus), the most abundant cervid species in Europe. Here, we investigated the evolutionary mechanisms by which pathogen‐mediated selection could shape this innate immunity genetic diversity by examining the relationships between Tlr (Tlr2, Tlr4 and Tlr5) genotypes (heterozygosity status and presence of specific alleles) and infections with Toxoplasma and Chlamydia, two widespread intracellular pathogens known to cause reproductive failure in ungulates. We showed that Toxoplasma and Chlamydia exposures vary significantly across years and landscape features with few co‐infection events detected and that the two pathogens exert antagonistic selection on Tlr2 polymorphism. By contrast, we found limited support for Tlr heterozygote advantage. Our study confirmed the importance of looking beyond Mhc genes in wildlife immunogenetic studies. It also emphasized the necessity to consider multiple pathogen challenges and their spatiotemporal variation to improve our understanding of vertebrate defence evolution against pathogens.

Abstract

Multiple pathogens may exert antagonistic selective pressures on toll‐like receptors gene polymorphism in wild ungulates.

antagonistic effects, balancing selection, habitat heterogeneity, roe deer, Toll‐like genes

INTRODUCTION

Balancing selection is a collective term for evolutionary processes that adaptively maintain variation in populations. It is a widespread process driving the diversity of genes involved in numerous biological functions including pathogen resistance (Sommer, 2005), colour polymorphism (Wellenreuther, 2017) or mate choice/male competitiveness (Johnston et al., 2013; Moore & Moore, 1999). Disentangling the mechanisms driving balancing selection in natural systems (e.g., “heterozygote advantage” (Lewontin & Hubby, 1966), “fluctuating selection” (Hill, 1991; Levene, 1953), “negative frequency‐dependent selection” (also called rare‐allele advantage) (Phillips et al., 2018; Takahata & Nei, 1990) is particularly challenging since they can operate in synergy and are not mutually exclusive (Spurgin & Richardson, 2010).

Immunity‐related genes are ideal candidates for studying how adaptive genetic diversity is maintained in wild populations because of their direct effect on survival (Møller & Saino, 2004) or reproductive success (Pedersen & Greives, 2008). The genes of the major histocompatibility complex (MHC) are often used as a model since they show an exceptional polymorphism in vertebrates (Spurgin & Richardson, 2010). Their structure and mechanisms of interaction with antigens are well described (Klein & Figueroa, 1986). Furthermore, their ability to trigger a rapid microevolutionary response to varying pathogen‐mediated selective pressures has been demonstrated experimentally (Eizaguirre et al., 2012; Phillips et al., 2018). However, Mhc genes, although important, only correspond to a fraction of the genetic variation in pathogen resistance (Acevedo‐Whitehouse & Cunningham, 2006; Jepson et al., 1997). Before the adaptive immunity has an opportunity to intervene in the host response, invading pathogens interact with different effectors of innate immunity that recognize specific structures (Jepson et al., 1997; Lam‐Yuk‐Tseung & Gros, 2003). Being evolutionary ancient and under strong functional constraints, genes involved in innate immune responses have long been expected to evolve primarily under purifying selection and to exhibit limited polymorphism (Parham, 2003). Yet, it is now quite clear that some of these genes may have relaxed selective constraints and show large nucleotide diversity that mediates population response to emerging diseases (Lundberg et al., 2020; Seabury et al., 2010). In particular, genes encoding Toll‐like receptors (TLR) play a key role in the host frontline defence against a wide range of microparasites such as fungi, bacteria or protozoa. This family of genes is conserved in both invertebrate and vertebrate lineages, and is critical for the stimulation of many immune pathways such as inflammation, innate or acquired immunity (Akira et al., 2001). A multitude of studies have revealed associations between Tlr gene polymorphism and infectious diseases in humans and domestic species (Mun et al., 2003; Yang & Joyee, 2008). However, very few (most focussing on birds or rodents) have investigated their role in host‐pathogen interactions in natural environments and their contemporary mechanisms of evolution at a population level (but Heni et al., 2020; Kloch et al., 2018).

Moreover, infection with multiple pathogens is a real‐world rule (Maizels & Nussey, 2013) and some selective mechanisms such as heterozygous advantage or fluctuating selection may better be detected when accounting for the multi‐pathogen context (Sin et al., 2014; Wegner et al., 2003). Most immunity‐related genes have pleiotropic effects on resistance against several pathogens (Hedrick, 1999) and the maintenance of “susceptible” (deleterious) alleles in a wild population can be due to antagonist effects of a given immunity‐related gene on several pathogens (i.e., antagonistic pleiotropy). To date, the great majority of studies examining the selective effects of multiple pathogens on host immunogenetic patterns focussed on Mhc genes (Tollenaere et al., 2008 on voles; Loiseau et al., 2008 on house sparrows; Kamath et al., 2014 in plain zebras; Froeschke & Sommer, 2012 on striped mouse but see Antonides et al., 2019 in bananaquit and Heni et al., 2020 in neotropical rodents).

Here, we investigated the selective mechanisms favouring the maintenance of innate immunity genetic diversity in a free‐ranging population of European roe deer (Capreolus capreolus) inhabiting an agricultural landscape in Southwestern France. We focussed on three genes encoding pattern‐recognition receptors (PRRs) that detect specific microbes or pathogen‐associated molecular patterns: the Toll‐like receptors genes (Tlr2, Tlr4 and Tlr5). Previous work showed that roe deer populations have maintained high levels of Tlr genetic diversity, similar to those observed at Mhc genes (Quéméré et al., 2015). In particular, Tlr2 and Tlr4 exhibited several haplotypes at moderate frequency and very low levels of differentiation among roe deer populations compared to the patterns obtained with neutral markers. These findings support the hypothesis that the polymorphism of these genes is shaped by balancing selection although the underlying mechanisms remain to be investigated. Support for the “heterozygote advantage” hypothesis has been found in many Mhc gene studies of mammals (Hedrick, 2012; Oliver et al., 2009; Sin et al., 2014), including wild ungulates (Brambilla et al., 2018; Paterson et al., 1998). The importance of this hypothesis for innate immunity‐related genes has rarely been assessed (but see Antonides et al., 2019; Grueber et al., 2013) but a study in birds (Grueber et al., 2013) has found evidence for survival advantage associated with Tlr4 heterozygote genotypes. Moreover, Tlr genes are known to have pleiotropic effects on the resistance to pathogens harbouring similar pathogen‐associated molecular patterns (Hedrick, 1999; Kaisho & Akira, 2004; Loiseau et al., 2008). Therefore, immunogenetic diversity may be maintained when different pathogens exert varying directional selection on concurrent alleles of the same locus (“host immunogenetic trade‐off” hypothesis) (Kamath et al., 2014).

The study population is heavily exposed to a range of pathogens that also infect livestock, including the abortive Toxoplasma gondii or Chlamydia abortus (Candela et al., 2014). These two pathogens have negative influence on the survival and reproductive performances of ungulates (Dubey, 2009; Pioz et al., 2008). It can therefore be assumed that they exert selective pressures on immunity‐related genes and in particular on Tlr genes that play a major role in host defence against protozoans and gram‐negative bacteria (Beckett et al., 2012; Gazzinelli & Denkers, 2006; Miller et al., 2009; Netea et al., 2012; Yarovinsky, 2008). To test this hypothesis, we first established whether there are biological (age, sex) and/or environmental patterns (year, landscape structure) in Toxoplasma and Chlamydia infectious status across hosts. Then, we analysed whether pathogen prevalence was associated with specific patterns of genetic variation observed at the three immunity‐related genes (Tlrs), while accounting for biological/environmental heterogeneity. Specifically, we first tested whether heterozygous individuals exhibit a lower prevalence than homozygotes (heterozygote advantage). Then, we tested whether particular alleles are associated with lower or higher Toxoplasma and Chlamydia prevalence (directional selection). Since the two pathogens employ similar strategies to invade host cells (Romano & Coppens, 2013), we expect that they should both exert selective pressures on the same genes, but with potentially antagonistic effects if the same allele confers resistance to one pathogen and susceptibility to the other one (i.e., “antagonistic pleiotropy,” Kubinak et al., 2012).

MATERIALS AND METHODS

Study population and data collection

The study focussed on a roe deer population inhabiting a heterogeneous agricultural landscape in Southern France (43°16′N, 0°52′E). This area, called “Vallons et Côteaux de Gascogne” (VCG), is part of the « Zone Atelier Pyrénées Garonne » (https://pygar.omp.eu/). It consisted of three sectors with contrasting landscape structure with regard to the proportion of woodland: a “closed” sector with two forest blocks, an “open” landscape including mainly meadows, crops and pastures with few fragmented woodlots, and an “intermediate” sector with inter‐connected woodland fragments (Hewison et al., 2009). Neutral genetic data (both microsatellite loci and SNPs) suggested that deer from these three sectors belong to the same panmictic population (Coulon et al., 2006; Gervais et al., 2019). Individuals were caught by drive‐netting during winter (from January to March) between 2008 and 2016. Individual sex, age class (juveniles: ≤1 year of age vs. adults: >1 year of age), capture sector and body mass were recorded. Blood samples were collected for pathogen screening and ear punches for genetic analyses. In total, we gathered samples from 433 annual captures corresponding to 328 different individuals (190 females and 138 males, 164 juveniles and 164 adults, 1 to 5 captures per individual). All applicable institutional and European guidelines for the care and use of animals were followed. The study was permitted by the land manager. All the procedures involving animals were approved by the Ethical Committee 115 of Toulouse and authorized by the French Ministry in charge of ethical evaluation (n° APAFIS#7880‐2016120209523619v5).

Tlr genotyping

The Tlr (Tlr2, Tlr4 and Tlr5) genes of 157 roe deer from the area (VCG) had been genotyped in a previous immunogenetic study (Quéméré et al., 2015). We completed this dataset by genotyping 171 new individuals using exactly the same procedure. DNA was extracted from alcohol‐preserved tissues using the DNeasy Blood and Tissue kit (QIAGEN). Tlr genes were genotyped using a two‐step approach: a pilot study on 32 individuals was first performed to identify polymorphic sites (SNPs) by screening almost the entire coding region of the three Tlr genes (82% in average) including the leucine‐rich extracellular region of receptors involved in antigen‐binding (between 2081 and 2,436 bps sequenced per gene). A total of 38 SNPs (6–16 SNPs per locus) were uncovered at this step. Details about primer sequences, SNP positions and codon changes are provided in Quéméré et al. (2015). Then, we applied the exact test of linkage disequilibrium (LD) implemented in arlequin 3.5 (Excoffier & Lischer, 2010). Based on LD scores and P‐values (after Bonferroni correction), we identified linkage disequilibrium (LD) groups. For each group, we selected one SNP (primarily targeting nonsynonymous sites) that was genotyped for all individuals using the KASPar allele‐specific genotyping system provided by KBiosciences (Hoddesdon, UK). A total of 13 SNPs (out of 38) were typed including 5, 3 and 5 SNPs for Tlr2, Tlr4 and Tlr5, respectively. Details about SNP position and codon change can be found in Table S1. Lastly, haplotypes were reconstructed from the phased SNPs using the procedure implemented in dnasp v5 (Librado & Rozas, 2009). All sequences have been submitted to NCBI GenBank (Accession nos. are in Table S2).

Pathogen screening

The serological status of roe deer for Toxoplasma gondii and/or Chlamydia abortus was analysed using classical enzyme‐linked immunosorbent assays (ELISA) with specific commercial kits (Sevila et al., 2014). The specificity and sensitivity of these kits were, respectively, 97.4% and 99.4% for T. gondii and 92.2 and 89% for C. abortus. Although the kits were developed for domestic ruminants, they were shown to be reliable and efficient in wild ungulates with a high concordance (0.8) between ELISA and the Modified Agglutination Test, a reference test for Toxoplasma in roe deer (Gotteland et al., 2014). According to the manufacturer's instructions, blood samples with antibody recognition level (ARL) >30% and >40% were considered positive for T. gondii and C. abortus, respectively (see Sevila et al., 2014 for further details). In total, we obtained the Toxoplasma and Chlamydia serological status of 277 (with 74 repeated measures) and 196 (36 repeated measures) roe deer, respectively (caught between 2008 and 2016, Table 1, Table S3). The individual repeatability (R) and its standard error (SE) for Toxoplasma and Chlamydia seroprevalence (i.e., individual consistency among capture events) were calculated using the r package “rptR” (Stoffel et al., 2017). The proportion of individuals seropositive for both Toxoplasma and Chlamydia was examined in the years where both pathogens were screened (between 2008 and 2013, Table 1).

TABLE 1

Summary of the pathogen data

Pathogen	N  SAMP	N  IND	Overall seroprevalence	Seroprevalence per site			Seroprevalence per year
				Closed	Open	Mixed	2008	2009	2010	2011	2012	2013	2014	2015	2016
Toxoplasma	351	277	0.35	0.29 (85)	0.45 (175)	0.27 (91)	0.23 (40)	0.83 (42)	0.56 (43)	0.24 (50)	0.21 (42)	0.24 (38)	0.16 (31)	0.3 (49)	0.43 (16)
Chlamydia	232	196	0.16	0.15 (51)	0.11 (121)	0.2 (60)	0.00 (40)	0.14 (41)	0.37 (43)	0.28 (50)	0.05 (42)	0.00 (16)	NA	NA	NA
Co‐infection	232	196	0.08	0.08 (51)	0.08 (121)	0.07 (60)	0 (40)	0.12 (41)	0.23 (43)	0.04 (50)	0 (42)	0 (16)	NA	NA	NA

N_SAMP and N_IND are the total number of samples and individuals, respectively. Number of samples per year and sector are in brackets. NA indicates the years for which no Chlamydia data were available.

Statistical analyses

In the first step, we investigated the effects of nongenetic biological and environmental factors that possibly affect T. gondii and C. abortus exposure (binary response variable–presence/absence) using generalized linear mixed‐effect models (GLMM) with a binomial family (using a logit link function). Analyses were performed using the glmer function implemented in the “lme4” package (Bates et al., 2012) and “Mumin” v1.7.7 (Barton, 2009) in r version 3.3.3 (R Core Team, 2017). The starting model included roe deer sex (male vs. female) and age class (juvenile vs. adult) because both behaviour and host susceptibility may vary between sexes and across an individual's lifetime (Klein, 2000) and because infection persistence is expected to result in a positive age‐prevalence relationship (Gotteland et al., 2014). We also included the capture year (nine levels) and sectors (three levels) as fixed effects to account for the inter‐annual and among‐landscape variation in environmental conditions that may affect pathogen survival and thus encounter probability (e.g. climate, food resources, population density). In all models, we included individual identity as a random effect on the intercept to account for correlation among different observations of the same individual. We used a model selection procedure using Akaike information criterion with a correction for small sample sizes (AIC_C) to determine which model best explained variation in pathogen infection (models with lowest AIC). We only examined models with ΔAIC < 2 relative to the best model (Burnham & Anderson, 1998). Significant variables were retained in a reduced model for use in step 2.

In the second step, we investigated the “heterozygote advantage hypothesis” by fitting a binary fixed effect (heterozygote versus homozygote) and the “allele advantage” hypothesis by considering associations between the infection status and the number of copies (0, 1 or 2) of the most or second most frequent haplotype of Tlr2, Tlr4, Tlr5. We used variance inflation factors (VIF) to measure collinearity among predictors and retained variables with VIF values <3 (Zuur et al., 2009). Haplotypes from different Tlr genes can be included in the same model. However, to avoid redundancy, the effects of the two most frequent haplotypes of a gene and its heterozygosity were evaluated in separate models. Considering the results of step 1, the identity of year and sector (for Toxoplasma) and year (for Chlamydia) were included as random effects on the intercept in all models to control for the effects of nongenetic predictors. In total, 64 candidate models were evaluated for both Toxoplasma and Chlamydia (see Table S4 for the full list of models). The performance of the best Toxoplasma and Chlamydia models was assessed using a hold‐out validation (k‐fold random subsampling, Berrar, 2019) repeated 10 times: for each replicate, we randomly split the dataset into a training set including 90% of the data (used to fit the model) and a testing set including the 10% data remaining. We then compared observed and predicted seroprevalence classes generated based on the typical 50% cut‐off and calculated the overall accuracy (ACC) by summing the number of correctly classified values dividing by the total number of values. We calculated the mean ACC across the 10 replicates.

RESULTS

Immunogenetic variation

Among the 328 tested individuals, we isolated six different haplotypes for Tlr2, all corresponding to different amino acids sequences. Two haplotypes (Tlr2‐2: “TGCCG” and Tlr2‐1: “CATCG”) showed particularly high frequency (51% and 35%, respectively) while the four others were rare (< 5%). The two most common haplotypes belonged to two highly divergent genetic clusters and were separated by seven substitutions (including four nonsynonymous ones) (Quéméré et al., 2015). The Tlr4 gene exhibited four alleles in 328 individuals including two alleles (“CGG” and “GAG”) with moderate to high frequency (26% and 55%, respectively) and two alleles with low frequencies (10% and 7%). Lastly, we isolated six haplotypes (in 324 individuals) for the Tlr5 among which two (“GCTCG,” “ACCTG”) had high frequencies (38% and 31%, respectively), one (“ATCCG”) had intermediate frequency (18%) and the three others were rare (<5%). Observed heterozygosity was relatively high for the three Tlr genes (H_O = 0.62, 0.64 and 0.69 for Tlr2, Tlr4 and Tlr5, respectively).

Ecological patterns of parasitism

Seroprevalence of Toxoplasma reached 35% and was moderately repeatable among years at the individual level (R = 0.21, SE = 0.06 p‐value = 0.002). By comparison, seroprevalence of Chlamydia was lower (16%) but showed a higher individual consistency among capture events (R = 0.70, SE = 0.05, p‐value <0.001). Among the 232 annual captures for which both pathogens were screened (196 different deer), only 18 individuals (7%) had anti‐Toxoplasma and anti‐Chlamydia antibodies simultaneously while 93 (40%) and 36 (15%) animals were seropositive solely for Toxoplasma and Chlamydia, respectively (Table 1). The best‐fitting model for Toxoplasma included an effect of age class with seroprevalence higher in adults. It also included the capture year with an increased probability of being seropositive for Toxoplasma in 2009 and 2010 (Seroprevalence = 83% and 56%, respectively) and the habitat type with higher seroprevalence in the “open” sector (Table S5, Table S6). The best‐fitting model for Chlamydia prevalence included only a “capture year” effect with a decreased seroprevalence observed in 2009 (16%) (Table 1).

Genetic effects on pathogen infection

The inclusion of Tlr2 genotype significantly improved the Toxoplasma infection nongenetic model. The best‐fitting model (ΔAICc = −5.71 compared to nongenetic model, AICcWt = 0.416) (Table 2) included a positive effect of the number of copies of the “TGCCG” haplotype of Tlr2 (Tlr2‐2). Roe deer carrying two copies of this frequent haplotype (51%) showed a decreased Toxoplasma seroprevalence (odd ratio OR = 0.29 [0.11–0.77]) compared to individuals without this haplotype (Figure 1a, S1a, Table S7). Hold‐out cross‐validation procedure indicated that this model has a high overall accuracy in predicting Toxoplasma prevalence (mean ACC = 0.76, CI 0.72–0.81 across the 10 replicates). The model including the second most frequent haplotype “CATCG” (Tlr2‐1) was not supported (ΔAICc = 3.58 with the best model): Toxoplasma seroprevalence substantially increased for deer carrying one copy of Tlr2‐1 (OR = 2.41 [1.15–5.07]) but this effect was not additive since individuals with two copies showed similar seroprevalence (OR = 2.21 [0.75–6.48]). In a post hoc analysis, we re‐ran this model by considering the presence/absence of Tlr2‐1 (0/1) (dominant effect) rather than the number of copies (additive effect) and this improved the model performance (ΔAICc = 1.56 with the best model) (Figure 1b, S1b). Deer carrying at least on Tlr2‐1 copy had a decreased Toxoplasma seroprevalence. We did not detect any association between Toxoplasma infection and Tlr4 or Tlr5 genetic variation.

TABLE 2

Performance of the best‐fitting generalized linear mixed‐effect models of Toxoplasma and Chlamydia seroprevalence (with ΔAIC < 2 relative to the best model)

	k	Toxoplasma (n = 351)				k	Chlamydia (n = 232)
		AICc	ΔAICc	AICcWt			AICc	ΔAICc	AICcWt
nb(Tlr2−2)+age	7	382.1	−5.71	0.416	nb(Tlr2−2)+nb(Tlr5−2)	6	166.5	−7.35	0.430
nb(Tlr2−2)+h(Tlr4)	8	383.6	−4.28	0.204	nb(Tlr2−1)+nb(Tlr5−2)	6	167.8	−5.96	0.215
nb(Tlr2−2)+nb(Tlr4−4)	9	383.7	−4.15	0.192	nb(Tlr2−2)+nb(Tlr5−2)+nb(Tlr4−4)	8	168.1	−5.71	0.190
nb(Tlr2−2)+ h(Tlr5)	8	383.7	−4.11	0.188	nb(Tlr2−2)+nb(Tlr5−2)+h(Tlr4−4)	7	168.4	−5.44	0.165
Nongenetic model	5	387.65	0		Nongenetic model	2	173.7	0

“Toxoplasma” models all included age as a fixed factor and individual ID, year and capture sectors as random factors. “Chlamydia” models included animal ID and year as random factors. “n” refers to the number of seroprevalence data. The selected models occur in bold. “k” refers to the number of model parameters. “nb(haplotype)” refers to the number of copies of the haplotype. h(gene) refers to the heterozygosity status of the gene (0 for homozygotes/1 for heterozygotes).

FIGURE 1

Genetic polymorphisms at Tlr2 are associated with both Chlamydia and Toxoplasma serological status: (a) Observed seroprevalence of Toxoplasma and Chlamydia in roe deer with the number of copies (0/1/2) of Tlr2‐2, (b) Observed seroprevalence of Toxoplasma with the presence/absence of Tlr2‐1, (c) Observed seroprevalence of Chlamydia with the number of copies (0/1/2) of Tlr2‐1. Error bars indicate standard errors

Similarly to Toxoplasma, we found a strong relationship between Tlr2 genotype and Chlamydia seroprevalence (Table 2) with the best (ΔAICc = −7.35 with nongenetic model, AICcWt = 0.430) and second‐best fitted models (ΔAICc = −5.96, AICcWt = 0.215) including the number of copies of the “TGCCG” (Tlr2‐2) and “CATCG” (Tlr2‐1) haplotypes, respectively. However, we observed the opposite pattern of association than for Toxoplasma: homozygous individuals for the most common haplotype (Tlr2‐2) had an increased probability to be seropositive (OR = 3.81 [1.13–12.89]), while individuals carrying Tlr2‐1 were less likely to be seropositive for Chlamydia (OR = 0.09 [0.01, 0.82] for “CATCG” homozygous deer) (Figure 1a, c, S1a, b, Table S8). The best‐fitting model also revealed an association between Tlr5 genotype and Chlamydia seroprevalence: roe deer with the “ACCTG” (Tlr5‐2, 31%) haplotype were less likely to be Chlamydia seropositive (OR = 0.28 [0.11, 0.72]). Hold‐out cross‐validation procedure indicated that the best‐fitting model has a high overall accuracy in predicting Chlamydia prevalence (mean ACC = 0.85, CI 0.80–0.89 across the 10 replicates). No association was observed when considering Tlr4 genotypes (neither the number of alleles nor the heterozygosity status).

DISCUSSION

Despite strong demographic fluctuations throughout the quaternary related to both natural and anthropogenic founder events (de Jong et al., 2020), roe deer has maintained similarly high levels of Tlr genetic diversity as Major Histocompatibility Complex (MHC). Previous work showed that this high sequence and allelic diversity could have been preserved through balancing selection (Quéméré et al., 2015). In this study, we developed a multi‐gene/multi‐pathogen association approach to bring new insights into the mechanisms underlying such balancing selection. Our results suggested a key role of pathogen‐mediated directional selection on specific Tlr haplotypes rather than heterozygote advantage. We showed that different Tlr genes may be associated with resistance to the same microparasite (here Tlr2 and Tlr5 with Chlamydia) in agreement with the general idea that pathogens interact with many immune recognition receptors. We also revealed that different pathogens may exert antagonistic selective pressures on the polymorphism of a given gene (here Toxoplasma and Chlamydia on Tlr2‐1 and Tlr2‐2 alleles, respectively), that may have favoured the maintenance of Tlr2 immunogenetic diversity at the population level through balancing selection.

High temporal and spatial variation in pathogen exposure

Our results revealed a high temporal heterogeneity in both pathogen seroprevalence at individual and population levels. In contrast with the traditional view of lifelong persistence of Toxoplasma infection (Tenter et al., 2000), we observed low individual repeatability of Toxoplasma seroprevalence at the within‐individual level (R = 0.21). This is in line with a previous study on the same population that reported frequent seroconversion with initially seropositive individuals becoming seronegative within 1 to 3 years (Sevila et al., 2014). Toxoplasma seroprevalence increased with roe deer age as frequently observed in wild ungulate populations, but it is most likely due to repeated exposure from the same environment rather than antibody persistence (Gotteland et al., 2014; Opsteegh et al., 2011). Overall, this suggests that the presence of antibodies in the study population would reflect relatively recent rather than long‐term infection. The observed annual variation in pathogen prevalence at the population level may be related to both climate and host factors. Meteorological conditions such as temperature or humidity are known to affect pathogen survival in the environment and may influence host transmission via the quantity of infected aerosols for Chlamydia (Tang, 2009) or oocysts for Toxoplasma (Gotteland et al., 2014). For example, Sevila et al. (2014) showed that Toxoplasma prevalence in roe deer was higher in mild and wet years and that Chlamydia prevalence increased in cold years. The high turnover of individuals in the population in relation to a strong harvest pressure may also partly explain why we observed such high inter‐annual variability in multi‐cohort samples (Candela et al., 2014). Moreover, we also observed strong variation in Toxoplasma exposure across sampling sectors that most likely results from landscape heterogeneity and in particular to varying proportions of human settlements. Indeed, exposure increases in the vicinity of human dwellings as a proxy of domestic cat presence, cats being the definitive host for T. gondii (Sevila et al., 2014).

Predominant role of directional selection on Tlr gene variation

We revealed that Tlr2 gene polymorphism is partly shaped by directional selection exerted by both Toxoplasma and Chlamydia. These two abortive pathogens may have significant impact on the annual reproductive success of ungulates (Pioz et al., 2008). TLRs have been shown to activate the complement system, a component of innate immunity known to be important for resistance to microparasites during early infection (Raby et al., 2011). A higher affinity of a Tlr haplotype to Toxoplasma or Chlamydia ligands may result in a stronger complement response and thus in an increased resistance to these pathogens (Tschirren et al., 2013).

Spatial and temporal variations in pathogen exposure suggest that the role of Tlr2 haplotypes in roe deer susceptibility to Toxoplasma and Chlamydia may vary among years and among landscape structures where females established their home range. This supports the general idea that different habitats in terms of biotic and abiotic factors are likely to support distinct pathogen communities and so distinct pathogen‐mediated‐selection regime as observed for Mhc (Alcaide et al., 2008; Eizaguirre & Lenz, 2010; Hedrick, 2002). Further data using longer time series are required to explicitly test whether Toxoplasma and Chlamydia prevalence are negatively correlated in space and/or time, which is a prerequisite for investigating whether selection pressures are fluctuating (Spurgin & Richardson, 2010; Heni et al., 2020).

Antagonistic pleiotropy and potential immunogenetic trade‐off

Interestingly, we observed that different haplotypes of the same Tlr gene (here Tlr2) may be associated with the infection status of the two pathogens. Pleiotropic effects in Tlr2 were expected since this receptor is known to be implicated in the recognition of many bacteria, fungus and protozoa (de Oliviera Nascimento et al., 2012) including Toxoplasma (Yarovinsky, 2008) and Chlamydia (Beckett et al., 2012; Darville et al., 2003). Here, the most common Tlr2 haplotype was associated with decreased Toxoplasma infection and increased Chlamydia infection while the concurrent haplotype showed the opposite trend. This result provides an example of antagonistic pleiotropy, which has been proposed as a widespread mechanism favouring the maintenance of genetic variation within populations (Roff & Fairbairn, 2007; Rose, 1982), in particular at immunity‐related genes (Aderem & Ulevitch, 2000; Carter & Nguyen, 2011). Antagonistic pleiotropy has rarely been evidenced in natural populations (Turner et al., 2011), and, to our knowledge, this study is the first demonstration with regard to a Tlr gene. Two nonexclusive mechanisms may be invoked to explain this antagonistic effect: decrease of intracellular pathogen competition (Loiseau et al., 2008) and host immunogenetic trade‐off (Kamath et al., 2014). Toxoplasma and Chlamydia use similar strategies to invade cells route and compete for the same nutrient pools in co‐infected cells (Romano & Coppens, 2013). Therefore, antagonistic effects might arise because the Tlr2 gene alters the competitive interactions between the two pathogens. In other words, because individuals carrying the most frequent Tlr2 haplotype are less infected by Toxoplasma, this may indirectly promote infection by Chlamydia without a direct role of the second most frequent haplotype (and vice versa). In a study of house sparrows, Loiseau et al. (2008) found antagonistic effects of a Mhc class I gene on multiple malarial parasite strains. In this particular case, deleterious “susceptibility alleles” could be maintained in the population due to within‐host competition between malaria parasite strains. Here, because we observed relatively few co‐infections, it is unlikely that antagonistic effects may result from the competitive interactions between Toxoplasma and Chlamydia.

The most likely hypothesis is the presence of a host immunogenetic trade‐off (Kamath et al., 2014). In this case, pathogens do not directly compete but Tlr2 genotype has a direct effect on the resistance/susceptibility to both pathogens. The allele inferred to be beneficial for decreasing Toxoplasma infection is also associated with increased susceptibility to Chlamydia (and vice versa). This scenario is reinforced by the fact that the two concurrent Tlr2 haplotypes exhibit very distinct DNA and amino‐acid sequences, which suggests marked functional differences. Similar patterns and processes have previously been described for tick and nematode infections (Kamath et al., 2014; Turner et al., 2011). To our knowledge, our study is the second one showing evidence of antagonistic selection at a non‐Mhc gene in a wild population (see for the other one Turner et al., 2011).

No evidence of “heterozygous advantage”

Our data provided weak support for the “heterozygous advantage hypothesis.” We noted that most authors revealing heterozygous advantage in immunity‐related genes (most often on Mhc genes) studied animals co‐infected with multiple pathogens (Hughes & Nei, 1992; Froeschke & Sommer, 2012; McClelland et al., 2003 but see Oliver et al., 2009). Roe deer classically inhabit a single home range of relatively small size (mean 0.76 km²) across their entire adult life (Morellet et al., 2013). We observed few co‐infection events that suggest that individuals are rarely exposed to both pathogens Toxoplasma and Chlamydia across their lifetime, explaining the low fitness benefit of heterozygous genotypes. A limitation of this work is that we could not investigate the importance of other mechanisms such as negative frequency‐dependent selection, because of insufficient statistical power to test the role of rare variants and the lack of long time series (see Charbonnel & Pemberton, 2005 for an example).

CONCLUSION

In conclusion, our study illustrates the importance of looking beyond Mhc genes (Acevedo‐Whitehouse & Cunningham, 2006) and single snapshot gene/pathogen association (Maizels & Nussey, 2013) in wildlife immunogenetic studies. In the present case, we highlighted the importance of antagonistic effects on the maintenance of innate immunogenetic diversity, which could not have been possible using a classical single gene‐single pathogen approach. Decreasing sequencing costs now offer new opportunities to examine multiple genes and diseases in a single study, leading to a more realistic view of the functioning of the immune system and of the ecological context in which it evolves (Becker et al., 2020; DeCandia et al., 2020).

ACKNOWLEDGEMENTS

This work was supported by the French National Institute for Agricultural and Environmental Research (INRAE). We thank the local hunting associations, the Fédération Départementale des Chasseurs de Haute‐Garonne. We thank numerous co‐workers and volunteers for their assistance in field data collection, in particular J‐M. Angibault, B. Cargnelutti, J‐L. Rames, J. Joachim, B. Lourtet, D. Picot and N. Cebe. We also thank Karine Laroucau and Rachid Aaziz from ANSES the laboratory for animal health in Maison‐Alfort (France) for their support and expertise on C. abortus results. We also thank Monica G Candela, Julie Sevila and co‐workers for their support on laboratory work (serology). Genetic data used in this work were produced through molecular genetic analysis technical facilities of the labex “Centre Méditerranéen de l'Environnement et de la Biodiversité.”

CONFLICT OF INTEREST

The authors have no conflict of interest to declare.

AUTHORS’ CONTRIBUTIONS

EQ, EGF and NC designed the study and interpreted the data; HV, NM and YC led the data collection process on the field; MG, JM and EQ carried out the Tlr and Mhc‐Drb genotyping; EQ, PH, NM and MF performed the statistical analyses; EQ led the writing, all authors contributed to the manuscript and gave final approval for publication.

PEER REVIEW

The peer review history for this article is available at https://publons.com/publon/10.1111/jeb.13876.

DATA AVAILABILITY STATEMENT

Haplotype DNA sequences: Primers and GenBank accessions are in Table S1. SNP positions and characteristics are in Table S2. Genotypes of all individuals at each Tlr gene, and environmental/parasitological data and R scripts can be found in Dryad at https://doi.org/10.5061/dryad.7m0cfxpth.