# The Value of Molecular vs. Morphometric and Acoustic Information for Species Identification Using Sympatric Molossid Bats

Yann Gager, Emilia Tarland, Dietmar Lieckfeldt, Matthieu Ménage, Fidel

Botero-Castro, Stephen J Rossiter, Robert H S Kraus, Arne Ludwig, Dina K

N Dechmann

<figure><img src="image_1.png" /></figure>

► To cite this version:

Yann Gager, Emilia Tarland, Dietmar Lieckfeldt, Matthieu Ménage, Fidel Botero-Castro, et al.. The Value of Molecular vs. Morphometric and Acoustic Information for Species Identification Using Sympatric Molossid Bats. PLoS ONE, 2016, 11 (3), pp.e0150780. 10.1371/journal.pone.0150780. hal-03043588

HAL Id: hal-03043588

https://hal.umontpellier.fr/hal-03043588v1

Submitted on 7 Dec 2020

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
<figure><img src="image_3.png" /></figure>

RESEARCH ARTICLE

# The Value of Molecular vs. Morphometric and Acoustic Information for Species Identification Using Sympatric Molossid Bats

Yann Gager¹,²,³*, Emilia Tarland⁴,⁵, Dietmar Lieckfeldt⁵, Matthieu Ménage¹,², Fidel Botero-Castro⁶, Stephen J. Rossiter⁷, Robert H. S. Kraus¹,², Arne Ludwig⁵, Dina K. N. Dechmann¹,²,⁸

¹ Department of Migration and Immuno-Ecology, Max Planck Institute for Ornithology, Radolfzell, Germany,

² Department of Biology, University of Konstanz, Konstanz, Germany, ³ International Max Planck Research

School for Organismal Biology, University of Konstanz, Konstanz, Germany, ⁴ Swedish University of

Agricultural Sciences, Department of Animal Breeding and Genetics, Uppsala, Sweden, ⁵ Department of

Evolutionary Genetics, Leibniz-Institute of Zoo and Wildlife Research, Berlin, Germany, ⁶ Institut des

Sciences de l'Evolution, UMR 5554-CNRS-IRD, Université de Montpellier, Montpellier, France,

⁷ School of Biological and Chemical Sciences, Queen Mary University of London, London, United Kingdom,

⁸ Smithsonian Tropical Research Institute, Panamá, Rep. of Panamá

* yann.gager@gmail.com

<figure><img src="image_4.png" /></figure>

OPEN ACCESS

**Citation:** Gager Y, Tarland E, Lieckfeldt D, Ménage M, Botero-Castro F, Rossiter SJ, et al. (2016) The Value of Molecular vs. Morphometric and Acoustic Information for Species Identification Using Sympatric Molossid Bats. PLoS ONE 11(3): e0150780.
doi:10.1371/journal.pone.0150780

**Editor:** Danilo Russo, Università degli Studi di Napoli Federico II, ITALY

**Received:** October 16, 2015

**Accepted:** February 17, 2016

**Published:** March 4, 2016

**Copyright:** © 2016 Gager et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.

**Data Availability Statement:** All relevant data are within the paper and its Supporting Information files.

**Funding:** This work was funded by the German Research Foundation (DFG-GZ DE 1807/3-1 to DKND, http://www.dfg.de/en/index.jsp) and by the International Max Planck Research School for Organismal Biology (to YG, http://www.orn.mpg.de/2453/Short_portrait). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

## Abstract

A fundamental condition for any work with free-ranging animals is correct species identification. However, in case of bats, information on local species assemblies is frequently limited especially in regions with high biodiversity such as the Neotropics. The bat genus *Molossus* is a typical example of this, with morphologically similar species often occurring in sympatry. We used a multi-method approach based on molecular, morphometric and acoustic information collected from 962 individuals of *Molossus bondae*, *M. coibensis*, and *M. molossus* captured in Panama. We distinguished *M. bondae* based on size and pelage coloration. We identified two robust species clusters composed of *M. molossus* and *M. coibensis* based on 18 microsatellite markers but also on a more stringently determined set of four markers. Phylogenetic reconstructions using the mitochondrial gene co1 (DNA barcode) were used to diagnose these microsatellite clusters as *M. molossus* and *M. coibensis*. To differentiate species, morphological information was only reliable when forearm length and body mass were combined in a linear discriminant function (95.9% correctly identified individuals). When looking in more detail at *M. molossus* and *M. coibensis*, only four out of 13 wing parameters were informative for species differentiation, with *M. coibensis* showing lower values for hand wing area and hand wing length and higher values for wing loading. Acoustic recordings after release required categorization of calls into types, yielding only two informative subsets: approach calls and two-toned search calls. Our data emphasizes the importance of combining morphological traits and independent genetic data to inform the best choice and combination of discriminatory information used in the field. Because parameters can vary geographically, the multi-method approach may need to be adjusted to local species assemblies and populations to be entirely informative.
**Introduction**

Molecular biology, with the study of mitochondrial and nuclear genomes, has revolutionized our understanding of the distribution and evolutionary history of worldwide species diversity. In the context of mammalian species diversity, the order Chiroptera (bats) constitutes an exceptional taxon, with over 1331 species listed in a recent systematic review [1] representing a fifth of all extant mammals. Molecular studies have also led to the discovery of many cryptic lineages and boosted the number of described bat species. For example, analyses of mitochondrial genes revealed several cryptic species in well-studied areas such as Europe [2–6]. The use of DNA barcoding [7] led to a reevaluation of the number of bat species in the tropics [8–10]. Based on their sequence similarity, the barcodes can be clustered into a Molecular Operational Taxonomic Unit (MOTU) [11]. One great advantage of DNA barcoding is the important database available for comparative purposes (BOLD: The Barcode of Life Data System, [12]). However, DNA barcodes present pitfalls linked to maternal inheritance (reviewed in [13]) and should always be considered in conjunction with other sources of data. For instance, nuclear microsatellite loci were used successfully to identify *Pipistrellus kuhlii* as one biological species with two mitochondrial barcodes [14]. The use of nuclear microsatellites is also powerful to detect potential interspecific hybridization, otherwise undetected via the sole use of mitochondrial barcodes [15]. Other taxonomic parameters, such as morphological characters or echolocation calls, should also be combined with molecular data, following, for example, the framework of Integrated Operational Taxonomic Units (IOTUs) [16]. Integrating traditional taxonomy to molecular taxonomy is seen as the future of taxonomy [17].

Despite this recent boost of bat diversity with molecular species identification, the status of many bat taxa is not yet firmly established. The bat genus *Molossus* (family Molossidae; E. Geoffroy, 1805) is a typical example of this. These Neotropical bats occur from Northern Mexico to Southern Argentina. A systematic review from 1913 described a total of 19 species [18]. Many of these species were later synonymized and seven or eight species, depending on the authors, were recognized in the latest taxonomic reviews [19–22]. In addition, one species, *M. alvarezi*, was newly described based on size, pelage coloration and morphological characteristics [23]. Despite broad agreement among systematic reviews, the taxonomic boundaries and names within the genus are not settled. For example, *M. bondae* (J.A. Allen, 1904) and *M. currentium* (O. Thomas, 1901) can be grouped under the name *M. currentium* [19] or considered as two species based on their distribution in Central or South America [22]. Similarly, *M. molossus* (Pallas, 1766) has been described as being “desperately in need of revision” [19] and probably represents a species complex; indeed, *M. coibensis* (J. A. Allen, 1904) was treated as a synonym of *M. molossus* [24,25] yet is now considered a full species based on recent systematic assessments [19,21,22].

To date, few studies have applied molecular information to address questions regarding the taxonomy of the genus *Molossus*. The first molecular investigation of the evolutionary relationships within the genus relied on allozymes [21]. A more recent study identified only higher-level relationships between genera of the family Molossidae using one mitochondrial gene and three nuclear genes [26]. Here we compare molecular data from DNA-based markers with more commonly used morphometric and bioacoustic information to assess the reliability of each type of information for the identification of several *Molossus* species in Panama. We distinguished the *Molossus* species at our study site with a set of newly developed microsatellite markers, sequence data from the mitochondrial gene *co1* (DNA barcode), and the mitochondrial region d-loop for *M. molossus* and matched them with common field identification methods, i.e. morphological measurements and echolocation call recordings. While we were able to identify the molossid species at our site in Panama with our methods, we also find that one or
even several field-based methods may not be sufficient for the proper identification of morphologically similar species whose traits may locally vary quite substantially.

# Materials and Methods

## Ethics statements

Capture and handling of animals was carried out with permission from the Autoridad Nacional del Ambiente in Panama with approval from the Institutional Animal Care and Use Committee of the Smithsonian Tropical Research Institute (2012-0505-2015). All animals were gently handled during measurements of morphological parameters, photographs of wings, genetic sampling and acoustic recording. All animals were released back in clearings in the same area in which they were captured. Heart tissue for genetic marker development came from a freshly deceased bat found in a private home.

## Sampling and data acquisition

During different fieldwork seasons between 2008 and 2013, we captured a total of 962 bats of the genus *Molossus* in Panama. Of these, 935 individuals were captured from various buildings in the village of Gamboa (Panama, 09°07′ N, 79°41′ W), 21 from the roof of the Smithsonian Tropical Research Institute's (STRI) laboratory building on Barro Colorado Island (09°09′ N, 79°50′ W) as well as a dead tree off the shore of BCI, and seven from the roof of STRI's dormitory at the Bocas del Toro research station (09°21′ N, 82°15′ W). We used mist-nets (Ecotone, Gdynia, Poland) to catch bats during their evening emergence. We determined sex, age, forearm length, body mass, reproductive status and marked each individual with unique subcutaneous passive integrated transponder (Trovan ID-100, Euro ID, Weilerswist, Germany). We also sampled wing membrane tissue using a biopsy punch (2 or 3 mm, Stiefel, U.S.A.) for genotyping purposes [27]. During some fieldwork seasons, we also collected wing photos and echolocation calls for some individuals. We selected data only for individuals that were genotyped later for microsatellites. We retained size-referenced wing photos for the 116 genotyped bats to obtain measurements for several wing parameters (see below for details on wing morphology evaluation). Finally, we selected echolocation calls for 80 genotyped bats. The recording protocol was as follows: bats were placed individually in a semi-open environment on a cloth wrapped over the end of a 2-meter pole to allow them to orientate and choose their moment of take-off freely. When the bat left the pole, acoustic recordings were made at a sampling rate of 448 kHz with an Acer Aspire One laptop computer (model KAV60, Acer Inc., Taiwan) using the Avisoft-UltraSoundGate 116H and the Avisoft-RECORDER USHG software (Avisoft Bioacoustics, Germany). Recordings were semi-automatic, with manual activation, a pre-trigger of 2 seconds and a post-trigger of 5 seconds to ensure the acquisition of full call sequences. The condenser microphone CM16/CMPA used (Avisoft Bioacoustics, Germany) had a sensitivity ranging from 10 to 200 kHz. The datasets of wing and echolocation calls overlapped for 35% of the analyzed individuals. The overlap of the datasets in terms of individuals is of minor concern here. We used microsatellite clusters (see further methods) to identify species for the individuals found in the different datasets. Our approach benefited from larger sample sizes that are more representative of the populations studied.

## Molecular analyses and species identification

A subset of the captured individuals (n = 27) was clearly identified as *M. bondae* based on their size and darker pelage coloration [22]. The species status of these 27 individuals was therefore not checked with molecular methods. For the remaining 935 individuals of *M. molossus* and
*M. coibensis*, we used molecular methods; specifically i) for genetic clustering of nuclear microsatellite markers and ii) for phylogenetic tree reconstruction using 659 base pairs of the mitochondrial gene cytochrome oxidase subunit 1 (co1) and 615 base pairs of the hyper variable fragment of the control region (d-loop). Laboratory work with these markers was initially targeted at different questions, i.e. a study of genetic population structure in *M. molossus* as well as an exploration of fur color variation. This explains the use of different markers as well as protocols and number of individuals in each analysis.

**Microsatellite development and genotyping.** The detailed laboratory protocol for the nuclear microsatellite markers is available in the S1 Text. Eighteen primer pairs successfully amplified; we report the sequences, accession numbers for the NCBI Probe database, the fluorescent dyes and the multiplex combinations in S1 Table. We used these 18 microsatellite markers to genotype 935 individuals.

**Microsatellite evaluation and clustering.** To identify the number of species captured, we performed microsatellite-based clustering of 935 genotyped individuals. This aim was achieved in three steps: i) genetic clustering of individuals based on the 18 microsatellite loci, ii) assertion of different assumptions for the genetic analysis software packages used in this study (Hardy-Weinberg Equilibrium, low frequency of null alleles and linkage equilibrium) and iii) genetic clustering based on a robust, filtered set of those loci that adhered closely to the respective genetic assumptions. We first determined the number of genetic clusters corresponding to the number of species (at least two). We used the 18 microsatellite loci using a two-step Discriminant Analysis of Principal Components (DAPC [28]), a clustering method that does not require specific genetic assumptions for the loci used (unlike other clustering software that typically make use of patterns, e.g., Hardy Weinberg and linkage equilibria [28]). The second step consisted of checking three genetic assumptions within each cluster defined by DAPC: Hardy-Weinberg Equilibrium (HWE), low frequency of null alleles and linkage equilibrium. Only loci following these three conditions in each cluster were used for a second stringent clustering analysis performed using the software STRUCTURE 2.3.4 [29,30].

For the first part of the microsatellite analysis, we selected the number of genetic clusters (corresponding to the different species) based on Bayesian Information Criterion (BIC), a measure of the trade-off between goodness of fit and complexity of the model. We calculated the BIC for 18 clusters (the number of buildings sampled) and 100 PCs with the *adegenet* package [31] in R 3.1.0 (R Development Core Team 2014). A two-step Discriminant Analysis of Principal Components (DAPC) [28] was used to infer the selected number of clusters. We retained the number of principal component axes corresponding to ~80% of the cumulative score in the Principal Component Analysis step and the number of axis corresponding to the optimized *a*-score in the Discriminant Analysis step [32].

For each cluster defined with the DAPC, we identified the number of alleles at each locus, the heterozygosity (observed and expected), tested for deviations of HWE and estimated the null allele frequency using CERVUS v3.0.3 [33]. For each cluster, we also tested for linkage disequilibrium between all pairs of loci using the log likelihood ratio statistic and default parameters implemented in GENEPOP ON THE WEB [34,35] and we applied a Bonferroni correction to the significance level of 0.05 (0.05: 9 loci at HWE = 0.00556) to correct for multiple testing. For the following steps, we selected only loci that were in HWE, had an estimated null allele frequency < 0.10 with CERVUS, and were in linkage equilibrium for all clusters. It has recently been shown that null allele estimation with CERVUS can be misleading [36]. We therefore additionally used the software ML-NULL, which has been shown to perform best among a number of methods [36,37], to obtain additional estimates and confirm frequencies < 0.10. The outcomes of both methods (i.e., CERVUS and ML-NULL) did not differ in our case (results not shown).
The last clustering analyses were based on the selected number of genetic clusters in the data and only those loci closely following the genetic assumptions of HWE, null alleles and linkage equilibrium. As a complementary method to the two-step DAPC (following the procedure described earlier), we ran an analysis with the software STRUCTURE 2.3.4 [29,30]. We used default parameters from the software with an admixture model, a length of burn-in period of 20,000 and a number of MCMC repetitions after burn-in of 80,000. We performed 10 replicate runs for the number of determined genetic clusters and averaged the results in CLUMPP 1.1.2 [38]. A few individuals showed lower membership probability to a genetic cluster with STRUCTURE (≤ 0.9), even though showing a strong assignment with DAPC. We excluded these admixed individuals, potentially hybrids or attributed to the wrong species, to avoid potential mistakes in subsequent analyses. Indeed, it is known that DAPC can be over-confident in making genetic cluster assignments and more than one method should be utilized to check for cluster assignment [39]. The pruned dataset was used to identify the number of alleles for each cluster.

**Sequencing and phylogenetic reconstructions.** We sequenced *co1* for 96 individuals and d-loop for 150 individuals. The detailed laboratory protocol for the mitochondrial genes is available in the S1 Text. The newly generated sequences are available on GenBank, respectively under the accession numbers KT721362—KT721412 for the 51 *co1* sequences and KT721413—KT721441 and KT721443—KT721563 for the 150 d-loop sequences.

We obtained 74 *co1* sequences from GenBank, including all sequences for *Molossus* and four outgroups from the molossid family (three species of *Cynomops* and one species of *Promops*). We aligned the 51 *co1* sequences from this study with the 74 GenBank sequences using MUSCLE [40] with default parameters as implemented in SeaView 4.5.4. [41]. We aligned the d-loop sequences from this study with MEGA 4.0 [42] and visually checked the alignment for repeated sequence arrays, a pattern already found in different bat species [43]. Three *Cynomops* species (GenBank accession numbers JF447634, JN312044 and EF080319) and *Promops centralis* (JF444936) were used as outgroups to root the tree inferred from *co1* sequences. The last sequence was labelled as *M. rufus* but we verified the species using the Barcode of Life Data System (more than 29000 sequences for the Order Chiroptera, [12]). The tree inferred from d-loop sequences was unrooted because we could not find publicly available sequences of close outgroups that could be satisfactorily aligned with our sequences. In order to find the best-fitting model for each gene, we compared 56 models of nucleotide evolution using jModelTest 2.1.7 and the Bayesian Information Criterion (BIC) [44]. The best-fitting model was then used in PAUP* 4b10 [45] to infer the respective phylogenies. Reliability of nodes was measured using 100 non-parametric bootstraps that were then mapped on the inferred trees using the *plotBS* option in the R package *phangorn* [46]. We validated the taxonomic identification of the sequences deposited in GenBank a posteriori (see discussion). The information on genetic clustering from the STRUCTURE analysis was also plotted on the tips of the final trees.

## Variation of fur color

We selected a set of eight individuals from the three species with pictures of fur color from the back of the individuals. This set of individuals was representative of the whole range of fur color observed in the field. This selection of pictures displayed the intra-species variation but also inter-species overlap in fur color. Our further use of the pictures to quantify colors was limited by the absence of camera calibration [47].

## Analyses of body parameters

We investigated morphological species differences based on two parameters: forearm length (mm) and body mass (g). We used these parameters to estimate a linear discriminant function
using the “lda” function (library Mass) in R 3.1.0 [48] to separate the three *Molossus* species.
Only adults, but not pregnant females, were used in the analysis. We calculated means and
95% confidence intervals (CI) for each combination of morphological parameter, species and
sex. We obtain the 95% Confidence Intervals by using the formula provided by the R book
[49]. We also assessed the classification rate of the species by the Lda function with the leave-
one-out cross-validation procedure.

Analyses of wing shape

We used the wing photos to extract a series of wing parameters and morphological traits relevant to flight performance and foraging strategy [50]. We followed an established procedure to define landmarks and obtain the following measurements [51] from wing photos (Fig 1): forearm length (mm), total area (mm²), total wing length (mm), arm wing area (mm²), arm wing length (mm), hand wing area (mm²), hand wing length (mm), wing aspect ratio (wing length²/wing area), wing loading (body mass*g/wing area), tip length ratio (hand wing length / arm wing length), tip area ratio (hand wing area / arm wing area), tip shape index (tip area ratio / tip length ratio—tip area ratio) and a circularity index (4*π*wing area / wing perimeter²). In order to minimize inter-observer error, all measurements were collected by a single individual observer. We calculated the mean and the 95% CI [49] for each combination of wing parameter, species and sex.

<figure><img src="image_5.png" /><figcaption>Fig 1. Right wing of a *Molossus molossus* showing areas used to analyze wing shape.</figcaption></figure>

doi:10.1371/journal.pone.0150780.g001
## Analyses of echolocation calls

We analyzed echolocation calls from a subset of individuals genetically identified as *M. molossus* or *M. coibensis* using the software Batsound 4.1.300 (Pettersson Electronik AB, Uppsala, Sweden). We randomly selected sequences of up to ten calls that contained a sufficient signal to noise ratio for each individual. Sampling frequency was configured at 44.1 kHz, with 16 bits per sample and a 512-point FFT with a Hamming window for analysis. A 112 Hz frequency resolution was obtained for spectrograms and power spectrum. In each call, we measured six echolocation parameters using the software Batsound (Pettersson Elektronik AB, Sweden). From the spectrogram, based on the fundamental call, we measured 1) the Start Frequency (SF; frequency measured at the beginning of the call), 2) the End Frequency (EF; frequency measured at the end of the call) and 3) the bandwidth (BW; difference between SF and EF) in kHz. From the maximal intensity in the power spectrum, we determined the 4) Peak Frequency (PF). From the oscillogram, we extracted 5) the Duration (D) and 6) the Pulse Interval (PI; time interval between two consecutive calls) in ms.

First, we analyzed all calls to examine the entire recorded acoustic diversity. We found a great range of variability in the calls, consistent with previous studies on *M. molossus* in Belize and Cuba [52,53]. We examined the Pearson's product moment correlation using R 3.1.0 [48]. Only two of the acoustic parameters (SF and PF) showed a strong correlation of 0.95 (all others ranging from -0.69 to 0.85). We excluded PF and ran a Principal Component Analysis (PCA) of all calls with the five remaining acoustic parameters. Secondly, we categorized our different sequences of calls into call types. A typical sequence of calls started at the release perch with short calls with a downward frequency modulation and a prominent second harmonic, similar to the approach call recorded for *M. molossus* in the vicinity of their roosts in Cuba [54]. We also recorded search flight calls with narrow bandwidths [54] when a bat was higher above the ground. Search flight calls were typically two-toned and alternating between a lower frequency pulse (SI) and a higher frequency pulse (SII) [53,54]. Some search flight calls were also irregularly alternating the SI and SII or were three-toned, a known pattern for this species [55,56]. For our purpose of species comparisons, we selected only sequences with a clear call structure: the approach calls where all calls had harmonics and the two-toned search flight calls consistently alternated with a lower and higher frequency pulse (SI and SII, Fig 2). For each combination of call type and species, we calculated mean and the 95% CI [49]. We disregarded sequences of calls that could not be firmly categorized such as sequences of approach calls that did not always show harmonics, sequences mixing approach calls and search flight calls as well as search flight calls irregularly alternating the tones or showing an uncertain number of tones.

## Results

### Microsatellite evaluation and clustering

We genotyped 935 individuals at 18 microsatellite loci (for dataset, see S1 Dataset). Based on the complete dataset with 18 loci, we selected *K* = 2 clusters because of the shape of the BIC curve as a function of the number of clusters (ranging from one to 18), showing a much better likelihood for *K* = 2 than for *K* = 1 and only little gain in likelihood for additional clusters (S1 Fig). In the two-step DAPC, we retained 60 axes (~80% of the cumulative variance) in the Principal Component Analysis step and one axis (optimized *a*-score) in the Discriminant Analysis step. From the 935 individuals, 841 were attributed to *Cluster One* and 94 to *Cluster Two*.

For these genetic clusters, we list the number of alleles, the range of allele size, the observed and expected heterozygosities and the estimated null allele frequency in Table 1. Two loci from *Cluster One* and three from *Cluster Two* significantly departed from Hardy-Weinberg
<figure><img src="image_6.png" /><figcaption>**Fig 2.** Sonograms of the two types of calls informative for identification found in *Molossus molossus* and *coibensis*. (A) Approach calls with harmonics. (B) Search calls alternating a lower and higher frequency pulse (respectively SI and SII).</figcaption></figure>

doi:10.1371/journal.pone.0150780.g002
<table><caption>**Table 1.** Cross-amplification and genetic tests for 18 *Molossus molossus* loci grouped in two genetic clusters. The columns respectively represent: A, Number of alleles; AS, range of allele sizes (bp); Ho, observed heterozygosity and He, expected heterozygosity; F(null), estimated null allele frequency. Loci or values highlighted in boldface departed significantly from HWE (following p-values from testing in CERVUS) or had high estimated null allele frequencies (> 0.10 in CERVUS).</caption><thead><tr><th rowspan="2">Genetic cluster<br>Locus</th><th colspan="4">Cluster One (n = 841)</th><th colspan="4">Cluster Two (n = 94)</th></tr><tr><th>A</th><th>AS (bp)</th><th>Ho / He</th><th>F(null)</th><th>A</th><th>AS (bp)</th><th>Ho / He</th><th>F(null)</th></tr></thead><tbody><tr><td>Mol_A2</td><td>9</td><td>191-211</td><td>0.74 / 0.76</td><td>0.02</td><td>7</td><td>189-203</td><td>0.65 / 0.68</td><td>0.01</td></tr><tr><td>Mol_A221</td><td>11</td><td>286-321</td><td>0.45 / 0.47</td><td>0.02</td><td>7</td><td>286-315</td><td>0.14 / 0.18</td><td><strong>0.16</strong></td></tr><tr><td>Mol_B233</td><td>4</td><td>198-209</td><td>0.51 / 0.54</td><td>0.02</td><td>3</td><td>198-205</td><td>0.02 / 0.08</td><td><strong>0.48</strong></td></tr><tr><td>Mol_C3</td><td>11</td><td>257-284</td><td>0.81 / 0.80</td><td>-0.01</td><td>11</td><td>261-286</td><td>0.82 / 0.8</td><td>-0.02</td></tr><tr><td>Mol_C6</td><td>7</td><td>101-117</td><td>0.52 / 0.60</td><td>0.06</td><td>6</td><td>100-107</td><td><strong>0.16 / 0.53</strong></td><td><strong>0.55</strong></td></tr><tr><td>Mol_C20</td><td>4</td><td>136-142</td><td>0.71 / 0.72</td><td>0.00</td><td>5</td><td>143-151</td><td>0.58 / 0.59</td><td>0.00</td></tr><tr><td>Mol_C27</td><td>23</td><td>270-320</td><td>0.83 / 0.87</td><td>0.02</td><td>7</td><td>268-304</td><td>0.33 / 0.40</td><td>0.09</td></tr><tr><td>Mol_C56</td><td>18</td><td>171-213</td><td>0.88 / 0.84</td><td>-0.02</td><td>12</td><td>186-210</td><td>0.83 / 0.84</td><td>0.00</td></tr><tr><td>Mol_C61</td><td>19</td><td>177-214</td><td><strong>0.78 / 0.89</strong></td><td>0.07</td><td>6</td><td>182-198</td><td>0.59 / 0.61</td><td>0.01</td></tr><tr><td>Mol_C77</td><td>11</td><td>198-225</td><td>0.71 / 0.73</td><td>0.02</td><td>5</td><td>198-231</td><td>0.27 / 0.30</td><td>0.05</td></tr><tr><td>Mol_C109</td><td>17</td><td>255-293</td><td>0.60 / 0.87</td><td><strong>0.18</strong></td><td>7</td><td>243-277</td><td>0.31 / 0.79</td><td><strong>0.44</strong></td></tr><tr><td>Mol_C109bis</td><td>18</td><td>218-247</td><td><strong>0.81 / 0.83</strong></td><td>0.01</td><td>6</td><td>214-226</td><td>0.61 / 0.63</td><td>0.02</td></tr><tr><td>Mol_C114</td><td>12</td><td>268-310</td><td>0.78 / 0.76</td><td>-0.02</td><td>10</td><td>268-313</td><td>0.73 / 0.81</td><td>0.05</td></tr><tr><td>Mol_C115</td><td>12</td><td>265-294</td><td>0.75 / 0.78</td><td>0.02</td><td>7</td><td>265-282</td><td>0.61 / 0.67</td><td>0.05</td></tr><tr><td>Mol_C117</td><td>15</td><td>294-348</td><td>0.28 / 0.85</td><td><strong>0.51</strong></td><td>6</td><td>298-340</td><td>0.57 / 0.55</td><td>-0.01</td></tr><tr><td>Mol_C118</td><td>10</td><td>213-224</td><td>0.55 / 0.81</td><td><strong>0.19</strong></td><td>6</td><td>214-224</td><td><strong>0.06 / 0.51</strong></td><td><strong>0.81</strong></td></tr><tr><td>Mol_C132</td><td>15</td><td>147-182</td><td>0.77 / 0.81</td><td>0.03</td><td>4</td><td>178-186</td><td>0.42 / 0.46</td><td>0.03</td></tr><tr><td>Mol_D109</td><td>20</td><td>291-324</td><td>0.86 / 0.90</td><td>0.02</td><td>6</td><td>296-317</td><td>0.03 / 0.07</td><td><strong>0.42</strong></td></tr></tbody></table>

doi:10.1371/journal.pone.0150780.t001

equilibrium (HWE). Three loci from *Cluster One* and six from *Cluster Two* showed high estimated null allele frequencies (over 10%). Two of the loci from *Cluster Two* departing from the HWE also had high estimated null allele frequency, potentially resulting from null amplification. Of the nine loci at HWE, many pairs showed significant linkage disequilibrium (22 for *Cluster One* and nine for *Cluster Two* out of 36). The only loci in HWE, in linkage equilibrium and with estimated null allele frequency < 0.10 across the two clusters were C56, C77, C115 and C132.

We consequently based all following clustering analyses on only four loci and two clusters. We also excluded two individuals with missing data for these specific loci. Some individuals retained also showed missing data at two loci (n = 6) and one locus (n = 87) of these four. The performance of the DAPC and STRUCTURE analyses on four loci matched that of the analyses with 18 loci resulting in two similar genetic clusters. The majority of individuals were clearly found in *Cluster One* and *Cluster Two* (respectively orange and blue on Fig 3). Only ten individuals out of 933 (1%) showed discrepancies in clustering, with clear assignment to a cluster in the DAPC but admixture in the STRUCTURE analysis (posterior assignment probability < 0.9). These individuals with uncertain assignment were removed from the subsequent analyses as explained in the methods. The pruned dataset was composed of 923 individuals: 833 in *Cluster One* that occurred in all 18 sampled buildings and 90 in *Cluster Two* found in five of the 18 sampled buildings. These two clusters were used to determine the number of alleles for each of the 18 loci (S2 Table).

## Phylogenetic results: the co1 tree

The final co1 alignment consisted of 659 nucleotides for 125 individuals (S2 Dataset). Phylogenetic reconstruction under a Tamura-Nei (TrN) model with rate variation among sites (Γ) and
<figure><img src="image_7.png" /><figcaption>**Fig 3.** Genetic clustering (K = 2) of 933 *Molossus* bats from the village of Gamboa. *Cluster One* is represented in orange and *Cluster Two* in blue. The two clusters were obtained with four microsatellite primers using the software STRUCTURE. The few individuals at the edge of the two clusters are admixed.</figcaption></figure>

doi:10.1371/journal.pone.0150780.g003

a proportion of invariable sites (I) displayed numerous polytomies with most nodes showing low statistical support (<70%) (S2 Fig). The majority of samples genotyped with co1 (96.1%) were also genotyped using microsatellites, the membership to the microsatellite clusters are represented in color at the nodes of the tree. In the *Molossus* crown, bootstrap values were >70% in only five nodes. The tree consisted mostly of i) deep rooted clades comprising newly sampled *M. molossus* from Panama together with published haplotypes from Ecuador, Suriname and Guyana, ii) a “floating” clade from Panama and iii) a tree crown with two *M. rufus*, four *M. coibensis* and a few more individuals from Panama. The tree crown was composed of a well-supported clade with two *M. rufus* (BS = 100) and a polytomy composed of a *M. sp.* from Venezuela (JF447833), a clade with four *M. coibensis* and a *Molossus* sp. (JF442201) from Ecuador and a clade with 11 of our individuals from Panama. As these 11 bats were also found in *Cluster Two* from the microsatellite clustering analysis, we later considered all these animals to be *M. coibensis*. The deeper branches of the tree consisted of *M. molossus* from Panama, Ecuador, Suriname and Guyana. At the roots of the tree, the 22 individuals from Panama were mixed together with the GenBank sequences of *M. molossus*, with no clear biogeographic pattern. We also found a well-supported clade with 17 bats from Panama (BS = 98), sister to the crown tree. However, the 22 individuals mixed with the GenBank sequences and the 17 individuals from this Panamanian clade were previously grouped in the same *Cluster One* from the microsatellite clustering analysis. We therefore did not consider the Panamanian clade as a new species, but rather a “floating” clade and defined all individuals from *Cluster One* as *M. molossus* in this tree and subsequent analyses.

## Phylogenetic results: the d-loop tree

The final alignment of d-loop was 615 base pair long and encompassed 150 individuals (S3 Dataset). The 150 new d-loop sequences showed high genetic diversity, with 42 different haplotypes. Most haplotypes (n = 28) were found in one individual, except for a common haplotype that was shared by 38 individuals. Fourteen haplotypes were shared between different roosts, one of them being shared between ten roosts. According to the BIC criterion, the Hasegawa-Kishino-Yano model (HKY) with a proportion of invariable sites (I) was the best fitting model for the tree reconstruction (S3 Fig). The d-loop tree presented a similar topology to the co1 tree, with several clades of *M. molossus* and a clade with *M. coibensis*. One of the individuals from this tree (KT721428) was previously identified as *M. coibensis* in the co1 tree (KT721364, transponder number EAC87), we therefore assigned its clade in the d-loop tree to the species *M. coibensis* and the rest of the individuals to *M. molossus*. Statistical support for the *M. coibensis'* clade and six subclades of *M. molossus* was high (BS ≥ 90). Most of the d-loop sequences
<figure><img src="image_8.png" /><figcaption>Fig 4. Variation of fur color (back) for eight individuals of *M. molossus*, *M. coibensis* and *M. bondae* from Panama.</figcaption></figure>

doi:10.1371/journal.pone.0150780.g004

(77.3%) were also genotyped for microsatellites, the membership to the microsatellite clusters are represented in color at the nodes of the tree.

## Variation of fur color

The subset was composed of two *M. bondae* identified morphologically as well as three *M. molossus* and three *M. coibensis* confirmed genetically. We observed a fur color on the back ranging from light brown to dark brown. The similarity of the fur color emerges clearly from this panel of three species (Fig 4).

## Analyses of body parameters

We obtained morphological data from 617 adults of both sexes: 526 *M. molossus*, 64 *M. coibensis* and 27 *M. bondae* (the first two genetically assigned to species; S4 Dataset). Forearm length (mm) and body mass (g) were normally distributed and overlapped between the three species (Fig 5). Forearm length was ranked in increasing order for *M. coibensis*, *M. molossus* and *M. bondae*. Body mass was ranked in increasing order for *M. molossus*, *M. coibensis* and *M. bondae*. Only *M. molossus* showed body mass below 9.5 g while only *M. bondae* had forearm length above 39.63 mm and body mass above 17 g. Twenty-five of the 617 bats (4.1%) were misclassified by the `lda` function using forearm length and body mass. Misclassification occurred for animals with extreme values for the species range. Thus, the lightest *M. bondae* (n = 3, range = 10.5–12.0 g) and the heaviest *M. molossus* (n = 3, range = 16.0–17.0 g) were wrongly identified as well the smallest *M. molossus* (n = 9, range = 31.6–36.6 mm) and the largest *M. coibensis* (n = 10, range = 34.5–37.02 mm).

We provide means, 95% CI and range of forearm length and body mass in Table 2 for each combination of species and sex. We found a significant male-biased sexual dimorphism in all three species at the intra-specific level, as the 95% CI did not overlap. Inter-specific values differed significantly for both parameters based on the 95% CI.

Based on the leave-one-out cross-validation procedure, the overall classification rate of the function was 95.9% (97.7% for *M. molossus*, 84.4% for *M. coibensis* and 88.9% for *M. bondae*). Only 25 of 617 individuals (4.1%) were misclassified based on the combined two morphological parameters alone (symbols outlined in black in Fig 4).

## Analyses of wing shape

We analyzed wing photos (S5 Dataset) of 104 *M. molossus* (87 females and 17 males) and 12 *M. coibensis* (8 females and 4 females). The means and 95% CI for each combination of wing parameter, species and sex are summarized in Table 3. Four parameters were significantly
<figure><img src="image_9.png" /><figcaption>**Fig 5.** Forearm length (mm) plotted against body mass (g) for three Panamanian *Molossus* species. The color code is as follows: *M. molossus* (orange dots), *M. coibensis* (blue diamonds) and *M. bondae* (green squares). Following the same color code, the frequency distribution of body mass is plotted above the graph and the frequency distribution of forearm length on the right side of the graph. Points outlined in black are misclassified individuals based on the linear discriminant function and the leave-one-out cross-validation procedure (4.1% of the individuals).</figcaption></figure>

doi:10.1371/journal.pone.0150780.g005

different between species: *Molossus molossus* had longer forearms (confirming the results outlined in the previous paragraph), larger hand wing area, and a longer hand wing, while wing loading was greater in *M. coibensis*. Wrongly classified by the Ida function as *M. molossus*, a female *M. coibensis* (E8472) had values for the wing parameters falling within the 95% CI of *M. molossus*.

<table><caption>**Table 2.** Sex-specific means, [95% CI], sample size and range of forearm length and body mass.</caption><thead><tr><th rowspan="2">Parameter</th><th colspan="2">Forearm length</th><th colspan="2">Body mass (g)</th></tr><tr><th>Female</th><th>Male</th><th>Female</th><th>Male</th></tr></thead><tbody><tr><td rowspan="3">M. coibensis</td><td>34.6</td><td>35.6</td><td>12.8</td><td>14.1</td></tr><tr><td>[34.3–34.8]</td><td>[35.3–35.9]</td><td>[12.3–13.3]</td><td>[13.5–14.8]</td></tr><tr><td>(n = 41, 33.0–37.0)</td><td>(n = 23, 34.0–37.0)</td><td>(n = 41, 9.5–16.2)</td><td>(n = 23, 11.7–17.0)</td></tr><tr><td rowspan="3">M. molossus</td><td>36.7</td><td>37.4</td><td>10.2</td><td>11.4</td></tr><tr><td>[36.6–36.8]</td><td>[37.2–37.6]</td><td>[10.1–10.3]</td><td>[11.1–11.6]</td></tr><tr><td>(n = 403, 31.9–38.8)</td><td>(n = 123, 31.4–39.6)</td><td>(n = 403, 7.5–16.0)</td><td>(n = 123, 7.7–17.0)</td></tr><tr><td rowspan="3">M. bondae</td><td>39.4</td><td>40.6</td><td>16.4</td><td>17.9</td></tr><tr><td>[39.0–39.8]</td><td>[39.8–41.3]</td><td>[14.9–18.0]</td><td>[15.8–20.0]</td></tr><tr><td>(n = 18, 38.3–41.4)</td><td>(n = 9, 39.0–42.2)</td><td>(n = 18, 10.5–25.0)</td><td>(n = 9, 13.5–21.5)</td></tr></tbody></table>

doi:10.1371/journal.pone.0150780.t002
## Analyses of echolocation calls

We recorded echolocation calls of 80 adult bats: 65 *M. molossus* and 15 *M. coibensis*. We measured 8 to 30 calls per individual, resulting in 834 for *M. molossus* and 255 pulses for *M. coibensis* (S6 Dataset). When analyzing unclassified calls, we found no clear species differences based on the Principal Component Analysis (Fig 6), with just 51.9% of the variance explained by the first axis and 21.6% by the second axis.

We found sequences of approach calls with harmonics in 35.4% of the *M. molossus* and 46.7% of the *M. coibensis*. We observed search calls that regularly alternated between the two tones SI and SII in 18.5% of the *M. molossus* and 40.0% of the *M. coibensis*. Average SF and EF were higher in approach calls of *M. coibensis* but lower in the two-toned calls. For the two-toned calls, call duration was shorter in *M. molossus*, and bandwidth was higher in *M. coibensis*. Finally, pulse interval in the SII of the search calls was higher in *M. coibensis*. Mean values and SD for the five acoustic parameters and the different call types are compiled in Table 4 (for results of t-tests, see S3 Table). None of the individuals misclassified by the Lda function had categorized calls to be compared to the table.

<table><caption>**Table 3.** Mean [95% CI] of wing parameters for species and sex. The four parameters highlighted in bold differed significantly between species. For each species, the intermediate column compares mean values between sexes.</caption><thead><tr><th rowspan="2">Species</th><th colspan="2">Molossus coibensis</th><th colspan="2">Molossus molossus</th></tr><tr><th>Female (n = 8)</th><th>Male (n = 4)</th><th>Female (n = 87)</th><th>Male (n = 17)</th></tr></thead><tbody><tr><td>Forearm length (mm)</td><td>33.8<br>[32.5-35.1]</td><td>35.3<br>[35.0-35.6]</td><td>36.7<br>[36.6-36.9]</td><td>37.1<br>[36.7-37.5]</td></tr><tr><td>Total area (mm²)</td><td>2292<br>[2064-2520]</td><td>2469<br>[2233-2706.0]</td><td>2792.1<br>[2739-2845]</td><td>2660.6<br>[2531-2789]</td></tr><tr><td>Total wing length (mm)</td><td>94.8<br>[89.3-100.2]</td><td>99.7<br>[92.5-107.0]</td><td>103.8<br>[102.7-104.9]</td><td>102.7<br>[100.1-105.3]</td></tr><tr><td>Arm wing area (mm²)</td><td>1086<br>[939-1233]</td><td>1263<br>[1103-1422]</td><td>1354<br>[1319-1388]</td><td>1321<br>[1226-1416]</td></tr><tr><td>Arm wing length (mm)</td><td>35.1<br>[32.0-38.2]</td><td>39.3<br>[34.4-44.2]</td><td>38.8<br>[38.3-39.4]</td><td>38.5<br>[36.9-40.1]</td></tr><tr><td>Hand wing area (mm²)</td><td>1179<br>[1092-1266]</td><td>1211<br>[1107-1315]</td><td>1432<br>[1404-1459]</td><td>1361<br>[1318-1403]</td></tr><tr><td>Hand wing length (mm)</td><td>59.9<br>[57.2-62.5]</td><td>60.0<br>[57.6-62.3]</td><td>64.9<br>[64.0-65.7]</td><td>64.1<br>[62.6-65.6]</td></tr><tr><td>Circularity</td><td>0.49<br>[0.48-0.51]</td><td>0.48<br>[0.47-0.49]</td><td>0.51<br>[0.50-0.51]</td><td>0.49<br>[0.48-0.50]</td></tr><tr><td>Tip area ratio</td><td>1.1<br>[1.0-1.2]</td><td>1.0<br>[0.8-1.1]</td><td>1.1<br>[1.0-1.1]</td><td>1.0<br>[1.0-1.1]</td></tr><tr><td>Tip length ratio</td><td>1.7<br>[1.6-1.8]</td><td>1.5<br>[1.4-1.7]</td><td>1.7<br>[1.6-1.7]</td><td>1.7<br>[1.6-1.7]</td></tr><tr><td>Tip shape index</td><td>1.8<br>[1.5-2.1]</td><td>1.7<br>[1.1-2.3]</td><td>1.9<br>[1.7-2.0]</td><td>1.7<br>[1.5-1.9]</td></tr><tr><td>Aspect ratio</td><td>3.9<br>[3.8-4.1]</td><td>4.0<br>[3.8-4.3]</td><td>3.9<br>[3.8-3.9]</td><td>4.0<br>[3.9-4.1]</td></tr><tr><td>Wing loading (Nm²)</td><td>53.9<br>[46.9-60.9]</td><td>56.1<br>[47.6-64.5]</td><td>39.3<br>[37.9-40.7]</td><td>43.2<br>[40.6-45.9]</td></tr></tbody></table>

doi:10.1371/journal.pone.0150780.t003
<figure><img src="image_10.png" /><figcaption>Fig 6. Principal Component Analysis of five acoustic parameters for *M. molossus* (orange) and *M. coibensis* (blue).</figcaption></figure>

doi:10.1371/journal.pone.0150780.g006

## Discussion

Our combination of morphometric and molecular data confirmed the sympatry of at least three species: *Molossus molossus*, *M. coibensis* and *M. bondae* at our study site in Panama. *Molossus molossus* was more abundant than *M. coibensis* in the sampled buildings while *M. bondae* was rarely captured.

## Microsatellite clustering

Microsatellite analyses were useful for detecting species but also to reveal potential interspecific hybrids. Both the DAPC BIC method based on 18 loci and the STRUCTURE method based on
<table><caption>**Table 4.** Comparison of acoustic parameters between *M. molossus* and *M. coibensis*. Values are means ± 1 standard deviation. Values in boldface represent significant differences between species based on a Student's t-test for the given acoustic parameter. The two figures for sample size indicate the number of individuals and the number of calls.</caption><thead><tr><th>Call type</th><th>Species</th><th>Sample size</th><th>Start Fr. (kHz)</th><th>End. Fr.(kHz)</th><th>Bandwidth (kHz)</th><th>Duration (ms)</th><th>Pulse interval (ms)</th></tr></thead><tbody><tr><td rowspan="2">Approach call</td><td><i>M. coibensis</i></td><td>7 (70)</td><td><b>56.1 ± 2.9</b></td><td><b>30.6 ± 5.5</b></td><td>25.5 ± 6.4</td><td>0.4 ± 0.05</td><td>62.3 ± 54.9</td></tr><tr><td><i>M. molossus</i></td><td>23 (226)</td><td><b>52.8 ± 3.7</b></td><td><b>29.0 ± 3.2</b></td><td>23.8 ± 4.9</td><td>0.4 ± 0.08</td><td>50.7 ± 37.3</td></tr><tr><td rowspan="2">Two-toned SI</td><td><i>M. coibensis</i></td><td>6 (29)</td><td><b>35.4 ± 1.3</b></td><td><b>29.8 ± 1.9</b></td><td><b>5.6 ± 1.3</b></td><td><b>0.4 ± 0.4</b></td><td>76.2 ± 23.6</td></tr><tr><td><i>M. molossus</i></td><td>12 (60)</td><td><b>39.1 ± 3.6</b></td><td><b>34.4 ± 3.8</b></td><td><b>4.7 ± 2.0</b></td><td><b>0.6 ± 0.3</b></td><td>75.4 ± 24.6</td></tr><tr><td rowspan="2">Two-toned SII</td><td><i>M. coibensis</i></td><td>6 (29)</td><td><b>39.7 ± 2.1</b></td><td><b>35.0 ± 1.8</b></td><td><b>4.7 ± 1.2</b></td><td><b>0.3 ± 0.3</b></td><td><b>153.8 ± 58.9</b></td></tr><tr><td><i>M. molossus</i></td><td>12 (60)</td><td><b>42.8 ± 2.9</b></td><td><b>39.1 ± 3.0</b></td><td><b>3.8 ± 1.4</b></td><td><b>0.6 ± 0.3</b></td><td><b>117.6 ± 44.5</b></td></tr></tbody></table>

doi:10.1371/journal.pone.0150780.t004

four loci, recovered two genetic clusters, with consistent cluster membership with few exceptions. Our clustering analyses are therefore very robust to the choice of markers. Only 1% of the individuals clearly assigned to one cluster with the DAPC method showed admixture based on STRUCTURE. Sequencing a subset of these individuals with the gene *co1* revealed that the two microsatellite clusters corresponded to *M. molossus* and *M. coibensis*. As we successfully genotyped both species using microsatellites, they offer the potential for cross-species amplification in the genus *Molossus*. Multiple and non-exclusive explanations such as non-random mating could be an explanation as to why so many of the loci diverged from HWE and/or showed high prevalence of null alleles [57]. We may have sampled a non-random subset of the males in the gene pool [58]; in particular, our sampling was biased towards harem social groups occupying buildings [59] whereas unsampled males were probably solitary and roosting elsewhere. The four loci that stringently followed the genetic assumptions for the STRUCTURE analyses revealed a low number of admixed individuals (n = 10), representing only 1% of all genotyped individuals. Admixture signatures in STRUCTURE can result from interspecific hybridization or retention of ancestral polymorphism [15,60–62]. Admixture can also result from microsatellite size homoplasy, a well-described phenomenon that remains infrequently tested [63,64]. Further investigation of admixture in this study, be it technical artifact or biological reality, is not of relevance here because of its low incidence, and therefore out of scope presently. Further analysis will be required to identify the mechanism leading to admixture and the taxonomic status of admixed individuals.

## Phylogenetic reconstructions

Sympatry of very similar-looking species is common in bats [3,65,66] and can make species identification in the field extremely difficult. We successfully clustered individuals according to species with the *co1* sequences. In addition, the *co1* tree allowed us to incorporate GenBank sequences from different species and countries. We thus matched our sequences from Panama to *M. molossus* from Guyana, Ecuador and Suriname, and to *M. coibensis* from Ecuador. Our phylogenetic reconstructions provided the second piece of molecular evidence that *M. coibensis* and *M. molossus* occur in sympatry in Panama following the allozyme study by Dolan [21]. We found that the two species occur in the same buildings where they probably form species-specific social groups. Sympatric individuals of these two species have previously been reported for the province of Napo in Ecuador [9,67]. Our phylogenetic tree also revealed a “floating” clade of *M. molossus* that we were not able to match to GenBank sequences. With low statistical support (BS = 23), this phylogenetic uncertainty may represent a soft polytomy that could be resolved with increased character sampling [68]. The same 17 individuals were assigned to *M. molossus* in the microsatellite clustering analyses. This molecular differentiation may result
from the occurrence of two distinct barcodes in the same species, as recently found in the bat *Pipistrellus kuhlii* [14]. The phylogenetic tree also revealed three GenBank sequences that were probably wrongly identified: JF442201 from Ecuador and JF447833 from Venezuela are probably not *M. molossus* but *M. coibensis*, and the inverse is true for JF442240 from Ecuador. Quality control of sequences using phylogenetic analyses [69] could easily avoid taxonomic misidentification of sequences submitted in public databases [70]. Despite its utility for species identification, our phylogenetic reconstructions using *co1* recovered a large polytomy with limited statistical support for the majority of nodes, and further phylogenetic reconstruction based on the fast-evolving mitochondrial region d-loop also recovered clades with low support in most cases. Future studies incorporating nuclear genes [68,71,72], combined datasets (i.e. morphology and genetics) [73] or even complete genomes [74–76] will thus be valuable in further elucidating the taxonomy of morphologically similar molossid bats. Whole genome analyses of thousands of species have been envisioned for many years [77] and this is becoming a reality with the development of next generation sequencing technologies, for example, birds with the B10K initiative [78,79].

## Variation of fur color

The panel of fur color (Fig 4) reveals the overlap between species, namely between *M. molossus* and *M. coibensis*.

## Analyses of body parameters

In contrast to the molecular methods, simple morphological parameters such as forearm length and body mass can easily be obtained in the field [80]. Used separately, these morphological parameters did not allow good discrimination of the three species due to the overlap in parameter distributions and the flip in ranks. Only when analyzed together in a linear discriminant function, the two parameters led to a high average rate of correct identification of the three species (95.9%). At the species level, the classification rate was ranked in decreasing order for *M. molossus* (97.7%), *M. bondae* (88.9%) and *M. coibensis* (84.4%). Similarly, for example, *Myotis* from Switzerland can be most reliably distinguished using a canonical discriminant function based on the forearm length and the ear length [81]. *Rhinolophus* from Bulgaria, Greece and Turkey can be correctly assigned using a canonical discriminant function of the forearm length and the first phalanx of the fourth finger [82] but for *Rhinolophus* from Tunisia the second phalanx of the fourth finger has to be added to the discriminant function [83]. The discriminant analysis constitutes a powerful approach to differentiate between morphologically similar species but only if the right combination of parameters from correctly assigned subsets of the individuals can be found. In addition to the species differences in body mass and forearm length, we also observed significant sexual size dimorphism with larger and heavier males in the three species. The sample size is low for *M. coibensis* (eight females and four males) but the reference values of the two parameters for each sex (Table 3) should allow correct species identification for most individuals of our three focal species in Panama.

## Analyses of wing shape

Wing shape is under strong selection for ecological niche use in bats because it underpins flight performance and foraging strategy [50]. Wing shape can also be useful for species identification [51,82]. However, all molossid bats have narrow wings well-adapted to hunting insects in open spaces [84] and, therefore, many of the wing parameters that we measured did not vary between species. We found higher values of forearm length (mm), hand wing area (mm²) and hand wing length (mm) in *M. molossus*, indicating longer wings at our study site. However, *M*.
*coibensis* had higher wing loading (Nm⁻²), suggesting higher flight speed and turning agility than *M. molossus* [50]. While such variation may potentially be ecologically significant for niche separation between the two species, significant intraspecific geographic variation in wing parameters can be found, for example in rhinolophids [82]. Until a dataset from a broader geographic range is available, our values should only be used for comparisons with individuals within Central America or even only Panama.

## Analyses of echolocation calls

Acoustic recording of free-flying bats is a widespread method for surveys and species differentiation including different species of *Molossus* [52,55,85]. Acoustic recordings after release such as ours are also commonly used but not, to date, in molossids. The method is generally criticized as being not representative of the environment and associated calls in free-flying animals [86] but remains invaluable to match acoustic and molecular data of the same individual.
Acoustic recordings after release proved useful only for a subset of our recordings after we categorized into different types of calls. The manual categorization of the calls confirmed the previously described call diversity: approach calls with harmonics described in *M. molossus* [53,54], two-toned search calls (described in the Molossidae and Vespertilionidae [87,88]) and three-toned search calls (*M. molossus* [55]). We only selected the two unmistakable categories: short calls with a downward frequency modulation and a prominent second harmonic (approach calls) and the alternating two-toned calls (search calls). Only a few call parameters showed significant differences between the two species, especially start and end frequency. Despite the apparent utility of these calls to discriminate species, several limitations to this method must be considered. For example, mean values and standard deviations for SI and SII strongly overlapped between *M. molossus* and *M. coibensis* (e.g. 0.3 ± 0.3 ms & 0.4 ± 0.4 ms vs. 0.6 ± 0.3 ms and 0.6 ± 0.3 ms). Following the recommendation of Barclay [89], our reference values should only be compared to individuals released under the open sky and away from background clutter.

## Comparison to other studies

Previous studies have reported reference values for different morphological parameters. However, several of these studies provided reference values using a different taxonomy from the one we used. For example, Simmons [19] considered *M. bondae* and *M. currentium* grouped under the name *M. currentium* and Reid [25] treated *M. coibensis* as a subspecies of *M. molossus*. Studies that followed the same taxonomy as we did provided values consistent with our results (based on molecular validation): for example, a range of 33.2–36.0 mm, 33.5–34.7 mm or 33.9–36.0 mm for *M. coibensis* [22,90,91], 35–40 mm for *M. molossus* [22] and 38.4–41.1 mm or 38–43 mm for *M. bondae* [22,25]. Measurements of body mass showed similar values too, with a range of 16–21 g for *M. bondae* [25]. However, larger values of forearm length in *M. coibensis* and *M. molossus* were found in a study from southeastern Brazil [92]. *M. coibensis* showed a range of 35.5–38.1 mm for females and 36.5–38.1 mm for males while *M. molossus* showed a range of 35.3–42.2 mm for females and 38.2–42.3 mm for males. The inconsistency between the values reported in our study in Panama and the one in Brazil could be a result of natural geographic variation—a common phenomenon in bats [93]—or from incorrect taxonomic attribution leading to wrong values in other studies. Problems in taxonomic attribution can result from unsettled taxonomy, for example *M. coibensis* from French Guyana referred to as a true species [22] or as the separate species *M. barnesi* [19,91]. To tackle these issues, additional comparative studies using molecular validation are required to provide reference values for these species in other regions of their biogeographic ranges.
## Conclusions

All methods we used, based on molecular, morphometric or acoustic data, provided useful information for species discrimination. However, all of these methods had their limitations too. While we were able to reliably identify *M. bondae* based on size and pelage coloration, the microsatellite analysis led to a reliable genetic clustering of *M. coibensis* and *M. molossus* using two different methods. Individuals were assigned correctly with just a few exceptions when using all 18 microsatellite loci as well as with the more stringently determined subset of markers (n = 4). However, developing microsatellite primers involves a considerable effort. Other molecular methods, PCR-based assays [94] or high resolution melting [95] are promising alternatives that should allow cheaper and faster results of similar quality. The phylogenetic reconstructions with the *co1* sequences were also useful for species identification but mitochondrial DNA markers alone did not provide strong clade support. Comparison of morphometric parameters, i.e. forearm length and body mass, was a simple and very useful tool for species discrimination. However, they were only discriminatory when combined in a linear discriminant analysis function or when the sex of the individuals was taken into account. Previous work on other species shows that a different trait combination may need to be found for each local species assembly, which may only allow correct species assignment after fieldwork has been completed, similar to the molecular methods. This may be particularly true for species-rich regions such as our study area where cryptic species are still being described and material for identification is patchy or even lacking. Only four of the 13 wing parameters we included in our analysis differed between species: forearm length, hand wing area, hand wing length and wing loading. Training the dataset on a subset of individuals was necessary to obtain reliable rules of thumb that can be used in the field and then again, potentially only locally. Finally, only a subset of the commonly used call recordings revealed species-specific differences in different acoustic parameters. This may be due to the artificial situation of a release, however as recording of free-flying bats cannot be matched to DNA or morphological measurements, it remains difficult to assess whether analysis of these calls would be more reliable even if species could be identified in such a situation.

Although any single morphological measurement proved to be unreliable for species identification, we found that by combining multiple measurements we could reliably identify the focal *Molossus spp.* in Panama, as verified by the genetic data. Based on these findings we emphasize the importance of combining morphological traits for field identification, as well as using independent genetic data to help decide upon the best combination of these traits in any given location. Proper species identification is the important basis for any work with wild animals and thus distinguishing a focal species within a local species complex may only be possible using a multi-method approach.

## Supporting Information

**S1 Dataset. Microsatellite dataset with 935 *Molossus spp.* from Panama genotyped at 18 loci.**

(CSV)

**S2 Dataset. Alignment of 659 nucleotides of the mitochondrial gene *co1* for 125 *Molossus spp.*** Consists of 51 individuals from Panama newly sequenced in this study and 74 molossid sequences extracted from GenBank.

(FAS)
**S3 Dataset. Alignment of the mitochondrial region d-loop (615 bp) for 150 *Molossus* spp. from Panama.** These were newly sequenced for this study.
(FAS)

**S4 Dataset. Dataset of morphological parameters for adults of three species of *Molossus* from Panama (.csv).** The columns respectively represent the transponder number of each individual (1), date of capture (2), roost building (3), species name (4), species' attribution based on the linear discriminant analysis (5), sex (6), mass in grams (7) and forearm length in mm (8).
(CSV)

**S5 Dataset. Dataset of wing parameters for *M. coibensis* and *M. molossus* from Gamboa, Panama (.csv).** The columns respectively represent the transponder number of each individual (1), date of capture (2), species name (3), sex (4), mass of the bat in g (5), forearm length in mm (6), total area in mm² (7), total wing length in mm (8), arm wing area in mm² (9), arm wing length in mm (10), hand wing area in mm² (11), hand wing length in mm (12), circularity (13), tip area ratio (14), tip length ratio (15), tip shape index (16), aspect ratio (17) and wing loading in Nm⁻² (18).
(CSV)

**S6 Dataset. Dataset of echolocation call parameters for *M. coibensis* and *M. molossus* from Gamboa, Panama (.csv).** The columns respectively represent the transponder number of each individual (1), species name (2), type of echolocation call (approach calls or search calls (3), subtype of call, either SI or SII for "search call" (4), start frequency of the pulse in kHz (5), end frequency of the pulse in kHz (6), bandwidth of the pulse in kHz (7), duration of the pulse in ms (8) and pulse interval in ms (9).
(CSV)

**S1 Fig. Values of Bayesian Information Criterion as a function of the number of clusters (one to 18).** This figure was used to determine the number of clusters to retain in the microsatellite clustering analysis. The dataset contained 935 *Molossus spp.* genotyped at 18 loci.
(PDF)

**S2 Fig. Maximum likelihood tree of the genus *Molossus* obtained from *co1* (659 bp).** The tree was created using the TrN+I+Γ substitution model and PAUP*. Three outgroups from the genus *Cynomops* were removed for visual display of the tree. Bootstrap percentages from ML analyses above 50, obtained from maximum likelihood analyses (see methods for the tree reconstruction), are indicated at the nodes. The orange and blue colors at tip labels correspond with the two genetic clusters identified with the STRUCTURE analysis and white tips indicate GenBank sequences.
(PDF)

**S3 Fig. Maximum likelihood tree of the genus *Molossus* obtained from d-loop (615 bp).** The tree was created using the HKY+I substitution model and PAUP*. Three outgroups from other bat genera were removed to display the tree. Bootstrap percentages from ML analyses above 50, obtained from maximum likelihood analyses (see methods for the tree reconstruction), are indicated at the nodes. The orange and blue colors at tip labels correspond with the two genetic clusters identified with the STRUCTURE analysis and white tips indicate sequences without species attribution.
(PDF)
**S1 Table.** Primer pairs for 18 microsatellite loci derived from *Molossus molossus*.
(XLSX)

**S2 Table.** Table summarizing allele amplification at 18 loci for 923 *Molossus* spp. divided into two genetic clusters. The four loci in bold represent the four loci retained for the final clustering analysis.
(XLSX)

**S3 Table.** Table with results from t-test comparisons of acoustic parameters between *M. coibensis* and *M. molossus*. A combination of three types of calls (rows) and five acoustic parameters (columns) were investigated. For each comparison, the t-test, the degree of freedom (df) and the p-value were provided.
(XLSX)

**S1 Text.** Laboratory protocol for microsatellites and mitochondrial sequences.
(DOCX)

## Acknowledgments

We thank Rachel Page and the members of the batlab as well as Julia Cramer, Marion Muturi, Teague O'Mara, Sebastian Rikker and Sebastian Stockmaier for helping in the field. We are grateful to James Borell, Evi Fricke, Mimoza Hoti and Monika Struebig for helping in the lab, Manuel Ruedi and M. Teague O'Mara for advising on the analyses, Karla Bauer for collecting the wing data as well as Andreas Kiefer, Kim G. Mortega, Conor Whelan and two anonymous reviewers for commenting on earlier versions of the manuscript. We also thank the World Bat Library (Museum of Geneva, Switzerland) for providing access to rare publications.

## Author Contributions

Conceived and designed the experiments: YG ET RHSK AL DKND. Performed the experiments: YG ET DL MM DKND. Analyzed the data: YG FBC MM. Contributed reagents/materials/analysis tools: AL SJR. Wrote the paper: YG ET DL MM FBC SJR RHSK AL DKND.

## References

1. 1331 and counting. In: BATS magazine [Internet]. 2015. Available: http://www.batcon.org/resources/media-education/bats-magazine/bat_article/1506

2. Mayer F, Dietz C, Kiefer A. Molecular species identification boosts bat diversity. Front Zool. 2007; 4: 4. doi: 10.1186/1742-9994-4-4 PMID: 17295921

3. Ibáñez C, García-Mudarra JL, Ruedi M, Stadelmann B, Juste J. The Iberian contribution to cryptic diversity in European bats. Acta Chiropterologica. 2006; 8: 277-297. doi: 10.3161/1733-5329(2006)8[277: TICTCD]2.0.CO;2

4. Puechmaille SJ, Allegrini B, Boston ESM, Dubourg-Savage M-J, Evin A, Knochel A, et al. Genetic analyses reveal further cryptic lineages within the *Myotis nattereri* species complex. Mamm Biol—Zeitschrift für Säugetierkd. 2012; 77: 224–228. doi: 10.1016/j.mambio.2011.11.004

5. Kiefer A, Mayer F, Kosuch J, Von Helversen O, Veith M. Conflicting molecular phylogenies of European long-eared bats (*Plecotus*) can be explained by cryptic diversity. Mol Phylogenet Evol. 2002; 25: 557–566. PMID: 12450759

6. Bogdanowicz W, Hulva P, Černá Bolfíková B, Buš MM, Rychlicka E, Sztencel-Jablonka A, et al. Cryptic diversity of Italian bats and the role of the Apennine refugium in the phylogeography of the western Palaearctic. Zool J Linn Soc. 2015; 174: 635–648. doi: 10.1111/zoj.12248

7. Hebert PDN, Gregory TR. The promise of DNA barcoding for taxonomy. Syst Biol. 2005; 54: 852–859. doi: 10.1080/10635150500354886 PMID: 16243770

8. Francis CM, Borisenko A V, Ivanova N V, Eger JL, Lim BK, Guillén-Servent A, et al. The role of DNA barcodes in understanding and conservation of mammal diversity in southeast Asia. PLoS One. 2010; 5: e12575. doi: 10.1371/journal.pone.0012575 PMID: 20838635
9. Clare EL, Lim BK, Fenton MB, Hebert PDN. Neotropical bats: estimating species diversity with DNA barcodes. PLoS One. 2011; 6: e22648. doi: 10.1371/journal.pone.0022648 PMID: 21818359

10. Wilson J-J, Sing K-W, Halim MRA, Ramli R, Hashim R, Sofian-Azirun M. Utility of DNA barcoding for rapid and accurate assessment of bat diversity in Malaysia in the absence of formally described species. Genet Mol Res. 2014; 13: 920–5. doi: 10.4238/2014.February.19.2 PMID: 24634112

11. Floyd R, Abebe E, Papert A, Blaxter M. Molecular barcodes for soil nematode identification. Mol Ecol. 2002; 11: 839–850. PMID: 11972769

12. Ratnasingham S, Hebert PDN. bold: The Barcode of Life Data System (http://www.barcodinglife.org). Mol Ecol Notes. 2007; 7: 355–364. doi: 10.1111/j.1471-8286.2007.01678.x PMID: 18784790

13. Rubinoff D, Cameron S, Will K. A genomic perspective on the shortcomings of mitochondrial DNA for "barcoding" identification. J Hered. 2006; 97: 581–594. doi: 10.1093/jhered/esl036 PMID: 17135463

14. Andriollo T, Naciri Y, Ruedi M. Two mitochondrial barcodes for one biological species: the case of European Kuhl's pipistrelles (Chiroptera). PLoS One. 2015; 10: e0134881. doi: 10.1371/journal.pone.0134881 PMID: 26241944

15. Berthier P, Excoffier L, Ruedi M. Recurrent replacement of mtDNA and cryptic hybridization between two sibling bat species *Myotis myotis* and *Myotis blythii*. Proc R Soc B. 2006; 3101–3109. doi: 10.1098/rspb.2006.3680 PMID: 17018432

16. Galimberti A, Spada M, Russo D, Mucedda M, Agnelli P, Crottini A, et al. Integrated operational taxonomic units (IOTUs) in echolocating bats: a bridge between molecular and traditional taxonomy. PLoS One. 2012; 7: e40122. doi: 10.1371/journal.pone.0040122 PMID: 22761951

17. Padial JM, Miralles A, De la Riva I, Vences M. The integrative future of taxonomy. Front Zool. 2010; 7: 16. doi: 10.1186/1742-9994-7-16 PMID: 20500846

18. Miller GS. Notes on bats of the genus *Molossus*. Proc United States Natl Museum. 1913; 46: 85–92.

19. Simmons NB. Order Chiroptera. In: Wilson DE, Reeder DM, editors. Mammal species of the world: a taxonomic and geographic reference. Third edit. Johns Hopkins University Press; 2005. pp. 312–529.

20. Nowak RN. Walker's bats of the world. John Hopkins University Press; 1994.

21. Dolan PG. Systematics of Middle American mastiff bats of the genus *Molossus*. Lubbock: Texas Tech University Press; 1989.

22. Eger JL. Family Molossidae. In: Gardner AL, editor. Mammals of South America Marsupials, Xenarthrans, Shrews, and Bats, vol 1. Chicago and London: The University of Chicago Press; 2008. pp. 399–440.

23. González-Ruiz N, Ramírez-Pulido J, Arroyo-Cabrales J. A new species of mastiff bat (Chiroptera: Molossidae: Molossus) from Mexico. Mamm Biol—Zeitschrift für Säugetierkd. 2011; 76: 461–469. doi: 10.1016/j.mambio.2010.06.004

24. Koopman KF. Chiroptera: Systematics. In: Niethammer J, Schliemann H, Starck D, editors. Handbuch der Zoologie VIII. Berlin: Walter de Gruyter; 1994. pp. 1–217.

25. Reid FA. A field guide to the mammals of Central America and Southeast Mexico. Oxford: Oxford University Press; 1998.

26. Ammerman LK, Lee DN, Tipps TM. First molecular phylogenetic insights into the evolution of free-tailed bats in the subfamily Molossinae (Molossidae, Chiroptera). J Mammal. 2012; 93: 12–28. doi: 10.1644/11-MAMM-A-103.1

27. Worthington Wilmer J, Barrett E. A non-lethal method of tissue sampling for genetic studies of chiropterans. Bat Res News. 1996; 37: 1–3.

28. Jombart T, Devillard S, Balloux F. Discriminant analysis of principal components: a new method for the analysis of genetically structured populations. BMC Genet. 2010; 11: 94. doi: 10.1186/1471-2156-11-94 PMID: 20950446

29. Pritchard JK, Stephens M, Donnelly P. Inference of population structure using multilocus genotype data. Genetics. 2000; 155: 945–959. PMID: 10835412

30. Falush D, Stephens M, Pritchard JK. Inference of population structure using multilocus genotype data: linked loci and correlated allele frequencies. Genetics. 2003; 164: 1567–87. PMID: 12930761

31. Jombart T. adegenet: a R package for the multivariate analysis of genetic markers. Bioinformatics. 2008; 24: 1403–5. doi: 10.1093/bioinformatics/btn129 PMID: 18397895

32. Jombart T, Collins C. A tutorial for Discriminant Analysis of Principal Components (DAPC) using adegenet 2.0.0. 2015.

33. Kalinowski ST, Taper ML, Marshall TC. Revising how the computer program CERVUS accommodates genotyping error increases success in paternity assignment. Mol Ecol. 2007; 16: 1099–106. doi: 10.1111/j.1365-294X.2007.03089.x PMID: 17305863
34. Raymond M, Rousset F. GENEPOP (Version 1.2): Population genetics software for exact tests and ecumenicism. J Hered. 1995; 86: 248–249.

35. Rousset F. *genepop*007: a complete re-implementation of the *genepop* software for Windows and Linux. Mol Ecol Resour. 2008; 8: 103–6. doi: 10.1111/j.1471-8286.2007.01931.x PMID: 21585727

36. Dąbrowski MJ, Bornelőv S, Kruczyk M, Baltzer N, Komorowski J. "True" null allele detection in microsatellite loci: a comparison of methods, assessment of difficulties and survey of possible improvements. Mol Ecol Resour. 2015; 15: 477–88. doi: 10.1111/1755-0998.12326 PMID: 25187238

37. Kalinowski ST, Taper ML. Maximum likelihood estimation of the frequency of null alleles at microsatellite loci. Conserv Genet. 2006; 7: 991–995. doi: 10.1007/s10592-006-9134-9

38. Jakobsson M, Rosenberg NA. CLUMPP: a cluster matching and permutation program for dealing with label switching and multimodality in analysis of population structure. Bioinformatics. 2007; 23: 1801–6. doi: 10.1093/bioinformatics/btm233 PMID: 17485429

39. Frosch C, Kraus RHS, Angst C, Allgöwer R, Michaux J, Teubner J, et al. The genetic legacy of multiple beaver reintroductions in Central Europe. PLoS One. 2014; 9: e97619. doi: 10.1371/journal.pone.0097619 PMID: 24827835

40. Edgar RC. MUSCLE: multiple sequence alignment with high accuracy and high throughput. Nucleic Acids Res. 2004; 32: 1792–7. doi: 10.1093/nar/gkh340 PMID: 15034147

41. Gouy M, Guindon S, Gascuel O. SeaView version 4: A multiplatform graphical user interface for sequence alignment and phylogenetic tree building. Mol Biol Evol. 2010; 27: 221–4. doi: 10.1093/molbev/msp259 PMID: 19854763

42. Tamura K, Dudley J, Nei M, Kumar S. MEGA4: Molecular Evolutionary Genetics Analysis (MEGA) software version 4.0. Mol Biol Evol. 2007; 24: 1596–9. doi: 10.1093/molbev/msm092 PMID: 17488738

43. Wilkinson GS, Mayer F, Kerth G, Petri B. Evolution of repeated sequence arrays in the D-loop region of bat mitochondrial DNA. Genetics. 1997; 146: 1035–48. PMID: 9215906

44. Darriba D, Taboada GL, Doallo R, Posada D. jModelTest 2: more models, new heuristics and parallel computing. Nature Methods. 2012. p. 772. doi: 10.1038/nmeth.2109

45. Swofford DL. PAUP*. Phylogenetic Analysis Using Parsimony (*and Other Methods). Sunderland, Massachusetts: Sinauer Associates; 2003.

46. Schliep KP. phangorn: phylogenetic analysis in R. Bioinformatics. 2011; 27: 592–3. doi: 10.1093/bioinformatics/btq706 PMID: 21169378

47. Stevens M, Párraga CA, Cuthill IC, Partridge JC, Troscianko TS. Using digital photography to study animal coloration. Biol J Linn Soc. 2007; 90: 211–237. doi: 10.1111/j.1095-8312.2007.00725.x

48. R Core Team. R: A language and environment for statistical computing. Vienna, Austria: R Foundation for Statistical Computing; 2015.

49. Crawley MJ. The R Book. John Wiley & Sons Ltd; 2007.

50. Norberg UM, Rayner JM V. Ecological morphology and flight in bats (Mammalia; Chiroptera): wing adaptations, flight performance, foraging strategy and echolocation. Philos Trans R Soc Lond B Biol Sci. 1987; 316: 335–427. doi: 10.1098/rstb.1987.0030

51. Schmieder DA, Benítez HA, Borissov IM, Fruciano C. Bat species comparisons based on external morphology: a test of traditional versus geometric morphometric approaches. PLoS One. 2015; 10: e0127043. doi: 10.1371/journal.pone.0127043 PMID: 25965335

52. O'Farrell MJ, Miller BW. Use of vocal signatures for the inventory of free-flying Neotropical bats. Biotropica. 1999; 31: 507–516. doi: 10.1111/j.1744-7429.1999.tb00394.x

53. Kossl M, Mora E, Coro F, Vater M. Two-toned echolocation calls from *Molossus molossus* in Cuba. J Mammal. 1999; 80: 929–932. doi: 10.2307/1383262

54. Mora EC, Macías S, Vater M, Coro F, Kössl M. Specializations for aerial hawking in the echolocation system of *Molossus molossus* (Molossidae, Chiroptera). J Comp Physiol A. 2004; 190: 561–74. doi: 10.1007/s00359-004-0519-2

55. Jung K, Molinari J, Kalko EK V. Driving factors for the evolution of species-specific echolocation call design in new world free-tailed bats (Molossidae). PLoS One. 2014; 9: e85279. doi: 10.1371/journal.pone.0085279 PMID: 24454833

56. Barataud M, Giosa S, Leblanc F, Rufray V, Disca T, Tillon L, et al. Identification et écologie acoustique des chiroptères de Guyane française. Le Rhinolophe. 2013; 19: 103–145.

57. Hartl DL, Clark AG. Principles of population genetics. Cornell University; 2007.

58. Storz JF. Genetic consequences of mammalian social structure. J Mammal. 1999; 80: 553–569. doi: 10.2307/1383301
59. McCracken GF, Wilkinson GS. Bat mating systems. In: Crichton EG, Krutzsch PH, editors. Reproductive biology of bats. San Diego: Academic Press; 2000. pp. 321–362.

60. Muir G, Schlötterer C. Evidence for shared ancestral polymorphism rather than recurrent gene flow at microsatellite loci differentiating two hybridizing oaks (*Quercus* spp.). Mol Ecol. 2005; 14: 549–61. doi: 10.1111/j.1365-294X.2004.02418.x PMID: 15660945

61. Randi E. Detecting hybridization between wild species and their domesticated relatives. Mol Ecol. 2008; 17: 285–93. doi: 10.1111/j.1365-294X.2007.03417.x PMID: 18173502

62. Brown RM, Nichols RA, Faulkes CG, Jones CG, Bugoni L, Tatayah V, et al. Range expansion and hybridization in Round Island petrels (*Pterodroma* spp.): evidence from microsatellite genotypes. Mol Ecol. 2010; 19: 3157–3170. doi: 10.1111/j.1365-294X.2010.04719.x PMID: 20618891

63. Adams RI, Brown KM, Hamilton MB. The impact of microsatellite electromorph size homoplasy on multilocus population structure estimates in a tropical tree (*Corythophora alta*) and an anadromous fish (*Morone saxatilis*). Mol Ecol. 2004; 13: 2579–88. doi: 10.1111/j.1365-294X.2004.02256.x PMID: 15315672

64. Rossiter SJ, Benda P, Dietz C, Zhang S, Jones G. Rangewide phylogeography in the greater horseshoe bat inferred from microsatellites: implications for population history, taxonomy and conservation. Mol Ecol. 2007; 16: 4699–4714. doi: 10.1111/j.1365-294X.2007.03546.x PMID: 17927706

65. Barratt EM, Deaville R, Burland TM, Bruford MW, Jones G, Racey PA, et al. DNA answers the call of pipistrelle bat species. Nature. 1997; 387: 138–9. doi: 10.1038/387138b0

66. von Helversen O, Heller KG, Mayer F, Nemeth A, Volleth M, Gombkötö P. Cryptic mammalian species: a new species of whiskered bat (*Myotis alcathoe* n. sp.) in Europe. Naturwissenschaften. 2001; 88: 217–23. PMID: 11482435

67. McDonough MM, Ferguson AW, Ammerman LK, Granja-Vizcaino C, Burneo SF, Baker RJ. Molecular verification of bat species collected in Ecuador: results of a country-wide survey. Occas Pap Museum Texas Tech Univ. 2011; 301: 1–28.

68. Lack JB, Van Den Bussche RA. Identifying the confounding factors in resolving phylogenetic relationships in Vespertilionidae. J Mammal. 2010; 91: 1435–1448. doi: 10.1644/09-MAMM-A-354.1

69. Botero-Castro F, Delsuc F, Douzery EJP. Thrice better than once: quality control guidelines to validate new mitogenomes. Mitochondrial DNA. 2015; 1–6. doi: 10.3109/19401736.2014.900666

70. Nilsson RH, Ryberg M, Kristiansson E, Abarenkov K, Larsson K-H, Köljalg U. Taxonomic reliability of DNA sequences in public sequence databases: a fungal perspective. PLoS One. 2006; 1: e59. doi: 10.1371/journal.pone.0000059 PMID: 17183689

71. Ruedi M, Stadelmann B, Gager Y, Douzery EJP, Francis CM, Lin L-K, et al. Molecular phylogenetic reconstructions identify East Asia as the cradle for the evolution of the cosmopolitan genus *Myotis* (Mammalia, Chiroptera). Mol Phylogenet Evol. 2013; 69: 437–449. doi: 10.1016/j.ympev.2013.08.011 PMID: 23988307

72. Velazco PM, Patterson BD. Diversification of the yellow-shouldered bats, genus *Stumira* (Chiroptera, Phyllostomidae), in the New World tropics. Mol Phylogenet Evol. 2013; 68: 683–98. doi: 10.1016/j.ympev.2013.04.016 PMID: 23632030

73. Springer MS, Teeling EC, Madsen O, Stanhope MJ, De Jong WW. Integrated fossil and molecular data reconstruct bat echolocation. Proc Natl Acad Sci U S A. 2001; 98: 6241–6246. PMID: 11353869

74. Delsuc F, Brinkmann H, Philippe H. Phylogenomics and the reconstruction of the tree of life. Nat Rev Genet. 2005; 6: 361–75. doi: 10.1038/nrg1603 PMID: 15861208

75. Murphy WJ, Pevzner PA, O'Brien SJ. Mammalian phylogenomics comes of age. Trends Genet. 2004; 20: 631–9. doi: 10.1016/j.tig.2004.09.005 PMID: 15522459

76. Song S, Liu L, Edwards S V, Wu S. Resolving conflict in eutherian mammal phylogeny using phylogenomics and the multispecies coalescent model. Proc Natl Acad Sci U S A. 2012; 109: 14942–7. doi: 10.1073/pnas.1211733109 PMID: 22930817

77. Genome 10K Community of Scientists. Genome 10K: a proposal to obtain whole-genome sequence for 10,000 vertebrate species. J Hered. 2009; 100: 659–74. doi: 10.1093/jhered/esp086 PMID: 19892720

78. Zhang G. Genomics: Bird sequencing project takes off. Nature. 2015; 522: 34–34. doi: 10.1038/522034d

79. Kraus RHS, Wink M. Avian genomics: fledging into the wild! J Ornithol. 2015; doi: 10.1007/s10336-015-1253-y

80. Dietz C, von Helversen O. Illustrated identification key to the bats of Europe. 2004.

81. Arlettaz R, Ruedi M, Hausser J. Field morphological identification of *Myotis myotis* and *Myotis blythii* (Chiroptera, Vespertilionidae): A multivariate approach. Myotis. 1991; 29: 7–16.
82. Dietz C, Dietz I, Siemers BM. Wing measurement variations in the five european horseshoe bat species (Chiroptera: Rhinolophidae). J Mammal. 2006; 87: 1241–1251. doi: 10.1644/05-MAMM-A-299R2.1

83. Puechmaille SJ, Hizem WM, Allegrini B, Abiadh A. Bat fauna of Tunisia: Review of records and new records, morphometrics and echolocation data. Vespertilio. 2012; 16: 211–239.

84. Voigt CC, Holderied MW. High manoeuvring costs force narrow-winged molossid bats to forage in open space. J Comp Physiol B. 2012; 182: 415–24. doi: 10.1007/s00360-011-0627-6 PMID: 22048527

85. Jung K, Kalko EK V. Adaptability and vulnerability of high flying Neotropical aerial insectivorous bats to urbanization. Divers Distrib. 2011; 17: 262–274. doi: 10.1111/j.1472-4642.2010.00738.x

86. O'Farrell MJ, Miller BW, Gannon WL. Qualitative identification of free-flying bats using the Anabat detector. J Mammal. 1999; 80: 11–23. doi: 10.2307/1383203

87. Fenton MB, Portfors C V, Rautenbach IL, Waterman JM. Compromises: sound frequencies used in echolocation by aerial-feeding bats. Can J Zool. 1998; 76: 1174–1182.

88. Kingston T, Jones G, Zubaid A, Kunz T. Alternation of echolocation calls in five species of aerial feeding insectivorous bats. J Mammal. 2003; 84: 205–215.

89. Barclay RMR. Bats are not birds—a cautionary note on using echolocation calls to identify bats: a comment. J Mammal. 1999; 80: 290–296. doi: 10.2307/1383229

90. Correa da Costa LJ, Gonçalves de Andrade FA, Uieda W, Gregorin R, Barroncas Fernandes ME. First record of *Molossus coibensis* (Chiroptera: Molossidae) in the Brazilian Amazon. Mastozool Neotrop. 2013; 20: 143–147.

91. Gregorin R, Tahara AS, Buzzato DF. *Molossus aztecus* and other small *Molossus* (Chiroptera: Molossidae) in Brazil. Acta Chiropterologica. 2011; 13: 311–317. doi: 10.3161/150811011X624794

92. Pimenta V, Fonseca BS, Hoppe JPM, Ditchfield AD. First occurrence of *Molossus coibensis* Allen, 1904 (Chiroptera, Molossidae) in Atlantic Forest. Chiropt Neotrop. 2014; 20: 1237–1242.

93. Storz JF, Balasing J, Bhat HR, Nathan PT, Doss DPS, Prakash AA, et al. Clinal variation in body size and sexual dimorphism in an Indian fruit bat, *Cynopterus sphinx* (Chiroptera: Pteropodidae). Biol J Linn Soc. 2001; 72: 17–31. doi: 10.1111/j.1095-8312.2001.tb01298.x

94. Boston ESM, Hanrahan N, Puechmaille SJ, Ruedi M, Buckley DJ, Lundy MG, et al. A rapid PCR-based assay for identification of cryptic *Myotis* spp. (*M. mystacinus*, *M. brandtii* and *M. alcathoe*). Conserv Genet Resour. 2011; 3: 557–563. doi: 10.1007/s12686-011-9404-9

95. Wittwer CT. High-resolution DNA melting analysis: advancements and limitations. Hum Mutat. 2009; 30: 857–9. doi: 10.1002/humu.20951 PMID: 19479960