Genetic epidemiology of the Alpine ibex reservoir of persistent and virulent brucellosis outbreak

While it is now broadly accepted that inter-individual variation in the outcomes of host–pathogen interactions is at least partially genetically controlled, host immunogenetic characteristics are rarely investigated in wildlife epidemiological studies. Furthermore, most immunogenetic studies in the wild focused solely on the major histocompatibility complex (MHC) diversity despite it accounts for only a fraction of the genetic variation in pathogen resistance. Here, we investigated immunogenetic diversity of the Alpine ibex (Capra ibex) population of the Bargy massif, reservoir of a virulent outbreak of brucellosis. We analysed the polymorphism and associations with disease resistance of the MHC Class II Drb gene and several non-MHC genes (Toll-like receptor genes, Slc11A1) involved in the innate immune response to Brucella in domestic ungulates. We found a very low neutral genetic diversity and a unique MHC Drb haplotype in this population founded few decades ago from a small number of individuals. By contrast, other immunity-related genes have maintained polymorphism and some showed significant associations with the brucellosis infection status hence suggesting a predominant role of pathogen-mediated selection in their recent evolutionary trajectory. Our results highlight the need to monitor immunogenetic variation in wildlife epidemiological studies and to look beyond the MHC.

The epidemiological dynamics of infectious diseases results from a complex interplay between the pathogen, the host and the environment 1 . Often neglected, inter-individual variation in host susceptibility plays a key role in this dynamics since, depending on their biology (e.g. age, sex), their spatial or social behaviours and/or their genetic characteristics, some individuals contribute far more than others to disease spread 2 . Therefore, identifying individuals and their characteristics that are most responsible for disease transmission is an important step for improving epidemiological surveillance and management of wildlife diseases 3 .
Host genetic and immunogenetic characteristics are rarely investigated in wildlife epidemiological studies despite evidence that inter-individual variation in the outcomes of host-pathogen interactions is at least partially genetically controlled 4 . Pathogen-mediated selective pressures shape the genetic components of host immunity and give rise to inter-individual variation in resistance to infectious diseases 5 . This immunogenetic variation is particularly well documented in humans and domestic animals but much less in wild animals. Furthermore, most studies focused solely on the adaptive branch of the immune system, particularly on the major histocompatibility complex (MHC) 6,7 . Although there is no doubt that MHC plays a key role in individual susceptibility to many diseases 8 , it accounts for only a fraction of the genetic variation in pathogen resistance 5,9 . Before the adaptive immunity intervenes in the host response, other genes involved in the innate immunity such as the Toll-like receptors (TLRs) genes are major players in the host frontline defense against a wide range of microparasites including vectors of zoonotic diseases such as Coxiella brunetii (Q fever) or Borrelia sp. (Lyme disease) [10][11][12] .
Genetic variation at the MHC and other immunity-related genes encoding pathogen recognition receptors is shaped by both neutral (i.e. genetic drift, mutation, migration) and selective processes. The relative importance of these evolutionary forces depends on both the population demographic history that influences the strength of the genetic drift, and the (multi)pathogen challenge that determines the type and intensity of selection Results neutral genetic diversity. The tests for Hardy-Weinberg equilibrium, after correction for multiple testing, revealed that SR-CRSP07 STR locus showed a significant excess of homozygotes. This marker was thus removed from further analysis. H E was 0.43 ± 0.17 (mean ± standard error) and mean N A was 3.12 over all microsatellite loci. Neutral variation did not vary significantly before (H E = 0.42 ± 0. 16 immunity-related gene polymorphism. The Illumina sequencing of the exon 2 of the Mhc-Drb gene on 144 individuals revealed a unique haplotype in Bargy corresponding to CaIb-DRB*1 haplotype, previously isolated by Schalschl et al. 39 and Grossen et al. 40 . The Sanger sequencing of the three Tlr genes generated an average of 2,169 bp of sequence covering almost the entire coding region of the genes (92% on average). Tlr genes exhibited seven SNPS among which five were non-synonymous. We isolated 1 SNPs for Tlr1 and revealed two functional haplotypes (i.e. encoding different amino acid sequences) including a rare haplotype (5%) observed in only 14 individuals. The Tlr2 gene exhibited three SNPs (two non-synonymous) and showed three functional haplotypes with a common haplotype (Tlr2a = 72%) and two haplotypes with low frequencies (Tlr2b = 15% and Tlr2c = 13%). Lastly, we found three SNPs (two non-synonymous) in Tlr4 gene and identified four functional haplotypes among which three haplotypes showed similar moderate frequencies ranging between 25% and 35% while the fourth one was rare (9%). Expected heterozygosity (H E ) (given differences in haplotype frequencies) varied greatly among Tlrs ranging from 0.09 for Tlr1 to 0.71 for Tlr4 (Table 1). The SLC11A1 gene showed two alleles ("A324" and "A330" hereafter) with frequencies of 75% and 25% respectively (H E = 0.37).
Genetic effects on brucellosis infection. Out of the 237 ibex from the Bargy massif that were first tested, 89 were found seropositive to brucellosis (37%). The probability to be seropositive was lower in males than in females and reached a maximum in the core area of the massif (SSU3, SSU4) (see Tables S4,S5, Figs. 1,2). We observed a slight temporal decrease of the seroprevalence over the monitoring period (only significant for the dataset covering the whole 2012-2017 study period). Brucellosis seropositivity also tended to increase with age but 95% confidence intervals of the effect size did not exclude zero (Average Effect Size = 0.41 [95% CI: −0.13,0.95]) ( Fig. 1, Table S4). We did not find any significant effect of the multi-locus heterozygosity (MLH) on the individual serological status (general effect) (Fig. 2, Table S5) but we observed associations with immunity-related genes. First, we found a strong relationship between SLC11A1 genotype and brucellosis status (Fig. 1). Prevalence decreased with the number of copies of the less frequent allele "A330". This effect was particularly marked for individuals carrying two copies of this allele (homozygous individuals) (AES = −2.36 [−4.60, −0.13]) but less clear for heterozygous individuals (AES = −0.24 [−0.86, 0.39]), which suggests that "A330" is recessive. We also revealed a significant association with Tlr1 genotype (Fig. 2, Table S5): homozygous individuals carrying two copies of the frequent haplotype (Tlr1a) had a lower probability to be seropositive than heterozygous individuals (AES = −1.90 [−3.52, −0.29]). Tlr1 heterozygous status was also associated with ibex brucellosis status (Table S5) but this most likely result from the systematic presence of the Tlr1b haplotype in Tlr1 heterozygous ibex. We did not find any effect of Tlr2 and Tlr4 haplotypes (directional selection) or of the heterozygous/ homozygous status (heterozygote advantage) (Table S5).

Discussion
The aim of this study was to assess the neutral and immunity-related genetic diversity of the Alpine ibex population of the Bargy massif, wildlife reservoir of a persistent and virulent outbreak of brucellosis. As expected, this population established a few decades ago from a very limited number of founders (<15 ibex) shows very low genetic variation at neutral microsatellite loci and exhibits a single MHC Class II haplotype similarly to most ibex populations of the Western Alps 40 . By contrast, other immunity-related genes involved in the innate immunity such as Tlr genes (Tlr1, Tlr2 and Tlr4) and Slc11A1 have maintained polymorphism despite serial founder events occurring throughout the Alpine ibex reintroduction history. We revealed significant associations between the polymorphism of some of these genes (Tlr1 and Slc11A1) and Brucella infection status, which may suggest a role of pathogen-mediated selection in their recent evolution.   17 ; H E = 0.58-0.67 in red deer Cervus elaphus, Zachos et al. 41 ). However, these differences may arise because the microsatellite loci used in the different species underwent different ascertainment bias. Our results are consistent with the levels previously observed in other Alpine ibex populations reintroduced in Switzerland during the XX th century 42 and in particular, for the Mont Pleureur population (H E = 0.42) from which the Bargy founders came from. Low genetic diversity might reduce individual fitness and population demographic performance due to reduced evolutionary potential 43,44 . For example, in the Gran Paradiso ibex population, there is a direct relationship between multi-locus microsatellite heterozygosity, male body mass and gastrointestinal nematode infestation (fecal egg counts) (i.e. Heterozygosity-fitness correlation) 35 . In our case, we did not observe any relationship between the genome-wide heterozygosity (MLH) of ibex and their Brucella infection status. This is in agreement with the recent study of Brambilla et al. 45 that did not find any effect of neutral heterozygosity on infectious keratoconjunctivitis symptoms, another bacterial disease recurrently affecting Alpine ibex populations. Similar low levels of neutral genetic diversity were observed in other brucellosis-infected ibex populations (Grand Paradiso) in which the prevalence level has remained low and infection vanished spontaneously 26 . Therefore, contrary to our hypothesis, the low level of neutral diversity at Bargy is in itself not sufficient to explain its potential higher sensitivity to brucellosis compared to other ibex populations. However, a growing number of studies illustrated the poor predictive power of microsatellite loci to evaluate functional diversity and showed the importance of accounting for immunogenetic variation to get a better understanding of the population's adaptive response potential to pathogen-mediated selection 46,47 .

Maintenance of immunity-related gene polymorphism in Alpine ibex.
Our study is the first to our knowledge to reveal immunity-related polymorphism in non-MHC genes in Alpine ibex. Previous work on ibex focused exclusively on MHC Class II genes 45,48 that accounts for only a fraction of the genetic variation in pathogen resistance 5,9 . MHC class II molecules principally bind exogenous antigens and are primarily involved in the immune response to extracellular pathogens 49 . By contrast, genes encoding toll-like receptors (Tlrs) play www.nature.com/scientificreports www.nature.com/scientificreports/ a fundamental role in vertebrate innate immune defense against intracellular and extracellular micropathogens including viruses, bacteria, protozoa or fungi 50 . In particular, these genes are involved in the recognition and immune response to pathogenic bacteria that recurrently cause outbreaks such as Mycoplasma sp. (causing pneumonia or infectious kerotoconjunctivitis infections) 51,52 or Brucella abortus 53,54 . Contrary to our expectation and despite serial founder events throughout the species reintroduction history, the three studied Tlr genes and the SlC11A1 have all maintained several haplotypes while the MHC Class II DRB (second exon) is monomorphic in Bargy. There are very few studies that investigate Tlr polymorphism in ungulates populations: levels of Tlr haplotype diversity in Bargy ibex (N H = 3-5, H D ) appears similar to those reported for the same genes in roe deer populations (N H = 3-5) 17 , domestic goat (N H = 1-4) 55 or cattle (N H = 2-5) 56 . Tlr4 showed four functional haplotypes (i.e. encoding different amino-acid sequences) including three alleles with moderate to high frequencies (>25%). This balanced polymorphism which contrasts with the very low neutral and MHC diversity in Bargy, suggests that balancing pathogen-mediated selection may have favoured the maintenance of genetic variation as showed in other free-ranging animal species 12,57,58 . Furthermore, it is worth noting that the presence of balanced polymorphism does not mean that the gene is involved in brucellosis resistance, which can be demonstrated only by studying associations between immunogenetic patterns and pathogen prevalence 59 .

Genetic factors of resistance to Brucella melitensis infection in Alpine ibex. In contrast to
Brambilla et al. 35 , we did not find evidence of correlation between multi-locus heterozygosity at neutral microsatellites and brucellosis infection status ("general effect" hypothesis). This means that a higher individual multi-locus heterozygosity (a proxy of the level of inbreeding) does not confer a resistance advantage regarding Brucella infection. In contrast, we found associations between single immunity-related genes (here Tlr1, Slc11A1) and ibex infection status hence supporting the "local effect" hypothesis in agreement with Brambilla et al. 45 . These authors observed a lower susceptibility of MHC Drb heterozygous to infectious keratoconjunctivitis (heterozygote advantage hypothesis) while we revealed here an effect of specific alleles. Ibex carrying two copies of the less frequent allele of the Slc11A1 gene were less likely to develop the Brucella infection, which suggests that this "resistance allele" is recessive. Associations between a specific Slc11A1 allele and Brucella infection has been already reported in domestic goat 36 and water buffalo 31 . Moreover, in an in-vitro experiment on water buffalo's monocytes infected by several Brucella species including B. melitensis, Borriello et al. 60 demonstrated that some Slc11A1 genetic variants confer a higher Slc11A1 mRNA expression and a higher ability in controlling the intracellular replication of the Brucella. The "resistant allele" (A330) showed similar frequency in the Bargy and Aravis populations (25% and 21% respectively) but is much more frequent in other Alpine populations (e.g. 54% in Belledonne massif (northern French Alps), 40% in Grand Paradiso National Park -results not shown). Further comparison with other recently reintroduced ibex populations should be performed to investigate whether the particularly low frequency of this resistant variant in Bargy may have favored the emergence of this brucellosis outbreak together with other, non-genetic factors such as spatial ibex distribution and abundance.
Our results also suggest that the most frequent haplotype of Tlr1 (96%) conferred a resistance advantage against Brucella infection. However, this resistance advantage is relative: although 92% of ibex in Bargy are homozygous for this allele, 37% are still seropositive. Therefore, individuals that carry out this allele have a lower probability to get infected after a contact with the bacteria, but do not show a complete resistance. Such relationship between Tlr1 genotype and Brucella infection has been already evidence in Cattle by Prakash et al. 32 . Tlr1 interact with Brucella sp. by recognizing specific pathogen-associated molecular patterns. Individuals carrying the rare allele (Tlr1b) may have a decreased ability to recognize Brucella and thus to trigger an adequate immune response. Further work is needed to clarify the functional role of these immunity-related genes against brucellosis in Alpine ibex and in particular how these polymorphisms affects the detection and the replication of the bacteria (see Borriello et al. 60 ). Particularly, in vitro studies of ibex white blood cells and their response to Brucella 36 would help to confirm a variation linked to immune genes and to identity cellular pathways to resistance.
In agreement with our hypothesis and the literature on domestic ungulates, we showed significant associations between immunity-related genes (specific alleles Slc11A1 and Tlr1) and brucellosis infection status in Alpine ibex. However, this does not in itself provide evidence for directional pathogen-mediated selection because we did not actually demonstrate that these alleles confer a selective advantage. Seropositive ibex were euthanized, so there is no data available on the impact of B. melitensis infection on survival and reproduction rates in this population. However, brucellosis is well-known to cause abortion and reproductive failures in other domestic and wild ungulates 37,38 . So we can reasonably hypothesize that B. melitensis infection actually affects the fitness of ibex in Bargy and may ultimately lead to selection on these genes.

Management implications.
Genetic monitoring is widely recognized as a valuable tool for the management and conservation of populations 61 . In particular, access to population genetic parameters may help to better understand the ecology of host/pathogen interaction and thus better manage wildlife diseases 62,63 . Most genetic studies used microsatellites and MHC Class II markers to investigate population adaptive diversity and tolerance or resistance to diseases 8 with very few exceptions 64 . However, Alpine ibex exhibit a very low neutral genetic diversity and a unique MHC Class II DRB exon 2 haplotype in most recently introduced populations. Furthermore, most potentially zoonotic diseases affecting wild ungulate populations and notably Alpine ibex (e.g. Q fever, tuberculosis, chlamydiosis) involve micropathogens (virus, bacteria, protozoa) and MHC Class II often plays a minor role in the immune response against these pathogens 49 . Our study highlights the need to look beyond MHC class II and to explore the diversity of other candidate immunity-related genes 5 , in particular MHC class I, cytokine and Toll-like receptors genes that are centrally involved in the innate immune response against micro pathogens in humans and domestic species. (2020) 10:4400 | https://doi.org/10.1038/s41598-020-61299-2 www.nature.com/scientificreports www.nature.com/scientificreports/ In a conservation framework, our results call for immediate management actions to increase the neutral and immunity-related diversity of the Bargy ibex population, whose census size had been halved (~300 individuals in 2014-2017 against ~600 in 2013) over the monitoring period. Translocations of unrelated individuals from the Gran Paradiso source population or from other genetic units of Alpine ibex 34,65 is recommended to prevent inbreeding depression and increase both genome-wide 35,45 and immunity-related diversity 47 . Further, these results raise the question of the interplay between sanitary management and population genetic structure. Here, the test-and-cull management procedure was expected to affect the genetic structure of the population, by selectively eliminating infected individuals carrying specific alleles. In the same line, some authors suggested that artificial selection of specific Slc11A1 or Tlr genotypes in breeding programs may help to increase population natural resistance and limit spread of brucellosis in domestic ungulates 31,32 . However, it is important to keep in mind that natural populations are exposed to multi-pathogen pressures (e.g. brucellosis and Q fever in the Bargy population) and that most immunity-related genes are involved in the immune response against several pathogens. Therefore, artificially selecting a specific genotype to increase population resistance to a pathogen may lead to the loss of key haplotypes in the host-control of other diseases, particularly in species with very low genetic diversity. Incorporating data on immunogenetic variation in epidemiological models should help to predict how variation in allele frequency may affect the spread dynamics of the disease, and how disease management may affect genetic structure of the population.  Microsatellite and immunity-related loci genotyping. DNA was extracted from blood using the DNeasy Blood and Tissue kit (QIAGEN). All individuals were genotyped at 25 polymorphic putatively neutral microsatellites 42 (see list and primers in Table S1). We carried out genotyping replicates for about 25% of the samples to ensure the reliability of the genotyping. Micro-Checker 2.2.3 69 was employed to assess the frequency of null alleles and scoring errors (stuttering or large allelic dropout). The same individuals were also typed at a microsatellite located at the 3′UTR of the Solute Carrier Family 11 member A1 (SLC11A1) gene using the following primers: Fw-TCTGGACCTGTCTCATCACC and Rv-ACTCCCTCTCCATCTTGCTG 70 . The SLC11A1 gene, formerly known as natural resistance-associated macrophage protein 1 (NRAMP1) gene, is involved in the innate resistance to intra-cellular pathogens 71 . Polymorphism in microsatellites of the 3′ UTR of the gene were associated with resistance to Brucella infection in water buffalo Bubalus bubalis 36 and in the domestic goat 31 , a species closely related to the Alpine ibex. In particular, some SLC11A1 genotypes have a higher ability of controlling the intracellular replication of Brucella in vitro 60 .

Methods
We also genotyped Mhc-Drb and Tlr genes for a subset of the individuals (N = 146) sampled between 2012 and 2014. The second exon of the Mhc-Drb class II gene encoding the ligand-binding domain of the protein was amplified and sequenced using Illumina MiSeq system following the procedure detailed in Quéméré et al. 17 . Haplotypes and individual genotypes were identified using the SESAME barcode software 72 . Tlr genes (Tlr1, Tlr2 and Tlr4) were genotyped using the two-step procedure described in Quéméré et al. 17 . A pilot study on 30 individuals was first performed to identify polymorphic sites (SNPs). We screened almost the entire coding region of the three Tlr genes (92% in average) including the leucine-rich extracellular region of receptors involved in antigen-binding. Details on primer sequences and accession numbers are provided in Table S2. SNPs were then genotyped for 146 individuals using the KASPar allele-specific genotyping system provided by KBiosciences (Hoddesdon, UK). Details on SNP position and codon change can be found in Table S3. Haplotypes were then reconstructed from the phased SNPs using the procedure implemented in DNASP v5 73 . Genetic diversity. The genetic diversity of both microsatellites and immunity-related genes was evaluated by calculating the expected heterozygosity (H E ) and the number of alleles/haplotypes (N A ) as implemented in GENETIX software v4.05.2 74 . For each locus, a test for Hardy-Weinberg proportions was performed based on 1000 permutations. To allow comparison with other ibex populations 42 and between the two monitoring periods (before and after culling operations), we used FSTAT version 2.9.3 75 to calculate the allelic richness (sNA), a standardized measure of number of alleles corrected for differences in samples size. Lastly, standardized multi-locus heterozygosity at neutral microsatellites (MLH) was calculated for each individual as the ratio of its heterozygosity to the mean heterozygosity of the population 76 . Our set of microsatellite loci showed significant identity disequilibrium (g2 ± SD = 0.011 ± 0.0065, p = 0.007), (e.g. a measure of covariance in heterozygosity), a condition for detecting heterozygosity-fitness correlations. Both MLH and g 2 were calculated using the R package inbreedR 77 . To assess whether management operations (i.e. removal of 40% of the population between 2013 and 2015) had led to a decrease of the genetic diversity, we performed a comparison (using Student's t-tests) of the allelic richness (sN A ) and expected heterozygosity (H E ) of ibex sampled before (between autumn 2012 and autumn 2013) and after (years 2016-2017) the slaughtering operations that occurred between autumn 2013 and autumn 2015.
Genetic association with brucellosis infection. We analyzed the relationships between the serological status of ibex and neutral or immunity-related gene variation using generalized linear models (GLM) of the binomial family (using a logit link function). The serological status of each individual was coded as 1 (positive, considered to have been exposed and not resistant) or 0 (negative, considered to have been either non-exposed or exposed but resistant).
We fitted models testing for the effect of genetic effects, while taking into account other biological and environmental factors that possibly affect pathogen exposure. We tested the "heterozygosity fitness-correlation" hypothesis (i.e. "general effect") by fitting individual standardized multi-locus microsatellite heterozygosity (MLH) as a fixed effect. For each immunity-related gene, we tested the "heterozygote advantage hypothesis" by fitting a binary fixed effect (heterozygote vs homozygote) and the "haplotype advantage" hypothesis by considering associations between the infection status and the number of copies (0, 1 or 2) of Tlr or Slc11A1 haplotypes with frequency >10%. Concerning non-genetic factors, following Marchand et al. 29 , we included the effects of ibex Scientific RepoRtS | (2020) 10:4400 | https://doi.org/10.1038/s41598-020-61299-2 www.nature.com/scientificreports www.nature.com/scientificreports/ sex as a categorical variable and age as a continuous variable in years (linear and quadratic terms). We also fitted a "year" effect as a categorical variable to test for varying epidemiological situations across years and a "location" effect because seroprevalence has been shown to vary among socio-spatial units 29 .
We used a multimodel inference approaches to establish which explanatory variables had an effect, averaged over plausible models [78][79][80] . The most complete models used in the inference included the effects of genetic effects hypotheses, age sex, year and location. We compared the maximal model with all its sub-models (See Tables S6  and S7 for the full lists of model tested). The effect of Slc11A1 and Tlr genes were evaluated in separate analyses because of different samples sizes (262 and 146 individuals respectively). We conducted an initial exploration of our data to ascertain their distribution and spread, to identify outliers and examine relationships between variables 81 . We used variance inflation factors (VIF) to assess which explanatory variables were collinear, and should be retained in the analyses (i.e. VIF values <3) 81 . To minimize multicollinearity, the effects of the most frequent haplotype of a gene and its heterozygosity status were evaluated in separate models.
Model averaged parameter estimates are provided with their unconditional standard errors (SE) and 95% confidence intervals, after averaging model with ΔAIC < 7 (relative to the best model) as recommended by Burnham et al. 79 and using the "zero method" (i.e. a zero is substituted into models where the parameter is absent). Analyses were performed using the package lme4 82 , MUMIN v1.7.7 83 and AICcmodavg v1.25 84 in R version 3.3.3 (R Development Core Team, 2017).

Data availability
Haplotype DNA sequences: primers and Genbank accessions are in Tables S1,S2. SNP positions and characteristics are in Table S3. Genotypes of all individuals at each microsatellite and immunity-related locus are available from the corresponding author upon request to anyone who wishes to repeat our analyses or collaborate with us.