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

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. Electronic supplementary material The online version of this article (doi:10.1186/s12862-016-0802-1) contains supplementary material, which is available to authorized users.


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 peptidebinding 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 [5][6][7][8]. 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 [13][14][15][16][17][18].
Local immunogenetic adaptation of hosts that live in different environments was associated with different parasite and pathogen pressures [19][20][21][22]. 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 [26][27][28]. Moreover, anthropogenic alterations of habitats induce changes in host-pathogenenvironment interactions and are consequently linked to the emergence of infectious zoonotic diseases [29][30][31]. 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.

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

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

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 (F IS ) 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 GENE-CONV 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.data monkey.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 d N and d S 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 d N /d S (ω) 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-l 1 -l 0 ) 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 d N and d S 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 a see Additional file: Figure S2, b Primer modified from Richman et al. [39] to amplify the exon 2 of MHC class II DRB gene 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 × 10 7 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 × 10 7 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.

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

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

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).  Table 3). More specifically, in putative antigen binding sites (ABS), d N /d S values were highly significant (p < 10 -3 ) and ranged from 4.081 for D. rotundus to 4.939 for C. perspicillata. We found significant d N substitutions also occurring in putative non-ABS for C. perspicillata (d N /d S = 2.792, Z = 3.283, p = 0.0013) but none in both D. rotundus and M. molossus samples (p > 0.40). 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.

Detection of selection signature
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.  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). 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

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   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 Fig. 7 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 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) [89][90][91].
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 (d N ) vs. synonymous (d S ) nucleotide substitutions were observed, especially in sites supposedly involved in antigen-binding. This excess of d N over d S 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 d N 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 speciesspecific 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 d N 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 subjectedindependently of their habitatsto 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.