The taxonomic status of Myotisnesopoluslarensis (Chiroptera, Vespertilionidae) and new insights on the diversity of Caribbean Myotis

Abstract Myotisnesopolus currently comprises two subspecies. The nominate subspecies (M.n.nesopolus) occurs on the Caribbean islands of Curaçao and Bonaire, Netherlands Antilles, whereas M.n.larensis is known from mainland South America in northeastern Colombia and northwestern Venezuela. Our Maximum Likelihood phylogenetic analyses of cytochrome-b gene sequences recovered M.nesopolus as a paraphyletic group, with M.n.nesopolus and M.n.larensis as non-sister lineages. The haplotype network indicates that these two subspecies do not share any haplotypes and are in different evolutionary trajectories. Additionally, these two subspecies can be distinguished on the basis of qualitative and quantitative morphological traits. This pattern supports the recognition of M.nesopolus and M.larensis as full species. Our results also reveal that the assemblage of Caribbean Myotis do not form a monophyletic group. Caribbean species are phylogenetically close to mainland species from northern South America and Central America, suggesting that colonization of Caribbean islands happened multiple times.


Specimens examined
Specimens of M. nesopolus used in this study are deposited in the American Museum of Natural History (AMNH, New York, USA), Carnegie Museum of Natural History (CM, Pittsburgh, USA), Smithsonian's National Museum of Natural History (USNM, Washington DC, USA), and Museum of Texas Tech University (TTU, Lubbock, USA). We examined the holotype of M. n. nesopolus (USNM 101849), two topotypes from Curaçao (CM 52432, USNM 105128), and nine specimens from Bonaire (Appendix 1). Material of M. n. larensis includes the holotype (AMNH 130709), and fifteen additional specimens from mainland Venezuela.

Molecular analyses
Phylogenetic analyses of complete cytochrome-b gene (cyt-b, 1,140 bp, no gaps) sequences were conducted for the Neotropical assemblage of Myotis. A total of 122 sequences, including outgroups, were retrieved from GenBank (Appendix 2). We used the palearctic species Myotis brandtii (Eversmann, 1845) and Myotis gracilis Ognev, 1927 as outgroups because they are sister to the Neotropical clade (see Ruedi et al. 2013). Multiple sequence alignment of full length cyt-b sequences were performed with MEGA X (Kumar et al. 2018), using MUSCLE algorithm with default settings (Edgar 2004). Subsequently, the Bayesian Information Criterion (BIC), as implemented in JModelTest2 (Darriba et al. 2012), was used to determine the best-fit models of nucleotide substitution. The Hasegawa-Kishino-Yano model (Hasegawa et al. 1985) was chosen to correct the heterogeneity rate using gamma-distribution with invariant sites (i.e., HKY + Γ + I).
The phylogenetic analysis was carried out using Maximum Likelihood (ML) method (Felsenstein 1981), in the software RAxML v8.0 (Stamatakis 2014). To assess the nodal support, we calculated a nonparametric bootstrap using 1000 replications. Genetic distance values for cyt-b sequences were calculated in MEGA X using the Kimura 2-parameter model (Kimura 1980).
To understand the population structure of M. n. nesopolus, M. n. larensis and other phylogenetically related population groups, we built a haplotype network (distribution of haplotypes by previously defined population groups) using the median-joining algorithm in the Network 4.6.1.3 software (Bandelt et al. 1999).

Morphological and morphometric analyses
We examined 284 specimens for the morphological comparisons, including M. n. nesopolus (N = 10), M. n. larensis (N = 9) and 14 species of Neotropical Myotis deposited in 11 collections in Brazil, Canada and United States (Appendix 1). Specimens were identified following Wilson (2008) and Moratelli et al. ( , 2013Moratelli et al. ( , 2017. The main qualitative morphological characters used in the comparisons were: (i) presence and height of sagittal crest; (ii) presence and height of lambdoidal crests; (iii) inclination shape of the frontal and parietal bones; (iv) presence of a fringe of hairs along the trailing edge of the uropatagium; (v) dorsal and ventral fur texture and height; (vi) pattern of fur coloring, with the capitalized color nomenclature following Ridgway (1912).
We took one external and 16 craniodental measurements (Table 1), using digital calipers to the nearest 0.01 mm. Measurements were made under binocular microscopes with low magnification (usually 6×). Measurements were recorded from adults and are reported in millimeters (mm). The length of ear and body mass were recorded from skin labels. We used a principal component analysis (PCA) to identify general trends of cranial size and shape variation among samples, and a discriminant function analysis (DFA), with a priori identification of samples, to compare skull size and shape of M. n. nesopolus (N = 9) and M. n. larensis (N = 9). For these analyses, we selected a subset of 11 craniodental dimensions representing different axes of the length and width of skull, rostrum, and mandible, as follows: greatest length of skull, including incisors (GLS), condylo-incisive length (CIL), mastoid breadth (MAB), braincase breadth (BCB), interorbital breadth (IOB), postorbital breadth (POB), breadth across canines (BAC), breadth across molars (BAM), maxillary toothrow length (MTL), molariform toothrow length (M1-M3), and mandibular toothrow length (MAN). PCA and DFA analyses were run in R software (R Development Core Team 2012) using the MASS and Lattice packages (Venables and Ripley 2002;Sarkar 2008). Because multivariate procedures require complete data sets, missing values (ca 1.5% of the total dataset) were estimated from the existing raw data using the Amelia II package (Honaker Table 1. Description of cranial, mandibular, and external dimensions (and their abbreviations). Lengths were measured from the anteriormost point or surface of the 1 st structure to the posteriormost point or surface of the 2 nd structure, except as specified. Subsequently, an analysis of variance using Mann-Whitney statistics was employed to test whether the population samples differ in cranial dimensions. The comparison was made using p-values and when less than 0.001 were considered as statistically significant. This analysis was run in the software PAST 3.3 (Hammer et al. 2001).

Molecular analyses
The ML phylogeny based on cyt-b sequences indicates that M. nesopolus, as currently recognized, is paraphyletic, with M. n. nesopolus more closely related to an eastern Peruvian unidentified lineage, whereas M. n. larensis was recovered more closely related to an unidentified lineage from western Ecuador ( Fig. 1), although this phylogeny and branching events has low nodal support. These unidentified species from Peru and Ecuador were originally designated as Myotis nigricans by the original collector due to morphological similarities. However, M. nigricans has been recovered as polyphyletic and considered a cryptic species complex in many studies (Moratelli et al. , 2013(Moratelli et al. , 2016(Moratelli et al. , 2017Larsen et al. 2012a). Therefore, we decided not to give a name to the lineages related to M. nesopolus and M. larensis. We emphasize that the previous identification of these specimens as M. nigricans by one of our authors (RJL) in a previous study (Larsen et al. 2012a) indicates that these populations are morphologically distinct from those considered here as M. nesopolus and M. larensis. The Caribbean Myotis species do not form a monophyletic group, being related to Myotis atacamensis (Lataste, 1892) and other mainland putative species. Nevertheless, the phylogenetic relationship of Caribbean Myotis clade is not fully resolved, since a polytomy was recovered among M. The average cyt-b pairwise distance between M. n. larensis and Myotis sp. 1 from western Ecuador is 2.1% ± 0.3; between M. n. nesopolus and Myotis sp. 2 from eastern Peru is 3.8% ± 0.4; and between M. n. nesopolus and M. n. larensis is 4.0% ± 0.3 ( Table 2). Levels of intraspecific variation were less than 0.8% for all recognized and putative species (Table 2).
The haplotype network indicates that there are no haplotypes shared between M. n. nesopolus, M. n. larensis, and phylogenetically close species (Fig. 2). The haplotypes were grouped into small clusters well-distributed among the populations, with no central haplotype. The network indicates spatial structuring with isolation among the population groups tested, agreeing with what was obtained by phylogenetic inference.

Morphological analyses
The first principal component (PC1) accounted for 87% of the total craniometric variation, and represents overall skull size (Fig. 3A, B). Along this axis, scores of M. n. larensis and M. n. nesopolus do not overlap. On the other hand, the two samples overlap broadly along the second principal component (PC2 = 5%) which represents overall skull shape. The distribution of M. n. larensis and M. n. nesopolus samples across size and shape axes in the discriminant analysis (Fig. 3C, D) is similar to that observed in Table 2. Average Kimura 2-parameter genetic distances within (along diagonal) and among (below diagonal) Myotis taxa based on cytochrome-b gene sequences. Boldface value indicates the distance between M. larensis and M. nesopolus. Hyphen indicates groups with a single sequence.  the PCA. Measurements associated with skull and mandible length (GLS, CIL, MAN) and skull width (IOB) were the most useful to discriminate samples (Table 3). Considering that skull axes are represented by the set of measurements used in the morphometric multivariate analysis, these results reveal that M. n. larensis and M. n. nesopolus have distinct skull size and shape. Populations from the Antilles and mainland South America do not overlap in measurements of several characters, which may be useful in distinguishing species: M. n. larensis forearm length ranges from 31.2 to 33.2 mm, and GLS from 13.6 to 14.5 mm; M. n. nesopolus forearm length ranges from 28.2 to 31.0 mm, and GLS from 12.9 to 13.4 mm. The Mann-Whitney test found significant differences in 11 of the 14 measurements tested (Table 4).

Description and comparisons
Myotis larensis is a small-sized bat (total length 78-82 mm; forearm length 31.2-33.2; body mass 3-5 g), morphologically similar to several Neotropical congeners. Ears are moderate in size (length 10-13 mm), and when laid forward extend halfway from eye to nostril. Antitragal notch is barely evident. Membranes are Mummy-brown. Fur on dorsal surface of uropatagium extends slightly past the knees. Plagiopatagium is attached to the foot at toes level by a broad band of membrane. Third metacarpal, tibia, and skull are long in relation to forearm (mean ratios 0. 96, 0.48, and 0.43, respectively;see LaVal (1973)).
Myotis larensis can be distinguished from all Caribbean and South American congeners by qualitative and quantitative traits. It differs from M. nesopolus by its larger size (no overlapping in forearm length and greatest length of skull), presence of sagittal crest, and dorsal fur longer and strongly bicolored. Considering the Myotis species that occurs in the northern South America, M. larensis differs from M. albescens (É. Geoffroy, 1806) by the absence of a fringe of hairs along the trailing edge of the uropatagium; from M. keaysi J. A. Allen, 1914, M. pilosatibialis LaVal, 1973, M. riparius Handley, 1960, and M. simus Thomas, 1901 Genoways and Williams (1979) determined that mainland and island specimens of M. larensis and M. nesopolus, respectively, were morphometrically similar, with Venezuelan specimens slightly smaller than those from Curaçao. As a result, they recognized M. larensis as a subspecies of M. nesopolus, which was followed by subsequent authors (e.g., Simmons 2005;Wilson 2008;Moratelli et al. 2019b). However, our results do not support this arrangement, indicating a morphometric discontinuity and qualitative morphological differences between M. larensis and M. nesopolus.

Discussion
Previous phylogenetic studies based on mitochondrial and nuclear DNA recovered M. nesopolus and M. larensis as sister lineages and questioned the subspecific status of M. larensis because the cyt-b genetic distance of 4% between mainland and Antilles populations suggests a potential for separation at the species level (see Bradley and Baker 2001;Larsen et al. 2012b). However, this study did not include the mainland samples from Ecuador and Peru. Our phylogenetic analyses revealed that M. nesopolus and M. larensis are not sister lineages and do not share haplotypes. The genetic distances between M. nesopolus, M. larensis and their sister species are greater than 2%. About this, Bradley and Baker (2001) indicate that genetic distance values between 2 and 11% from cyt-b sequences had a high probability of being indicative of conspecific populations or valid species and merit additional study concerning specific status. Our investigation found a conspicuous phenotypic discontinuity in variation of both the size and shape of the skull and other external characters. Thus, the strong congruence between the morphological, morphometric and molecular evidence presented here supports the hypothesis that M. larensis represents a full species.
Nevertheless, it is important to mention the limitation of cyt-b gene for establishing species boundaries in the Caribbean clade, particularly between M. larensis and M. sp. 1 from Ecuador and between M. nesopolus and M. sp. 2 from Peru. Although widely used (e.g., Larsen et al. 2012a, b;Moratelli et al. 2016Moratelli et al. , 2017Carrión-Bonilla and Cook 2020), the application of cyt-b data to species delimitation and inference of phylogenetic relationships in Myotis from the Caribbean clade was insufficient. This demonstrates the need to expand the use of new genetic markers for future systematic studies with the Caribbean Myotis assemblage.
With the recognition of M. larensis at the species level hierarchy, M. nesopolus is restricted to Bonaire and Curaçao and is the only species of the genus found in these islands (Fig. 6) Barbados and Grenada (LaVal 1973;Larsen et al. 2012a;Moratelli et al. 2017). However, the taxonomic status of some populations of these species needs to be reassessed. For example, Myotis nyctor was described from Barbados and subsequently recorded from Grenada (LaVal 1973;LaVal and Schwartz 1974;Moratelli et al. 2017). Although our phylogenetic analysis grouped the samples of M. nyctor from Barbados (N = 5) and Grenada (N = 1) in the same clade (Fig. 1), and with low genetic distance between them (ca 0.2%; Table 2), there are qualitative and quantitative morphological differences between specimens from these two islands (see Larsen et al. 2012a). The similarity in the cyt-b sequences between Grenada and Barbados specimens may be explained by the retained ancestral polymorphism due to the very recent separation (Stadelmann et al. 2007;Larsen et al. 2012a).
The biogeographic interpretations made by Larsen et al. (2012b) suggest at least two independent Myotis invasions into the Lesser Antilles, and reverse colonization by Caribbean Myotis to mainland Central and South America-the latter being a welldocumented pattern in other Caribbean bat lineages (Dávalos 2005(Dávalos , 2006(Dávalos , 2010Genoways et al. 2005;Pavan and Marroig 2017;Tavares et al. 2018). In addition, some biogeographic and ecological aspects suggest the need for taxonomic revision of some species. The distance and geographic isolation between Barbados and Grenada (ca 255 km) are greater than between Dominica and Martinique (ca 42 km), each one having a unique Myotis species. Moreover, Barbados and Grenada are separated by the Tobago Basin, with an ocean depth of approximately 2500 m and no ridges that may have connected these two populations during glaciation periods (Speed 1981;Humphrey 1997;Graham 2003). Considering the apparent low vagility and the small home range of Myotis in general (e.g., LaVal and Fitch 1977;Castella et al. 2001;Moratelli et al. 2019b), it is possible that the populations of M. nyctor from these two islands are isolated and on different evolutionary trajectories. The same rationale might be valid for M. dominicensis, where the populations from Guadeloupe and Martinique are isolated by approximately 42 km of sea. However, there are several oceanic ridges between these two islands, which may have served as bridges connecting these two populations during the last glaciation (Speed 1981;Humphrey 1997;Graham 2003). Thus, we suggest that future studies on systematics and biogeography of Caribbean Myotis should focus on the definition of the taxonomic status of island populations from Grenada and Guadeloupe.