Skip to main content

Advertisement

Research article | Open | Published:

Spatial pattern of genetic diversity and selection in the MHC class II DRB of three Neotropical bat species

Abstract

Background

Although bats are natural reservoirs of many pathogens, few studies have been conducted on the genetic variation and detection of selection in major histocompatibility complex (MHC) genes. These genes are critical for resistance and susceptibility to diseases, and host–pathogen interactions are major determinants of their extensive polymorphism. Here we examined spatial patterns of diversity of the expressed MHC class II DRB gene of three sympatric Neotropical bats, Carollia perspicillata and Desmodus rotundus (Phyllostomidae), and Molossus molossus (Molossidae), all of which use the same environments (e.g., forests, edge habitats, urban areas). Comparison with neutral marker (mtDNA D-loop) diversity was performed at the same time.

Results

Twenty-three DRB alleles were identified in 19 C. perspicillata, 30 alleles in 35 D. rotundus and 20 alleles in 28 M. molossus. The occurrence of multiple DRB loci was found for the two Phyllostomidae species. The DRB polymorphism was high in all sampling sites and different signatures of positive selection were detected depending on the environment. The patterns of DRB diversity were similar to those of neutral markers for C. perspicillata and M. molossus. In contrast, these patterns were different for D. rotundus for which a geographical structure was highlighted. A heterozygote advantage was also identified for this species. No recombination or gene conversion event was found and phylogenetic relationships showed a trans-species mode of evolution in the Phyllostomids.

Conclusions

This study of MHC diversity demonstrated the strength of the environment and contrasting pathogen pressures in shaping DRB diversity. Differences between positively selected sites identified in bat species highlighted the potential role of gut microbiota in shaping immune responses. Furthermore, multiple geographic origins and/or population admixtures observed in C. perspicillata and M. molossus populations acted as an additional force in shaping DRB diversity. In contrast, DRB diversity of D. rotundus was shaped by environment rather than demographic history.

Background

Regarding adaptation to infectious diseases, genes of the major histocompatibility complex (MHC) are the most commonly studied in mammals. MHC genes play a significant role in the functionality and effectiveness of the immune response of vertebrates and are directly involved in fitness and adaptation of hosts to their pathogens [1]. Therefore, studying MHC gene diversity in non-model animals, known to be reservoirs of pathogens, is a critical tool to assess its evolutionary process and implication in (1) the variation of individual fitness, (2) population viability and (3) evolutionary potential in a changing environment [2]. Studies have essentially focused on the genetic variability of exon 2 of the MHC class II DR beta (DRB) gene given that this exon encodes the peptide-binding region (PBR) [1]. Hence, most of the variability displayed is recorded on this exon. Variations within the PBR define the repertoire of antigenic determinants and subsequently the ability to recognize a variable number of circulating pathogens [3]. Optimum resistance to pathogens is, therefore, the result of all intrinsic and extrinsic factors that influence the genetic variability of individuals’ immune system [4]. Most of the vertebrate populations studied so far exhibited high levels of MHC diversity both in the number of alleles and in their heterozygosities and allelic variation [58]. Pathogen-driven balancing selection and sexual selection seem to play a fundamental role in the preservation of this high level of polymorphism [2]. Pathogen-driven balancing selection acts through antagonistic host–parasite coevolution via several mechanisms: overdominant and frequency-dependent selection [9], as well as spatial and temporal variation in host pathogens [10]. Sexual selection pressures encompass mechanisms such as maternal–fetal interactions [11] and mate selection [12]. However, when assessing the genetic variability of both neutral markers and MHC genes, studies have highlighted the role of past demographic processes (e.g., fragmentation, bottlenecks, geographic isolation) in shaping the pattern of MHC variability that can sometimes surpass that of natural selection [1318].

Local immunogenetic adaptation of hosts that live in different environments was associated with different parasite and pathogen pressures [1922]. Indeed, differences in the diversity of pathogens (inducing different selection pressures on the hosts) are directly linked to environmental components. These latter (e.g., vegetation cover and density, landscape fragmentation, human occupation) modulate parasite and pathogen species richness, their survival and adaptability, as well as their distribution, transmission, developmental success and their ability to induce diseases [23]. Environmental components likewise impact the richness, population dynamics, immunocompetence and nutritional status of host species, all of which subsequently determine resistance or susceptibility to disease [24, 25]. A strong correlation was also shown between host and parasite species richness, their life history and ecological traits [2628]. Moreover, anthropogenic alterations of habitats induce changes in host–pathogen–environment interactions and are consequently linked to the emergence of infectious zoonotic diseases [2931]. Therefore, considering the role of the environment is critical for the assessment of MHC gene variability.

Bats (Chiroptera) are the second largest species-rich mammalian order, accounting for 20 % of all mammal species in the world. Compared to other mammals, their uncommon biological features (e.g., ability to fly, long lifespan, complex social structures) and their long-term evolutionary history with pathogens make them excellent reservoirs and spreaders of viruses (e.g., rabies virus, henipaviruses, coronaviruses), which have a significant impact on both human and animal health [32]. However, bats rarely exhibit clinical or pathological signs of diseases and appear to coexist with their pathogens in a disease-free state, characteristic of reservoir hosts [33]. At the order level, the monophyly of Chiroptera was evidenced via phylogenetic relationships inferred from intron sequences of DRB [34]. The monophyletic origin of DRB genes was also evidenced at the family level via phylogenetic relationships inferred from DRB sequences from Saccopteryx bilineata, Myotis spp., Noctilio spp. and Carollia perspicillata [35]. Other studies, investigating the diversity of MHC DRB, showed significant differences in the polymorphism of MHC genes between species. This polymorphism was essentially influenced by (1) a pathogen-driven selection [36], (2) a social structure driven by MHC-mediated post-copulatory mechanisms [37], (3) diversifying selection and recombination events [35, 38] and (4) geographical constraints resulting in spatial variation of pathogen-mediated selection and enhanced susceptibility to environmental changes [39, 40].

The Amazon, a major biodiversity hotspot in South America, possesses a large diversity of bat species and pathogens in a wide variety of climates and vegetation formations [41, 42]. French Guiana, a tropical Amazonian region near the equatorial zone, harbors a great diversity of bats, with 103 species registered [43]. In this region, the impact of deforestation on the composition and dynamics of bat communities was assessed [44, 45] as well as, the role of habitats on rabies virus circulation and maintenance [46].

In this study, we analyzed the genetic variability of the expressed MHC class II DRB exon 2 in three bat species: two Phyllostomidae, Carollia perspicillata and Desmodus rotundus, and one Molossidae, Molossus molossus. These three widely distributed species have ecological plasticity and tolerance to disturbance. In French Guiana, these species are sympatric and use the same habitats (e.g., forests, edge habitats, and urban areas).

The key hypothesis of this study was that composition and distribution of MHC DRB alleles were specific to the environments (forests vs. disturbed areas), rather than randomly distributed in space. Consequently, we should observe local immunogenetic adaptation to the contrasting pathogen pressures or equally adapted alleles. To assess which are the best factors that predict the MHC diversity, pathogen-mediated selection, recombination, gene conversion, demographic history and population structure were investigated.

There is a higher diversity of microorganisms in forest environments, compared to disturbed environnments, due to greater host species richness and better transmission-promoting parameters [47, 48]. For this reason, we expect higher levels of MHC diversity in forest environments facing lower disturbance pressures, where higher parasite and pathogen diversities imply a higher selection pressure. Furthermore, assuming that bats using the same roosting area and/or the same foraging areas would be subjected to similar pathogen pressures, we should observe similar trends in intra- and inter-specific MHC diversity. In contrast, once the demographic neutral genetic histories—which may also influence MHC diversity—are controlled, we should observe differences in selective histories between bats inhabiting different environments.

To identify different signatures of selection in the DRB exon 2, which would imply area-specific recognition capabilities, conformation of the HLA-DRB1 was used as a reference to detect putative sites directly involved in antigen binding (ABS) and consequently subjected to selection. Additionally, to reduce the misidentification bias of ABS, the selection signature within the PBR was investigated without a priori identifying species-specific ABS for further comparisons.

MHC genes exhibit high levels of allele similarity within species as well as between related species and the occurrence of identical MHC alleles in related species is frequent. Convergence and trans-species polymorphism are thought to be responsible for this trans-species evolution. To highlight which of these two mechanisms acts predominantly on the evolutionary history of the DRB gene in the three species investigated, phylogenetic relationships were inferred from the DRB sequences obtained here and with other available chiropteran sequences. Finally, MHC spatial diversity was compared to that of neutral markers (mtDNA D-loop) to highlight the impact of demographic processes and population structure on the diversity pattern in the three bat species investigated.

Methods

Ethics statement

Animals were captured, handled, sampled and, whenever necessary, euthanized following ASM guidelines [49] under the supervision of researchers who had been granted the French “Expérimentation animale niveau 1” diploma. Bats are not protected by law in French Guiana; however, the project was submitted to the Conseil Scientifique Régional pour le Patrimoine Naturel de la Guyane and approved. Captures that occurred within protected areas (nature reserves) received approval by the Conseil Scientifique Régional du Patrimoine Naturel on 26 January 2010 and ad-hoc authorizations (no. 2011-35 dated 05/30/2011, no. 35 and 59 obtained on 03/21/2013 and 04/17/2013, respectively, and delivered by the Prefecture of French Guiana).

Study areas, sample collection, DNA and RNA extraction

The three bat species (C. perspicillata, D. rotundus and M. molossus) were sampled at 14 sites. Three sites were located in pristine lowland forests (PF sites) with low disturbance pressure, four sites were located in edge habitats (EH sites), five in anthropized areas (AA sites), and two in urban and periurban areas (UR sites) (Fig. 1, Additional file 1: Table S1 and Additional file 2: Figure S1).

Fig. 1
figure1

Map of sampling sites of all individuals used in this study across four environments. For clarity, nearby collecting sites within 15 km were grouped. Sites are numbered and labeled according to the type of environment to which they correspond: edge habitats (dark square), anthropized areas (light circle), pristine primary lowland forests (light square), urban and periurban areas (dark circle). Pie charts indicate the proportion of the most frequent alleles (alleles with the highest overall frequencies: Cape-DRB*16 for C. perspicillata (a); Dero-DRB*02 and 19 for D. rotundus (b); Momo-DRB*08 for M. molossus (c)), other shared alleles and private alleles (those detected only in one individual) in each collecting sites/environments, for each species: C. perspicillata (a), D. rotundus (b) and M. molossus (c)

Bats were trapped with mist nets erected near breeding sites, roosts, at forest edges, around livestock, or through putative foraging courses. Bat species were identified in the field using external morphology and, prior to release, a brachial vein puncture and wing biopsies were performed. A few individuals captured were euthanized at the laboratory to collect organs. A total of 45 C. perspicillata, 59 D. rotundus and 42 M. molossus were sampled.

Nucleic acids were extracted from blood, liver or spleen depending on sampling. Extractions were performed using the NucliSENS easyMAG bio-robot (bioMérieux®).

Amplification and genotyping of the MHC class II DRB gene

All amplifications of MHC DRB genes for the three species were performed using cDNA obtained from liver and spleen RNAs. All cDNAs were synthesized using SuperScript® III Reverse Transcriptase (Invitrogen) and random hexamers (Roche) following the manufacturer’s recommendations. Expressed MHC DRB genes of the three bat species were amplified by a minimum of two independent polymerase chain reactions (PCRs) using a standard procedure with the AmpliTaq Gold DNA Polymerase PCR kit (Thermo Fisher Scientific). To optimize the detection of different alleles, a minimum of two distinct PCRs per individual were performed. Species-specific primers were designed for C. perspicillata and D. rotundus to avoid occurrence of non-amplifying alleles (Table 1). To amplify the M. molossus DRB gene, primers designed by Richman et al. [39] were modified. Primers located in conserved parts of exon 1 (JSex1.5, [37]; EX1Fd, modified from [39]), exon 2 (JS2Cape, [37]), exon 3 (EX3Rd, modified from [39]) and exon 4 (L729, [50]) were used to amplify functional MHC class II DRB alleles from the cDNA of three bats from each species, C. perspicillata, D. rotundus and M. molossus. The sequences obtained were used to design specific primers complementary to conserved parts of exon 1 (BatCMHIIF (F) and BatCMHIIF.1 (F1), Table 1) and exon 4 (BatCMHIIR (R) and BatCMHIIR.1 (R1), Table 1). The primer combinations F-R and F1-R1 were used as the first PCR for D. rotundus and C. perspicillata (Additional file 3: Figure S2). For the second PCR, all of the newly designed primers were used together with JSex1.5, EX1Rd, EX3Rd and L729 in different combinations to screen the cDNA of a total of 24 C. perspicillata and 39 D. rotundus (Additional file 3: Figure S2). Two different primer combinations with JS2Cape, EX3Rd, R and R1 were used to amplify the DRB alleles from the cDNA of 40 M. molossus (Additional file 3: Figure S2). The PCR products were cloned using the pCR™4-TOPO® TA Cloning® Kit (Invitrogen). Ten clones per PCR amplification were sent for sequencing to Beckman Coulter Genomics (Takeley, UK) using T7 and T3 primers.

Table 1 PCR primers used to amplify MHC class II DR beta loci in the three species investigated: C. perspicillata, D. rotundus and M. molossus

Amplification and sequencing of the mtDNA D-loop

All amplifications of mitochondrial DNA (mtDNA) control region (D-loop) for the three species were performed using genomic DNA. The primers F(mt) and P(mt) [51] were used to amplify the hypervariable domain HVI of mtDNA D-loop. PCRs were performed using a standard procedure with BIOTAQ™ DNA Polymerase PCR kit (Bioline). Sequencing was carried out by Beckman Coulter Genomics (Takeley, UK) using F(mt) and P(mt) primers.

Sampling analysis

To assess the environmental impact on the MHC class II DRB diversity, the data were analyzed according to the location of the capture sites and the type of environment at these sites. For D. rotundus and M. molossus, each capture site corresponded to one environment type. D. rotundus sampling sites corresponded to two different roots (Cave F (CF) and Cave M (CM)), considered as PF sites, and one foraging area (Saut Athanase (SA)), an EH site (Fig. 1, Additional file 1: Table S1). M. molossus individuals were captured in one foraging area (Cacao (CC)), an AA site, and two different roosts located at the EH site with different disturbance levels: Paracou (PA) with a high level of disturbance and the SA site with a low level of disturbance (Fig. 1, Additional file 1: Table S1). C. perspicillata individuals were captured in the four types of environments PF, EH, AA and UR (Fig. 1, Additional file 1: Table S1). To explore the potential bias induced by the spatial structure of collecting sites in C. perspicillata samples, spatial structure was compared with environment structure.

Data analysis

Sequences were edited using Geneious R6 [52]. DNA sequences were aligned using the Mafft alignment tool [53] included in the software. Sequences were confirmed as MHC class II DRB and mtDNA D-loop by homology analysis using the NCBI BLAST search [54]. A DRB clone sequence was regarded as valid if the following criteria were met: (1) incidence in at least two independent PCRs either from the same individual or different ones, or amplification by two different primer pairs, and (2) identified by at least three identical clones [55]. Haplotypes were reconstructed for each bat species using the DnaSP 5 [56]. As proposed by [57], the nomenclature of MHC alleles for non-human species was used: to each allele a prefix (the first two letters of the genus and the species names) was given, with the serial number attached as follows: Cape-DRB for C. perspicillata, Dero-DRB for D. rotundus and Momo-DRB for M. molossus. Cape-DRB alleles identified in this study were named after the 15 alleles previously described by [35].

Sequence diversity

Allele frequencies (F), proportion of segregating sites (S), nucleotide diversity (π), observed (H O ) and expected (H E ) heterozygosities, and gene diversity (H D ) were estimated using Arlequin 3.5 [58]. Allelic richness (R), with a correction for sample size using a rarefaction method and inbreeding coefficients (FIS) were calculated using Fstat 2.9. [59]. The mean number of differences between nucleotide and amino acid alleles were counted using MEGA 6.06 [60]. Analyses were performed per species for the entire dataset, and for each environment subgroup separately.

Gene conversion and recombination events

Evidence of genetic recombination or gene conversion events between DRB sequences was assessed using GENECONV 1.81a [61]. This program detects lengthy patterns shared between sequences despite the existence of a pronounced polymorphism. P-values were determined from both global and pairwise permutation tests, with 10,000 replicates. No mismatches between fragments were accepted, and global p-values were corrected for multiple comparisons (gscale = 0). Recombination breakpoints were identified using the GARD and SBP algorithms [62] implemented in the HyPhy package [63], available in the standard Datamonkey analysis library (http://www.datamonkey.org/, Delport et al. [64]).

Codon-based analysis of positive selection

Putative antigen-binding sites (ABS) were identified by comparison with the human ABS of the HLA-DRB1 [65]. The relative rates of dN and dS substitutions were calculated for putative ABS, non-ABS and all codons following the method of Nei and Gojobori with the Jukes and Cantor correction for multiple hits [66, 67]. A Z-test was performed to assess the deviation of the dN/dS (ω) ratio from 1.0. These tests were performed using MEGA 6.06.

Codons potentially subjected to positive selection were further identified with no a priori assumption of ABS. For this purpose, two maximum-likelihood (ML) frameworks were carried out as proposed by [68]. First, two alternative models comprised in the CodeML Program included in the PAML package version 4.7 were executed [69]. The first one (M7) allows codons to evolve either neutrally (ω = 1) or under purifying selection (ω < 1), while the second one (M8) considers a class of sites under positive selection (ω > 1). These models were compared with a likelihood ratio test (LTR; 2Δl = 2-l1-l0) with two degrees of freedom (α = 0.05) [70]. As previously described, each analysis was run twice, with starting ω values of 0.5 and 1.5, to ensure convergence [68]. The F3 × 4 model of codon frequencies was assumed for all analyses. The amino acids positively selected were identified using the Bayes empirical Bayes approach (BEB), as recommended by Yang [71], with the cutoff posterior probability set at 90 %. Secondly, five ML methods proposed in the HyPhy package were used. For all analyses, the best fitting nucleotide substitution model was assessed with the automatic model selection tool available on the server. A single likelihood ancestor counting (SLAC) model was used to detect evidence of the non-neutral evolution of the DRB gene [72]. Based on ancestral reconstruction, this model counts dN and dS changes at each codon position in a phylogeny. Fixed (FEL) and random effects likelihood (REL) models were performed to estimate the ω ratio on a site-by-site basis [72]. For the FEL model, ω ratios were estimated with an a priori distribution of rates across sites, while this distribution was determined by the present data for the REL model. An internal FEL (IFEL) model was carried out to highlight selection pressure at the population level (i.e., along internal branches). Finally, a mixed effects model of evolution (MEME) was performed to highlight both diversifying and importantly episodic selection at individual sites [73]. Codon positions were regarded as candidates for positive selection if the following two criteria were met: (1) p-values < 0.1 for SLAC, FEL, IFEL and MEME and (2) Bayes Factor > 50 for REL. Accordingly, only sites identified as being under positive selection by at least two approaches were considered [68]. Analyses were performed per species for the entire dataset and for each environment subgroup separately.

Genetic differentiation and past population dynamics

For both neutral and functional markers, the genetic differentiation (F ST) among environment subgroups and collecting sites was calculated through an implementation of analysis of molecular variance (AMOVA; 10,000 permutations) using Arlequin 3.5.

Isolation-by-distance (IBD) was examined with the web service IDBWS v.3.23 (http://ibdws.sdsu.edu/~ibdws/; [74]) by assessing the statistical significance of the correlation in a Mantel test with 30,000 permutations. The quantity F ST/(1 – F ST) was used as a genetic distance according to Rousset [75]. Geographic distances between populations were determined from GPS coordinates recorded at each collecting site.

To add a measurement of relationships between alleles at the intraspecific level, haplotype networks were constructed for each species using the parsimony statistical algorithm TCS implemented in PopArt 1.7 [76].

A Bayesian skyline plot was used to investigate changes in effective population size (Ne) over time with Beast 1.8.1 [77], using mtDNA D-loop data only. Markov Chain Monte Carlo samples were based on 2 × 107 generations, logging every 200 steps, with the first 2,000,000 generations discarded as the burn-in. The best-fitted model HKY + G was used with a relaxed molecular clock. The number of stepwise changes in Ne was set at five.

MHC class II DRB phylogeny

Phylogenetic relationships and trans-species polymorphism between DRB alleles were inferred using a Bayesian inference approach implemented in MrBayes 3 [78]. The HLA-DRB1 sequence (accession number: NG_029921) was used as outgroup. The GTR + G + I model was selected as the best-fitted model of nucleotide substitution using jModelTest 2 [79, 80] under corrected Akaike information criteria (AICc). The program was run 10 × 107 generations, with a sampling frequency of 500 and a 25 % burn-in. Validation of the inference was assessed based on the standard deviation of split frequencies, less than the expected threshold value of 0.01.

Results

DRB sequence characterization

One hundred and three bats (24 C. perspicillata, 39 D. rotundus and 40 M. molossus) were tested for the amplification of MHC DRB using cDNA. PCR product sizes ranged from 357 bp to 600 bp, for both C. perspicillata and D. rotundus and covered all of exon 2 (267 bp in length) (Additional file 3: Figure S2). For M. molossus, we successfully amplified 244 bp, excluding primers, of exon 2. We successfully obtained 674 MHC DRB-like fragments from 97 bats (148 MHC DRB-like fragments for 22 C. perspicillata, 311 for 36 D. rotundus and 215 for 39 M. molossus, Additional file 1: Table S2). Based on BLAST searches, sequences showed homologies to only DRB loci.

Considering validation criteria, 480 sequences were validated, corresponding to 82 bats (Additional file 1: Table S2). One hundred and ten sequences were validated for 19 C. perspicillata, corresponding to 23 Cape-DRB alleles (previously reported Cape-DRB*05, GenBank accession number: JQ388834, and Cape-DRB*16-37; GenBank accession numbers: KU896612–KU896633; Additional file 1: Table S3). Two hundred and thirty-three sequences were validated for 35 D. rotundus, corresponding to 30 Dero-DRB alleles (Dero-DRB*0130; GenBank accession numbers: KU896562–KU896591; Additional file 1: Table S4). One hundred and forty-seven sequences were validated for 28 M. molossus, corresponding to 20 Momo-DRB alleles (Momo-DRB*0120; GenBank accession numbers: KU896592–KU896611; Additional file 1: Table S5). All amino acid sequences were unique.

A maximum of three transcribed alleles per individual was detected in one C. perspicillata and two D. rotundus individuals (Additional file 1: Tables S3 and S4). No more than two transcribed alleles were identified in M. molossus samples (Additional file 1: Table S5).

All sequences showed BLAST homology with the HLA-DRB1 locus (accession number: NG_029921), with a maximum of 85 % nucleotide identity and 75 % amino acid identity. The percentage of nucleotide and amino acid identity with other bat species (S. bilineata, M. davidii, C. perspicillata and N. albiventris) was above 87 % and 75 % for all transcripts, respectively.

Detection of indels

Two alleles (Cape-DRB*20 and 32) presented a 3-bp deletion at position 217–219 of the nucleotide alignment (Fig. 2). These single-codon indels did not induce a frame-shifting mutation or a stop codon. Bats carrying these alleles were heterozygotes and found in EH and UR environments. Two other alleles (Dero-DRB*05 and Dero-DRB*21) presented a 6-bp insertion at position 115–120 of the alignment. These two-codon indels – one serine (S) and one asparagine (D) – did not alter the reading frame. Dero-DRB*05 was carried by four bats, all captured in SA, with one homozygote. Dero-DRB*21 was carried by one heterozygote bat captured in CF. No indel event was detected for M. molossus.

Fig. 2
figure2

Amino acid sequence variation and overall frequency of MHC class II DR beta exon 2 alleles of C. perspicillata, D. rotundus and M. molossus. DRB sequences of Homo sapiens (accession number: NG_29921), C. perspicillata (Cape-DRB*01, accession number: JQ388830), N. albiventris (Noal-DRB*01, accession number: HM347941) and A. jamaicensis (Arja-DRB*01, accession number: KJ010995) are given for comparison. Antigen-binding sites (ABS) of the HLA-DRB1 molecule are shadowed [65]. Dots mark identity with the top sequence. Numbers in italics indicate the amino acid positions within the beta 1 domain of the HLA-DRB1 molecule. Asterisks indicate the positively selected sites (PSS) identified in each species, according to acceptance criteria described in the Methods section [68]. † indicates PSS identified by [35], and ‡ indicates PSS identified by [38]

mtDNA D-loop sequence characterization

We amplified the first hypervariable segment (HVI) of the mitochondrial DNA (mtDNA) control region (D-loop) for the 146 bats tested. We obtained a 291-bp alignment for M. molossus, a 340-bp alignment for C. perspicillata and a 359-bp alignment for D. rotundus. Thirty-nine mtDNA D-loop haplotypes were identified for C. perspicillata (45 individuals, Genbank accession numbers: KU896735–KU896779), 29 for M. molossus (42 individuals, Genbank accession numbers: KU896693–KU896734) and only five for D. rotundus despite the high number of sequenced individuals (59 individuals, Genbank accession numbers: KU896634–KU896692) (Table 2).

Table 2 The MHC class II DR beta (exon 2) allelic diversity in 19 C. perspicillata, 35 D. rotundus and 28 M. molossus

DRB allelic diversity

Overall, C. perspicillata presented the highest values of allelic richness (R = 23.000), number of segregating sites (S = 114), nucleotide diversity (π = 0.140 ± 0.070), gene diversity (H D  = 0.986) and mean nucleotide and amino acid differences (37.387 ± 3.154 and 25.249 ± 2.477, respectively), while M. molossus exhibited the lowest values (Table 2). The highest observed heterozygosity (H O ) was recorded in D. rotundus and, as for the two other species, H O was lower than expected heterozygosities (H E ). Considering environment subgroups, the highest value of R was detected in UR for C. perspicillata, (R = 6.576), in CF for D. rotundus (R = 7.0), and in PA for M. molossus (R = 7.375). S ranged from 51 in CC for M. molossus to 97 in EH for C. perspicillata. The highest values of π were found in PF for C. perspicillata, CF for D. rotundus and CC for M. molossus (Table 2). FIS values were high in all species as well as in all the subgroups per species, with values ranging from 0.330 (p <0.001) for D. rotundus to 0.777 (p <0.001) for M. molossus.

Recombination tests did not reveal fragments involved in gene conversion, recombination events or breakpoints among individuals from the same species or within a given environment or between species.

mtDNA D-loop allelic diversity

C. perspicillata had the highest allelic richness values (R = 38.873), number of segregating sites (S = 69), gene diversity (H D  = 0.993) and mean nucleotide differences: 16.439 ± 7.460 (Table 2). Nucleotide diversity, for each species, was low (<0.05 in C. perspicillata and M. molossus, <0.01 in D. rotundus). Considering environment subgroups, π values ranged from 0.001 ± 0.001 in CF (D. rotundus) to 0.051 ± 0.027 in AA (C. perspicillata). The highest R values were detected in SA (M. molossus and D. rotundus), AA and EH (C. perspicillata). S ranged from one in CF (D. rotundus) to 58 in AA (C. perspicillata), and H D ranged from 0.200 in CF (D. rotundus) to 1.0 in AA and EH (C. perspicillata).

Distribution of DRB alleles across the different environments

C. perspicillata

Three alleles (Cape-DRB*16, 17, and 25) were shared either between bats trapped in the same environment (EH) or between bats trapped in distinct environments (EH vs. PF). These three alleles were the most common, with an overall frequency (F ± SD) of 0.111 ± 0.062 for Cape-DRB*16 and 0.074 ± 0.051 for both Cape-DRB*17 and Cape-DRB*25 (Fig. 2). Cape-DRB*16 was shared between one bat trapped in PF (n = 3) and two bats trapped in EH. Cape-DRB*17 was shared between two bats trapped in EH (n = 11). Nine individuals in this region were homozygotes (Additional file 1: Table S3). Cape-DRB*25 was shared between two bats trapped in PF and EH. Alleles found in UR were private alleles and three of the five trapped bats were heterozygotes.

D. rotundus

Up to 14 alleles were shared either between bats roosting in the same area or between bats foraging in the same area. Out of these alleles, one (Dero-DRB*19) was shared between six bats from the three sampling sites (CF, CM, and SA sites) (Additional file 1: Table S3). Five alleles (Dero-DRB*13, 17, 18, 27 and 30) were shared between bats from SA and CM sites, and two alleles (Dero-DRB*02 and 15) between the SA and CF sites. Dero-DRB*02 and 19 had the highest overall frequency of 0.085 ± 0.037 (Fig. 2). Four alleles (Dero-DRB*05, 06, 14 and 29) were exclusively shared between individuals trapped in SA. Two alleles (Dero-DRB*24 and 26) were solely shared between individuals trapped in CM. We did not find shared alleles between bats trapped in CF (Additional file 1: Table S4). Sixteen of the 23 bats captured in SA were heterozygotes, five in CM (n = 8) and three in CF (n = 4).

M. molossus

Up to ten shared alleles were identified but none was shared between the three sampling sites (Additional file 1: Table S5). Three alleles (Momo-DRB*02, 04 and 19) were shared between bats captured in the SA and PA sites. Two alleles (Momo-DRB*09 and 10) were shared between individuals from the SA and CC sites. One allele (Momo-DRB*07) was shared between two bats captured in the CC and PA sites. Three alleles (Momo-DRB*08, 11 and 20) were exclusively shared within SA. Momo-DRB*08 had the highest overall frequency: 0.147 ± 0.062 (Fig. 2). Four of the five bats trapped in CC were homozygotes, seven in PA (n = 8), and 11 in SA (n = 15).

Detection of selection signature

An excess of dN substitutions was detected for the three species (dN/dS = [2.049–3.548], Z = [2.600–5.619], p = [<10-5-0.01] – Table 3). More specifically, in putative antigen binding sites (ABS), dN/dS values were highly significant (p < 10-3) and ranged from 4.081 for D. rotundus to 4.939 for C. perspicillata. We found significant dN substitutions also occurring in putative non-ABS for C. perspicillata (dN/dS = 2.792, Z = 3.283, p = 0.0013) but none in both D. rotundus and M. molossus samples (p > 0.40).

Table 3 Non-synonymous (dN) and synonymous (dS) substitutions (± standard deviation) for C. perspicillata, D. rotundus and M. molossus

PAML analyses (performed by species and within subgroups per species) showed significant support for M7 and M8 (log likelihood estimates (LTR), 2Δl = [63.985–123.085], df = 2, p < 0.0001), indicating an effect of positive selection acting on the DRB gene. Overall, 22 codons appeared to be under positive selection for C. perspicillata, eight for D. rotundus, and six for M. molossus (BEB posterior probability p ≥ 90 %, SLAC, FEL, IFEL and MEME p-values < 0.1, and REL Bayes Factor > 50; Table 4). Three positively selected sites (PSS 52, 65 and 66; Fig. 2) were common to the three species investigated. These codons were the only common between D. rotundus and M. molossus. Six of the eight PSS identified in D. rotundus were also found in C. perspicillata, and all of the PSS identified in M. molossus were identified in C. perspicillata.

Table 4 Identification of species-specific positively selected sites (PSS) by maximum likelihood (ML) analysis

Regarding environment subgroups for all species, a maximum of 13 PSS was identified in UR for C. perspicillata, eight in CF for D. rotundus and 11 in SA for M. molossus. Three PSS were found in all environment subgroups for D. rotundus, four in C. perspicillata and four in M. molossus. A maximum of six unique PSS (sites detected as positive only within the defined subgroup) was observed for C. perspicillata, four for M. molossus and two unique PSS were observed in the two roosts (CF and CM) for D. rotundus (Table 4). None of the PSS detected in C. perspicillata was shared between UR and PF. The same result was obtained between SA and CM for D. rotundus and between PA and CC or SA for M. molossus. Thirteen of the 22 PSS identified were located within the human ABS region for C. perspicillata, seven for D. rotundus and five for M. molossus (Fig. 2). All the PSS identified in M. molossus were described as PSS for N. albiventris, 13 for C. perspicillata and six for D. rotundus. Six PSS in C. perspicillata were also found in A. jamaicensis, three for M. molossus and one for D. rotundus. The sites where PSS differed were located within one to six amino acid residues (3–18 bp) of the human ABS or other bat species PSS, depending on the bat species investigated.

DRB population structure

We found no significant genetic differentiation between DRB alleles of the environment subgroups defined for each species (C. perspicillata: FST range: [(−0.068) to 0.042]; D. rotundus: FST range: [(−0.047) to 0.012]; M. molossus: FST range: [0.013–0.055], Additional file 1: Table S6(A) − 8). Taking into account the collecting sites, pairwise differentiation tests did not reveal significant genetic differentiation between DRB alleles (C. perspicillata: FST range: [(−0.048) to 1.000]; D. rotundus: FST range: [(−0.047) to 0.012]; M. molossus: FST range: [0.013–0.055], (Additional file 1: Table S6(B) − 8). Minimum spanning trees showed complex connections between the haplotypes, marked by numerous nucleotide mutations, without clustering of haplotypes either by site or by environment (Figs. 3, 4, 5 and 6). We did not find relevant isolation-by-distance in our samples (Additional file 1: Table S9).

Fig. 3
figure3

Haplotype networks of MHC class II DR beta exon 2 alleles and of the mitochondrial DNA (mtDNA) control region (D-loop) at the nucleotide level for C. perspicillata. Nodes are proportional to the number of bats carrying each haplotype, colored by the environments where the bats were trapped (see legend). Hatch marks represent mutations. Interruptions in lines indicate the presence of more than ten mutations

Fig. 4
figure4

Haplotype networks of MHC class II DR beta exon 2 alleles and of the mitochondrial DNA (mtDNA) control region (D-loop) at the nucleotide level for C. perspicillata based on collecting sites. Nodes are proportional to the number of bats carrying each haplotype, colored by the capture site where the bats were trapped (see legend). Hatch marks represent mutations

Fig. 5
figure5

Haplotype networks of MHC class II DR beta exon 2 alleles and the mitochondrial DNA (mtDNA) control region (D-loop) at the nucleotide level for D. rotundus. Nodes are proportional to the number of bats carrying each haplotype, colored by the capture site where the bats were trapped (see legend). Hatch marks represent mutations. Interruptions in lines indicate the presence of more than ten mutations

Fig. 6
figure6

Haplotype networks of MHC class II DR beta exon 2 alleles and of the mitochondrial DNA (mtDNA) control region (D-loop) at the nucleotide level for M. molossus. Nodes are proportional to the number of bats carrying each haplotype, colored by the capture site where the bats were trapped (see legend). Hatch marks represent mutations. Interruptions in lines indicate the presence of more than ten mutations

mtDNA D-loop population structure

We found significant genetic differentiation between some populations of D. rotundus (CF vs. CM FST: 0.475, p-value <0.0001; CM vs. SA FST: 0.617, p-value < 0.0001; but CM vs. SA FST: 0.016, p-value not significant). No significant genetic differentiation was observed between populations for the two other species (C. perspicillata FST range: [(−0.029) to (−0.016)]; M. molossus FST range: [(−0.040) to 0.017], Additional file 1: Table S6(A)–8). Pairwise computations performed between collecting sites did not reveal significant genetic differentiation (C. perspicillata FST range: [0.072 to 1.000]; M. molossus FST range: [(−0.040) to 0.017], Additional file 1: Table S6(B)–8). The haplotype networks for C. perspicillata revealed two main groups of haplotypes that possessed all the haplotypes found in the different environments or collecting sites (Figs. 3 and 4). The same result was obtained for M. molossus, with a maximum of four mutations between haplotypes (Fig. 5). The lowest number of mutations was found for D. rotundus for which the five haplotypes were linked by one connection, marked by only one or two mutations (Fig. 6). No population expansion or collapse was observed for the three species. We did not find relevant isolation-by-distance in our samples (Additional file 1: Table S7).

DRB phylogenetic relationships

The analysis included the 73 DRB sequences identified in this study (174 bp; nucleotide position 48–219), as well as the 63 published chiropteran DRB sequences selected randomly from A. jamaicensis [38], 15 from C. perspicillata [35], 15 from Myotis spp. (M. velifer and M. vivesi [39]), 28 from Noctilio spp. (N. albiventris and N. leporinus [35, 37]) and 17 from S. bilineata [35, 36]. The phylogenetic tree showed six major clades (Fig. 7). We observed a clustering by species or genus supported by medium and high posterior probabilities – C. perspicillata (0.61), D. rotundus (0.67), M. molossus (1), Noctilio spp. (0.50) – with an intermingled clustering of N. albiventris and N. leporinus DRB alleles. Within the C. perspicillata clade, DRB sequences clustered with the previously published sequences, with a high posterior probability (0.85). One DRB sequence of C. perspicillata was related to a D. rotundus clade with a posterior probability of 0.79. One DRB sequence from S. bilineata clustered with the clade Noctilio spp., with a low posterior probability value. The 63 DRB sequences of the A. jamaicensis species clustered with five DRB sequences of D. rotundus and two DRB sequences of C. perspicillata, with a posterior probability of 0.70. Within this phyllostomid clade, we observed a species-dependent clustering with high posterior probabilities. Two DRB sequences of C. perspicillata clustered with a low posterior probability with two sequences from A. jamaicensis. We did not observe a characteristic clustering for S. bilineata and Myotis spp. However, one DRB sequence from M. velifer clustered together with six DRB sequences from S. bilineata with a low posterior probability.

Fig. 7
figure7

Phylogenetic relationships of the MHC class II DR beta sequences based on part of exon 2 (174 bp) with values of posterior probabilities for nodes. Species designation follows the GenBank accession numbers for previously described DRB sequences. Arja A. jamaicensis, Cape C. perspicillata, Dero D. rotundus, Momo M. molossus, Myve M. velifer, Myvi M. vivesi, Noal N. albiventris, Nole N. leporinus, Sabi S. bilineata. HLA-DRB1 was used as the outgroup. TSP indicates trans-species polymorphism of the newly described DRB alleles with previously described DRB alleles

Discussion

In this study, we explored the genetic variability of the expressed MHC DRB genes of three sympatric Neotropical bat species, looking for selection signatures within the PBR and investigating the role of the environment and the population structure on MHC diversity.

Diversity patterns in the MHC class II DRB gene

In the three species, identified alleles (23 in C. perspicillata, 30 in D. rotundus and 20 in M. molossus) were unique in amino acid sequences, suggesting a non-redundant nucleotide polymorphism. The number of alleles detected per species and per environment subgroup within each species indicated a relatively high level of MHC DRB variation. These results are in agreement with previously reported allelic copy number variation between or within bat species [1, 3540].

This high DRB polymorphism was observed in the three species whatever the degree of habitat disturbance (heavily disturbed, slightly impacted or pristine). In addition, no loss of polymorphism induced by low host and/or pathogen diversity, expected to occur in disturbed environments [47, 48] was observed. Furthermore, mtDNA D-loop sequences, considered as neutral, revealed high levels of genetic diversity for C. perspicillata and M. molossus. High levels of neutral diversity have been associated with either multiple geographic origin and/or population admixture [81, 82], suggesting that C. perspicillata and M. molossus populations originated from one of these two phenomena or from both. The high level of DRB polymorphism can be related to population history for these two species with the effects of additional selection factors such as pathogen-mediated selection.

A high number of DRB homozygotes was also found in M. molossus in all habitats as well as in C. perspicillata but in edge habitats only. For these two species, the presence of a specific allele (allele-specific overdominance) ensuring an adequate immunity can be hypothesized [1]. Considering D. rotundus, a high number of heterozygotes for DRB was observed in all habitats. This result suggests a heterozygote advantage related to the forest (pristine and edge habitats) environment known to encompass higher species and pathogen richness [83, 84].

The most frequent MHC DRB alleles encountered in each bat species (Dero-DRB*02, Dero-DRB*15, Dero-DRB*19, Cape-DRB*16 and Momo-DRB*08) were all identified in forest environments. These alleles are likely ancient and play a significant role in the immune response, conferring a selective advantage to the bats in sympatry in these habitats. Indeed, given that forest environments are characterized by a high level of pathogen diversity, these alleles most likely possess enhanced recognition capabilities and might be able to recognize a broader range of pathogens. They thus allow a more efficient immune response in line with the ecological diversity. Moreover, the geographical structure highlighted in D. rotundus, using mtDNA D-loop results, did not affect the allelic distribution. Indeed, most of the alleles identified in D. rotundus were shared between sites, and one allele (Dero-DRB*19) was present at all sites. In contrast, no geographical structure was found in C. perspicillata even if the alleles identified in urban and periurban areas were rare and not shared with other groups. These results suggest that environment rather than population structure drives the allelic distribution, corroborating a local adaptation hypothesis. Although edge and forest habitats may differ in pathogenic pressures (see above), the results observed for C. perspicillata may be related to the high tolerance of this species to habitat disturbance and to its dominance in bat communities in edge habitats [45], likely allowing continuous exposure to the communities of pathogens of both pristine and disturbed forests, and pathogen transfer between these habitats.

Duplication detected in Phyllostomidae

Gene conversion and recombination events were not detected for any species. The identification of three expressed DRB alleles in individuals from C. perspicillata and D. rotundus supported the existence of at least two functional DRB loci. However, underestimation of the number of alleles and expressed loci can be suspected, considering the limited number of recombinant clones per primer combinations and the use of cDNA only allowing the detection of the most frequently expressed alleles. The occurrence of multiple MHC DRB loci varying in the number of loci between individuals and species was reported in Neotropical bat species such as C. perspicillata (three DRB loci), S. bilineata (ten DRB loci) [35], A. jamaicensis (three DRB loci) [38], as well as in other animals [1, 7, 85, 86]. Together with environmental and biological patterns, gene conversion and recombination processes are known to drive the extensive MHC polymorphism [18, 87].

Variable gene duplication between closely related taxa and among individuals from the same species is characteristic of MHC genes and plays a critical role in the adaptive evolution of the host [88]. Gene duplication can occur over different time scales and therefore predates or follows speciation events. Furthermore, the absence of recombination between alleles of different species was reported. This result was regarded as evidence of recent duplication events between loci that occurred after speciation [35]. Our findings support this hypothesis since some recent radiation events occurred in Neotropical bat species during the Pleistocene era (0.01–1.8 Mya) [8991].

Despite the indels detected in some C. perspicillata and D. rotundus sequences, the reading frames remained unaltered, suggesting that these sequences encode functional proteins (Fig. 2). Nevertheless, the heterozygote-dominant character of individuals with these alleles suggests a turnover of these loci generated by evolutionary processes such as the silencing of duplicate genes by mutation or deletion [88]. Under the birth-and-death model, genes created by duplication can be maintained for long periods, deleted or transformed as pseudogenes and, therefore, contribute to host fitness and adaptability [92]. The phenomenon is especially significant among heterozygote D. rotundus, for which one of the two alleles is among the most widespread in the population studied (Dero-DRB*02, 05 and 15), and supports a shift towards an allele with better recognition capabilities.

Evidence of historical positive selection

Higher rates of non-synonymous (dN) vs. synonymous (dS) nucleotide substitutions were observed, especially in sites supposedly involved in antigen-binding. This excess of dN over dS in ABS is consistent with balancing selection acting on DRB loci, as a sign of historical positive selection in polymorphic MHC genes [93]. Furthermore, an unexpected excess of dN in putative non-ABS for C. perspicillata suggests strong positive selection acting on these sites and highlights species-specific divergence on sites involved in antigen recognition.

A high congruence was observed between species-specific positively selected sites (PSS) in the three species and human ABS demonstrating the functional homology of these sites. In contrast, PSS not identified as human ABS might play a role in MHC DRB recognition capability or in the molecule stability for the species studied. Here, nine of the PSS were detected in C. perspicillata, while only one was detected in both D. rotundus and M. molossus. Therefore, the observed excess of dN in C. perspicillata was attributed to PSS identified as non-ABS in humans. This result demonstrates the critical role of positive selection in shaping MHC diversity and corroborates that positive selection is driven by a species-specific immune response. Focusing on C. perspicillata, comparisons of PSS identified here with those previously described revealed high congruence, even among PSS not identified as human ABS [35]. A similar result was observed between C. perspicillata and A. jamaicensis, two Phyllostomidae with the same diet [38]. Taken together, these findings reveal family and species-specific selection pressures acting on MHC genes. Moreover, the high similarity of PSS between A. jamaicensis and C. perspicillata highlights that bat species sharing similar feeding strategies might be subjected – independently of their habitats – to similar pathogenic pressures, mainly from their microbiota [94]. Comparisons between each environment subgroup, for each present case, revealed the existence of private PSS depending on the area, suggesting a habitat-specific selection pressure process. Indeed, while six private PSS were counted in urban and periurban areas for C. perspicillata, only two were detected in edge habitats and one in pristine forest.

Phylogeny, trans-species polymorphism and demographic process

Two mechanisms are thought to be responsible for the trans-specific similarity of MHC genes: convergent evolution [93] and trans-species polymorphism [95, 96]. Evidence of trans-species polymorphism in bats was reported in both Myotis spp. and Noctilio spp. [35] and at the family level between the two phyllostomids A. jamaicensis and C. perspicillata [38]. Phylogenetic analyses in our study did not reveal any characteristic clustering of MHC alleles depending on the habitat, but rather clades at the family, genus and species levels (Fig. 3). These results suggest either similar pathogenic pressures or that these alleles possess a larger antigenic recognition capability that constitutes a major asset to the immune response of these species. Intermingled clusterings observed in this survey (between DRB sequences from A. jamaicensis, C. perspicillata and D. rotundus) are in agreement with trans-species polymorphism reported in bats [35, 38]. The observed trans-species polymorphism for Myotis spp. [35] was not detected in the present study. The results highlighted clustering between S. bilineata and M. velifer. However, this clustering was no longer observed by analyzing full-length exon 2. Observation of trans-species polymorphism within bat families (e.g., Phyllostomidae) strengthens the scenario of independent modes of evolution of MHC DRB alleles, allowing balancing selection to retain substantial allelic lineages before speciation events [95, 96]. Extended periods of host–pathogen coevolution is thought to contribute to trans-species polymorphism given that host-sharing pathogens and exposure to similar pathogen pressures are believed to induce similar changes in different host species [97].

Although pathogen-driven selection is known to be crucial in MHC diversity [98], demographic processes such as genetic drift, bottlenecks, expansions, fragmentation and geographic isolation have been demonstrated to also participate in the shaping of this extensive polymorphism [15, 17, 99, 100]. Evolutionary analyses of MHC DRB alleles for the three species did not reveal any area-dependent clustering effect. Furthermore, analyses of partial sequences of the mtDNA D-loop revealed no structuration for C. perspicillata and M. molossus, whereas it was observed for D. rotundus. Demographic analysis performed for the three species did not detect any substantial bottleneck or expansion. While the patterns observed for the first two species corroborate the hypothesis of a similarity between the intrinsic genetic diversity of the species and that of MHC, the lack of similarity between these patterns of diversity for D. rotundus contradicts this statement. These results suggest that the signature of an area-limited pathogen-driven selection, supported by the presence of private PSS rather than by demographic processes such as distance isolation, is observed in this case.

Conclusions

This study was the first to investigate MHC DRB polymorphism in three sympatric bat species, M. molossus, C. perspicillata and D. rotundus, from the Amazonian region. Our results revealed a high genetic variability in the MHC DRB gene. Natural selection, as well as local adaptation driven by different parasite and pathogen exposures across environments, contribute to the maintenance of an extensive MHC polymorphism. The richness of pathogens in the different habitats should be investigated to strengthen our assumptions on potential local adaptation. Further studies using finer neutral markers such as microsatellites will be necessary to detect possible confounding effects, such as hidden structuring and demographic history, and thus to better understand the drivers of MHC allelic diversity.

Abbreviations

ABS:

Antigen-binding sites

AICc:

Corrected Akaike information criteria

BEB:

Bayes empirical Bayes

DRB:

Exon 2 of the MHC class II DR beta

FEL:

Fixed effects likelihood

IFEL:

Internal fixed effects likelihood

MEME:

Mixed effects model of evolution

MHC:

Major histocompatibility complex

ML:

Maximum-likelihood

mtDNA D-Loop:

Mitochondrial DNA control region

PBR:

Peptide-binding region

PCR:

Polymerase chain reaction

PSS:

Species-specific positively selected sites

REL:

Random effects likelihood

SLAC:

Single likelihood ancestor counting

References

  1. 1.

    Sommer S. The importance of immune gene variability (MHC) in evolutionary ecology and conservation. Front Zool. 2005;2:16.

  2. 2.

    Piertney SB, Oliver MK. The evolutionary ecology of the major histocompatibility complex. Heredity. 2006;96:7–21.

  3. 3.

    Kloch A, Babik W, Bajer A, SińSki E, Radwan J. Effects of an MHC-DRB genotype and allele number on the load of gut parasites in the bank vole Myodes glareolus. Mol Ecol. 2010;19:255–65.

  4. 4.

    Penn DJ, Ilmonen P. Major Histocompatibility Complex (MHC). In: John Wiley & Sons, Ltd, editor. Encycl. Life Sci. [Internet]. Chichester, UK: John Wiley & Sons, Ltd; 2005 [cited 2015 May 20]. Available from: http://doi.wiley.com/10.1038/npg.els.0003986

  5. 5.

    Meyer-Lucht Y, Otten C, Püttker T, Sommer S. Selection, diversity and evolutionary patterns of the MHC class II DAB in free-ranging Neotropical marsupials. BMC Genet. 2008;9:39.

  6. 6.

    Castro-Prieto A, Wachter B, Sommer S. Cheetah Paradigm Revisited: MHC diversity in the world’s largest free-ranging population. Mol Biol Evol. 2011;28:1455–68.

  7. 7.

    Rivero-de Aguilar J, Schut E, Merino S, Martínez J, Komdeur J, Westerdahl H. MHC class II B diversity in blue tits: a preliminary study. Ecol Evol. 2013;3:1878–89.

  8. 8.

    Yao G, Zhu Y, Wan Q-H, Fang S-G. Major histocompatibility complex class II genetic variation in forest musk deer (Moschus berezovskii) in China. Anim Genet. 2015;46:535–43.

  9. 9.

    Hughes AL, Yeager M. Natural selection at major histocompatibility complex loci of vertebrates. Annu Rev Genet. 1998;32:415–35.

  10. 10.

    Hedrick PW. Pathogen resistance and genetic variation at MHC loci. Evolution. 2002;56:1902–8.

  11. 11.

    Hedrick PW. HLA-sharing, recurrent spontaneous abortion, and the genetic hypothesis. Genetics. 1988;119:199–204.

  12. 12.

    Ditchkoff SS, Lochmiller RL, Masters RE, Hoofer SR, Bussche RA. Major-histocompatibility-complex-associated variation in secondary sexual traits of white-tailed deer (Odocoileus virginianus): evidence for good-genes advertisement. Evolution. 2001;55:616–25.

  13. 13.

    Oosterhout C, Joyce DA, Cummings SM, Blais J, Barson NJ, Ramnarine IW, et al. Balancing selection, random genetic drift, and genetic variation at the major histocompatibility complex in two wild populations of guppies (Poecilia reticulata). Evolution. 2006;60:2562–74.

  14. 14.

    Babik W, Pabijan M, Arntzen J w, Cogalniceanu D, Durka W, Radwan J. Long-term survival of a urodele amphibian despite depleted major histocompatibility complex variation. Mol Ecol. 2009;18:769–81.

  15. 15.

    Alcaide M. On the relative roles of selection and genetic drift in shaping MHC variation. Mol Ecol. 2010;19:3842–4.

  16. 16.

    Marsden CD, Woodroffe R, Mills MGL, McNUTT JW, Creel S, Groom R, et al. Spatial and temporal patterns of neutral and adaptive genetic variation in the endangered African wild dog (Lycaon pictus): Spatial and temporal diversity in wild dogs. Mol Ecol. 2012;21:1379–93.

  17. 17.

    Zeisset I, Beebee TJC. Drift rather than selection dominates MHC class II allelic diversity patterns at the biogeographical range scale in natterjack toads Bufo calamita. PLoS ONE. 2014;9:e100176.

  18. 18.

    Yu F-J, Zhu Y, Xiong T-Y, Wan Q-H, Zhang H-M. Balancing selection and recombination drive genetic variation at MHC class I genes in the giant panda. Sci Bull. 2015;60:136–8.

  19. 19.

    Dionne M, Miller KM, Dodson JJ, Caron F, Bernatchez L. Clinal variation in MHC diversity with temperature: evidence for the role of host-pathogen interaction on local adaptation in Atlantic salmon. Evol Int J Org Evol. 2007;61:2154–64.

  20. 20.

    Lenz TL, Eizaguirre C, Rotter B, Kalbe M, Milinski M. Exploring local immunological adaptation of two stickleback ecotypes by experimental infection and transcriptome-wide digital gene expression analysis. Mol Ecol. 2013;22:774–86.

  21. 21.

    Froeschke G, Sommer S. Role of selection versus neutral processes determining genetic variation in a small mammal along a climatic gradient in southern Africa. Evol Ecol. 2014;28:1169–90.

  22. 22.

    Clozato CL, Mazzoni CJ, Moraes-Barros N, Morgante JS, Sommer S. Spatial pattern of adaptive and neutral genetic diversity across different biomes in the lesser anteater (Tamandua tetradactyla). Ecol Evol. 2015;5:4932–48.

  23. 23.

    Froeschke G, Matthee S. Landscape characteristics influence helminth infestations in a peri-domestic rodent-implications for possible zoonotic disease. Parasit Vectors. 2014;7:393.

  24. 24.

    Fuller CA, Postava-Davignon MA, West A, Rosengaus RB. Environmental conditions and their impact on immunocompetence and pathogen susceptibility of the Caribbean termite Nasutitermes acajutlae. Ecol Entomol. 2011;36:459–70.

  25. 25.

    Hayman DTS, Pulliam JRC, Marshall JC, Cryan PM, Webb CT. Environment, host, and fungal traits predict continental-scale white-nose syndrome in bats. Sci Adv. 2016;2:e1500831.

  26. 26.

    Kamiya T, O’Dwyer K, Nakagawa S, Poulin R. Host diversity drives parasite diversity: meta-analytical insights into patterns and causal mechanisms. Ecography. 2014;37:689–97.

  27. 27.

    Poulin R. Parasite biodiversity revisited: frontiers and constraints. Int J Parasitol. 2014;44:581–9.

  28. 28.

    Zukal J, Bandouchova H, Bartonicka T, Berkova H, Brack V, Brichta J, et al. White-nose syndrome fungus: a generalist pathogen of hibernating bats. PLoS One. 2014;9:e97224.

  29. 29.

    Johnson JS, Reeder DM, McMichael III JW, Meierhofer MB, Stern DWF, Lumadue SS, et al. Host, pathogen, and environmental characteristics predict white-nose syndrome mortality in captive little brown myotis (Myotis lucifugus). PLoS One. 2014;9:e112502.

  30. 30.

    Zarlenga DS, Hoberg E, Rosenthal B, Mattiucci S, Nascetti G. Anthropogenics: human influence on global and genetic homogenization of parasite populations. J Parasitol. 2014;100:756–72.

  31. 31.

    Becker DJ, Streicker DG, Altizer S. Linking anthropogenic resources to wildlife-pathogen dynamics: a review and meta-analysis. Ecol Lett. 2015;18:483–95.

  32. 32.

    Chomel BB, Stuckey MJ, Boulouis H-J, Aguilar- Setién A. Bat-Related Zoonoses. In: Sing A, editor. Zoonoses - Infect. Affect. Hum. Anim. [Internet]. Dordrecht: Springer Netherlands; 2015 [cited 2015 May 20]. p. 697–714. Available from: http://link.springer.com/10.1007/978-94-017-9457-2_28

  33. 33.

    Brook CE, Dobson AP. Bats as “special” reservoirs for emerging zoonotic pathogens. Trends Microbiol. 2015;23:172–80.

  34. 34.

    Kupfermann H, Satta Y, Takahata N, Tichy H, Klein J. Evolution of Mhc–DRB introns: implications for the origin of primates. J Mol Evol. 1999;48:663–74.

  35. 35.

    Schad J, Voigt CC, Greiner S, Dechmann DKN, Sommer S. Independent evolution of functional MHC class II DRB genes in New World bat species. Immunogenetics. 2012;64:535–47.

  36. 36.

    Mayer F, Brunner A. Non-neutral evolution of the major histocompatibility complex class II gene DRB1 in the sac-winged bat Saccopteryx bilineata. Heredity. 2007;99:257–64.

  37. 37.

    Schad J, Dechmann DK, Voigt CC, Sommer S. MHC class II DRB diversity, selection pattern and population structure in a neotropical bat species, Noctilio albiventris. Heredity. 2011;107:115–26.

  38. 38.

    Real-Monroy MD, Martínez-Méndez N, Ortega J. MHC-DRB Exon 2 Diversity of the Jamaican Fruit-Eating Bat (Artibeus jamaicensis) from Mexico. Acta Chiropterologica. 2014;16:301–14.

  39. 39.

    Richman AD, Herrera MLG, Ortega-García S, Flores-Martínez JJ, Arroyo-Cabrales J, Morales-Malacara JB. Class II DRB polymorphism and sequence diversity in two vesper bats in the genus Myotis: DRB polymorphism in Myotis. Int J Immunogenet. 2010;37:401–5.

  40. 40.

    Palmer JM, Berkman LK, Marquardt PE, Donner DM, Jusino MA, Lindner DL. Preliminary characterization of little brown bats (Myotis lucifugus) immune MHC II DRB alleles using next-generation sequencing. Peer J Prepr. 2016;4:e1662v1.

  41. 41.

    Myers N, Mittermeier RA, Mittermeier CG, da Fonseca GAB, Kent J. Biodiversity hotspots for conservation priorities. Nature. 2000;403:853–8.

  42. 42.

    Fenton MB, Simmons NB. It’s a bat! In: Bats: a world of science and mystery. University of Chicago Press; 2015 [cited 20 mai 2015].

  43. 43.

    Catzeflis F, Dewynter M, Pineau K. Liste taxonomique commentée des chiroptères de Guyane. Le Rhinolophe. 2013;19:89–102.

  44. 44.

    Brosset A, Charles-Dominique P, Cockle A, Cosson J-F, Masson D. Bat communities and deforestation in French Guiana. Can J Zool. 1996;74:1974–82.

  45. 45.

    Delaval M, Charles-Dominique P. Edge effects on frugivorous and nectarivorous bat communities in a neotropical primary forest in French Guiana. Revue d'écologie. 2006;61:343–52.

  46. 46.

    de Thoisy B, Bourhy H, Delaval M, Pontier D, Dacheux L, Darcissac E, et al. Bioecological drivers of rabies virus circulation in a Neotropical bat community. PLoS Negl Trop Dis. 2016;10:e0004378.

  47. 47.

    Jones KE, Patel NG, Levy MA, Storeygard A, Balk D, Gittleman JL, et al. Global trends in emerging infectious diseases. Nature. 2008;451:990–3.

  48. 48.

    Pedersen AB, Davies TJ. Cross-species pathogen transmission and disease emergence in primates. Ecohealth. 2010;6:496–508.

  49. 49.

    Sikes RS, Gannon WL. Guidelines of the American society of mammalogists for the use of wild mammals in research. J Mammal. 2011;92:235–53.

  50. 50.

    Bowen L, Aldridge BM, Gulland F, Van Bonn W, DeLong R, Melin S, et al. Class II multiformity generated by variable MHC-DRB region configurations in the California sea lion (Zalophus californianus). Immunogenetics. 2004;56:12–27.

  51. 51.

    Wilkinson GS, Chapman AM. Length and sequence variation in evening bat D-loop mtDNA. Genetics. 1991;128:607–17.

  52. 52.

    Kearse M, Moir R, Wilson A, Stones-Havas S, Cheung M, Sturrock S, et al. Geneious Basic: An integrated and extendable desktop software platform for the organization and analysis of sequence data. Bioinformatics. 2012;28:1647–9.

  53. 53.

    Katoh K, Misawa K, Kuma K, Miyata T. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 2002;30:3059–66.

  54. 54.

    Johnson M, Zaretskaya I, Raytselis Y, Merezhuk Y, McGinnis S, Madden TL. NCBI BLAST: a better web interface. Nucleic Acids Res. 2008;36:W5–9.

  55. 55.

    Kennedy L, Ryvar R, Gaskell R, Addie D, Willoughby K, Carter S, et al. Sequence analysis of MHC DRB alleles in domestic cats from the United Kingdom. Immunogenetics. 2002;54:348–52.

  56. 56.

    Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25:1451–2.

  57. 57.

    Klein J, Bontrop RE, Dawkins RL, Erlich HA, Gyllensten UB, Heise ER, et al. Nomenclature for the major histocompatibility complexes of different species: a proposal. Immunogenetics. 1990;31:217–9.

  58. 58.

    Excoffier L, Lischer HEL. Arlequin suite ver 3.5: a new series of programs to perform population genetics analyses under Linux and Windows. Mol Ecol Resour. 2010;10:564–7.

  59. 59.

    Goudet J. FSTAT, a program to estimate and test gene diversities and fixation indices (version 2.9.3).

  60. 60.

    Tamura K, Stecher G, Peterson D, Filipski A, Kumar S. MEGA6: molecular evolutionary genetics analysis version 6.0. Mol Biol Evol. 2013;30:2725–9.

  61. 61.

    Sawyer S. Statistical tests for detecting gene conversion. Mol Biol Evol. 1989;6:526–38.

  62. 62.

    Pond SLK, Posada D, Gravenor MB, Woelk CH, Frost SDW. Automated phylogenetic detection of recombination using a genetic algorithm. Mol Biol Evol. 2006;23:1891–901.

  63. 63.

    Pond SLK, Frost SDW, Muse SV. HyPhy: hypothesis testing using phylogenies. Bioinformatics. 2005;21:676–9.

  64. 64.

    Delport W, Poon AFY, Frost SDW, Pond SLK. Datamonkey 2010: a suite of phylogenetic analysis tools for evolutionary biology. Bioinformatics. 2010;26:2455–7.

  65. 65.

    Brown JH, Jardetzky TS, Gorga JC, Stern LJ, Urban RG, Strominger JL, et al. Three-dimensional structure of the human class II histocompatibility antigen HLA-DR1. Nature. 1993;364:33–9.

  66. 66.

    Jukes TH, Cantor CR. CHAPTER 24 - Evolution of protein molecules. In: Munro HN, editor. Mamm. Protein Metab. [Internet]. Academic Press; 1969 [cited 2015 Jul 25]. p. 21–132. Available from: http://www.sciencedirect.com/science/article/pii/B9781483232119500097

  67. 67.

    Nei M, Gojobori T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol. 1986;3:418–26.

  68. 68.

    Wlasiuk G, Nachman MW. Adaptation and constraint at toll-like receptors in primates. Mol Biol Evol. 2010;27:2172–86.

  69. 69.

    Yang Z. PAML 4: Phylogenetic Analysis by Maximum Likelihood. Mol Biol Evol. 2007;24:1586–91.

  70. 70.

    Yang Z, Nielsen R, Goldman N, Pedersen A-MK. Codon-substitution models for heterogeneous selection pressure at amino acid sites. Genetics. 2000;155:431–49.

  71. 71.

    Yang Z. Bayes empirical bayes inference of amino acid sites under positive selection. Mol Biol Evol. 2005;22:1107–18.

  72. 72.

    Pond SLK, Frost SDW. Not so different after all: a comparison of methods for detecting amino acid sites under selection. Mol Biol Evol. 2005;22:1208–22.

  73. 73.

    Murrell B, Wertheim JO, Moola S, Weighill T, Scheffler K, Kosakovsky Pond SL. Detecting individual sites subject to episodic diversifying selection. PLoS Genet. 2012;8:e1002764.

  74. 74.

    Jensen JL, Bohonak AJ, Kelley ST. Isolation by distance, web service. BMC Genet. 2005;6:13.

  75. 75.

    Rousset F. Genetic differentiation and estimation of gene flow from F-Statistics under isolation by distance. Genetics. 1997;145:1219–28.

  76. 76.

    Leigh JW, Bryant D. popart: full-feature software for haplotype network construction. Methods Ecol Evol. 2015;6(9):1110–6.

  77. 77.

    Drummond AJ, Suchard MA, Xie D, Rambaut A. Bayesian phylogenetics with BEAUti and the BEAST 1.7. Mol Biol Evol. 2012;29:1969–73.

  78. 78.

    Ronquist F, Huelsenbeck JP. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics. 2003;19:1572–4.

  79. 79.

    Guindon S, Gascuel O. A simple, fast, and accurate algorithm to estimate large phylogenies by maximum likelihood. Syst Biol. 2003;52:696–704.

  80. 80.

    Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nat Methods. 2012;9:772–72.

  81. 81.

    Consuegra S, Phillips N, Gajardo G, de Leaniz CG. Winning the invasion roulette: escapes from fish farms increase admixture and facilitate establishment of non-native rainbow trout. Evol Appl. 2011;4:660–71.

  82. 82.

    Monzón-Argüello C, Garcia de Leaniz C, Gajardo G, Consuegra S. Less can be more: loss of MHC functional diversity can reflect adaptation to novel conditions during fish invasions. Ecol Evol. 2013;3(10):3359–68.

  83. 83.

    Dunn RR, Davies TJ, Harris NC, Gavin MC. Global drivers of human pathogen richness and prevalence. Proc R Soc Lond B Biol Sci. 2010;14:rspb20100340.

  84. 84.

    Gay N, Olival KJ, Bumrungsri S, Siriaroonrat B, Bourgarel M, Morand S. Parasite and viral species richness of Southeast Asian bats: Fragmentation of area distribution matters. Int J Parasitol Parasites Wildl. 2014;3:161–70.

  85. 85.

    Winternitz JC, Wares JP. Duplication and population dynamics shape historic patterns of selection and genetic variation at the major histocompatibility complex in rodents. Ecol Evol. 2013;3:1552–68.

  86. 86.

    Pechouskova E, Dammhahn M, Brameier M, Fichtel C, Kappeler PM, Huchard E. MHC class II variation in a rare and ecological specialist mouse lemur reveals lower allelic richness and contrasting selection patterns compared to a generalist and widespread sympatric congener. Immunogenetics. 2015;67:229–45.

  87. 87.

    Schaschl H, Wandeler P, Suchentrunk F, Obexer-Ruff G, Goodman SJ. Selection and recombination drive the evolution of MHC class II DRB diversity in ungulates. Heredity. 2006;97:427–37.

  88. 88.

    Hughes AL, Yeager M. Natural selection and the evolutionary history of major histocompatibility complex loci. Front Biosci. 1998;3:d509–16.

  89. 89.

    Martins FM, Templeton AR, Pavan AC, Kohlbach BC, Morgante JS. Phylogeography of the common vampire bat (Desmodus rotundus): Marked population structure, Neotropical Pleistocene vicariance and incongruence between nuclear and mtDNA markers. BMC Evol Biol. 2009;9:294.

  90. 90.

    Pavan AC, Martins F, Santos FR, Ditchfield A, Redondo RA. Patterns of diversification in two species of short-tailed bats (Carollia Gray, 1838): the effects of historical fragmentation of Brazilian rainforests. Biol J Linn Soc. 2011;102:527–39.

  91. 91.

    Rojas D, Warsi OM, Dávalos LM. Bats (Chiroptera: Noctilionoidea) challenge a recent origin of extant neotropical diversity. Syst Biol. 2016;65:432–48.

  92. 92.

    Nei M, Gu X, Sitnikova T. Evolution by the birth-and-death process in multigene families of the vertebrate immune system. Proc Natl Acad Sci. 1997;94:7799–806.

  93. 93.

    Hughes A. Balancing selection: The major histocompatibility complex. Adaptive Evolution of Genes and Genomes. New York: Oxford University Press. 1999; pp. 54–89.

  94. 94.

    Carrillo-Araujo M, Taş N, Alcántara-Hernández RJ, Gaona O, Schondube JE, Medellín RA, et al. Phyllostomid bat microbiome composition is associated to host phylogeny and feeding strategies. Evol Genomic Microbiol. 2015;6:447.

  95. 95.

    Klein J. Origin of major histocompatibility complex polymorphism: The trans-species hypothesis. Hum Immunol. 1987;19:155–62.

  96. 96.

    Těšický M, Vinkler M. Trans-species polymorphism in immune genes: general pattern or MHC-restricted phenomenon? J Immunol Res. 2015;2015:1–10.

  97. 97.

    de Bellocq JG, Suchentrunk F, Baird SJE, Schaschl H. Evolutionary history of an MHC gene in two leporid species: characterisation of Mhc-DQA in the European brown hare and comparison with the European rabbit. Immunogenetics. 2008;61:131–44.

  98. 98.

    Prugnolle F, Manica A, Charpentier M, Guégan JF, Guernier V, Balloux F. Pathogen-driven selection and worldwide HLA Class I diversity. Curr Biol. 2005;15:1022–7.

  99. 99.

    Miller HC, Allendorf F, Daugherty CH. Genetic diversity and differentiation at MHC genes in island populations of tuatara (Sphenodon spp.). Mol Ecol. 2010;19:3894–908.

  100. 100.

    Fumagalli M, Sironi M, Pozzoli U, Ferrer-Admettla A, Pattini L, Nielsen R. Signatures of environmental genetic adaptation pinpoint pathogens as the main selective pressure through human evolution. PLoS Genet. 2011;7:e1002355.

Download references

Acknowledgments

We are grateful to François Catzeflis and Marguerite Delaval for useful discussions on the biology, ecology and dynamics of bats. All field volunteers and owners and/or managers of capture sites are warmly acknowledged for their assistance in captures.

Funding

A. Salmier was supported by a grant from European funds (PO FSE 2007-2013) and "Investissement d’Avenir" managed by Agence Nationale de la Recherche (CEBA, Ref. ANR-10- LABEX-25-01). This study was conducted within the CAROLIA program supported by European funds (ERDF/FEDER) and assistance from Région Guyane and Direction Régionale pour la Recherche et la Technologie. It also received a European Commission "REGPOT-CT-2011-285837- STRonGer" grant within the FP7 and "Investissement d’Avenir" grants managed by the Agence Nationale de la Recherche (CEBA, Ref. ANR-10- LABEX-25-01). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Availability of data and materials

Nucleotide sequences for the 73 DRB alleles and for mtDNA D-loop are available on Genbank (Accession numbers: KU896562–KU896633 and KU896634–KU896779, respectively).

Authors’ contribution

Conceived and designed the experiments: AS AL. Analyzed the data: AS AL. Contributed reagents/materials/analysis tools: AS BdT VL AL. Wrote the paper: AS BdT BCR VL AL. Field work and acquisition data: AS AL BdT. All authors read and approved the final manuscript.

Competing interests

The authors declare that they have no competing interests.

Author information

Correspondence to Anne Lavergne.

Additional files

Additional file 1: Table S1.

Supplementary methods, including the characteristics of all capture sites (Table S1), the list of samples successfully amplified for the MHC class II DRB loci (Table S2), the distribution of DRB alleles per species (Table S3-5), the pairwise differentiation computations per species (Table S6-8), as well as the correlations of genetics and geographic distance (Table S9). (XLS 179 kb)

Additional file 2: Figure S1.

Map of French Guiana showing the capture sites of C. perspicillata, D. rotundus and M. molossus. For clarity, nearby sites within 15 km were grouped. Sites are numbered and labeled according to the type of environment to which they correspond: edge habitats (dark square), anthropized areas (light circle), pristine primary lowland forests (light square), urban and periurban areas (dark circle). Pie chart indicates the proportion of bat species sampled, with C. perspicillata in light grey, D. rotundus in black and M. molossus in dark grey. Small charts indicate a total number of individuals caught ≤ 10, medium-sized charts (16-22 individuals) , large charts (34-61 individuals) . Characteristics of the different sites are given in Additional file 1: Table S1. (PDF 432 kb)

Additional file 3: Figure S2.

Positions of PCR primers used to amplify the indicated fragments of the MHC class II loci in C. perspicillata, D. rotundus and M. molossus, based on cDNA. The structure of the cDNA of the MHC class II DRB gene is based on [50]. According to each species, boxes indicate the amplified region using each couple of primers. The shaded region represents the region of interest, namely exon 2. Sequences and references of all primers used are given in Table 1. (PDF 26 kb)

Rights and permissions

Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Keywords

  • MHC
  • DRB polymorphism
  • Selection
  • Carollia perspicillata
  • Desmodus rotundus
  • Molossus molossus