The Taxonomic Status of Mazama bricenii and the Significance of the Táchira Depression for Mammalian Endemism in the Cordillera de Mérida, Venezuela

Abstract

We studied the taxonomy and biogeography of Mazama bricenii, a brocket deer classified as Vulnerable by the IUCN, drawing on qualitative and quantitative morphology and sequences of the mitochondrial cytochrome-b gene. We used Ecological Niche Modeling (ENM) to evaluate the hypothesis that M. bricenii of the Venezuelan Cordillera de Mérida (CM) might have become isolated from populations of its putative sister species, Mazama rufina, in the Colombian Cordillera Oriental (CO). This hypothesis assumes that warm, dry climatic conditions in the Táchira Depression were unsuitable for the species. Our analyses did not reveal morphological differences between specimens geographically attributable to M. bricenii and M. rufina, and phylogenetic analyses of molecular data recovered M. bricenii nested within the diversity of M. rufina. These results indicate that M. bricenii should be regarded as a junior synonym of M. rufina. ENM analyses revealed the existence of suitable climatic conditions for M. rufina in the Táchira Depression during the last glacial maximum and even at present, suggesting that gene flow between populations in the CO and CM may have occurred until at least the beginning of the current interglacial period and may continue today. Because this pattern might characterize other mammals currently considered endemic to the CM, we examined which of these species match two criteria that we propose herein to estimate if they can be regarded as endemic to the CM with confidence: (1) that morphological or molecular evidence exists indicating that the putative endemic taxon is distinctive from congeneric populations in the CO; and (2) that the putative endemic taxon is restricted to either cloud forest or páramo, or both. Only Aepeomys reigi, Cryptotis meridensis, and Nasuella meridensis matched both criteria; hence, additional research is necessary to assess the true taxonomic status and distribution of the remaining species thought to be CM endemics.

About the Authors:

Eliécer E. Gutiérrez

* E-mail: [email protected]

Affiliations Division of Mammals, National Museum of Natural History, Smithsonian Institution, Washington DC, United States of America, Center for Conservation and Evolutionary Genetics, Smithsonian Conservation Biology Institute, National Zoological Park, Washington DC, United States of America

Jesús E. Maldonado

Affiliations Division of Mammals, National Museum of Natural History, Smithsonian Institution, Washington DC, United States of America, Center for Conservation and Evolutionary Genetics, Smithsonian Conservation Biology Institute, National Zoological Park, Washington DC, United States of America

Aleksandar Radosavljevic

Affiliations Plant Biology and Conservation, Northwestern University, Evanston, Illinois, United States of America, Department of Botany, National Museum of Natural History, Smithsonian Institution, Washington DC, United States of America

Jesús Molinari

Affiliation: Departamento de Biología, Facultad de Ciencias, Universidad de Los Andes, Mérida, Venezuela

Bruce D. Patterson

Affiliation: Integrative Research Center, Field Museum of Natural History, Chicago, Illinois, United States of America

Juan M. Martínez-C.

Affiliation: Instituto de Biología, Grupo Mastozoología & Colección Teriológica, Universidad de Antioquia, Medellín, Colombia

Amy R. Rutter

Affiliation: Department of Ecosystem Science and Management, Pennsylvania State University, University Park, Pennsylvania, United States of America

Melissa T. R. Hawkins

Affiliations Division of Mammals, National Museum of Natural History, Smithsonian Institution, Washington DC, United States of America, Center for Conservation and Evolutionary Genetics, Smithsonian Conservation Biology Institute, National Zoological Park, Washington DC, United States of America, Department of Environmental Science & Policy, George Mason University, Fairfax, Virginia, United States of America

Franger J. Garcia

Affiliation: Laboratorio Museo de Zoología, Departamento de Biología, Universidad de Carabobo, Valencia, Carabobo, Venezuela

Kristofer M. Helgen

Affiliation: Division of Mammals, National Museum of Natural History, Smithsonian Institution, Washington DC, United States of America

Introduction

Deer of the northern Andes are among the least studied groups of medium- to large-sized mammals worldwide. They include members of the traditionally recognized but polyphyletic genus Mazama [1, 2, 3], as well as the genera Odocoileus [4] and Pudu [5]. The phylogenetic affinities, biogeography, ecology, behavior and even taxonomy of these deer are poorly understood. The Briceño’s brocket deer, Mazama bricenii, represents one of the best examples of this confusion. Oldfield Thomas [6] named the species after Salomón Briceño Gabaldón, a professional collector who provided him with specimens from the Venezuelan Andes, resulting in the description of several new species of mammals [7], including the single specimen of M. bricenii available for its description. This specimen was collected in the Páramo La Culata, in the Venezuelan Cordillera de Mérida (hereafter referred to as “CM”). This locality is situated at an elevation of 3000 m and is separated from the Cordillera Oriental of Colombia (hereafter referred to as “CO”) by the Táchira Depression, a relatively low (maximum elevation, 960 m), warm, and dry corridor [8]. In the original description, Thomas noted that while among the species of Mazama recognized at that time, M. bricenii was “undoubtedly most nearly allied to Mazama tema” from Central America (= Mazama temama; [9, 10, 11]), another skull from Ecuador more closely reassembled M. bricenii. Thomas mentioned that this similarity with the Ecuadorian specimen “indicates the existence in Ecuador of a highland Brocket allied to, and perhaps identical with, the animal now described from Venezuela” (= M. bricenii). This specimen from Ecuador was likely a Mazama rufina, which is the species of brocket that is currently known to occur in the Ecuadorean Andes [12]. Since the description of M. bricenii, most authors have treated it as a valid species, e.g., [13, 14, 15, 16,17], but some have regarded it either as a subspecies of Mazama rufina, e.g., [18, 19, 20, 21], or as a junior synonym of M. rufina [22]. This discrepancy is not surprising, because no modern taxonomic work has been carried out on Andean Mazama. In addition, the geographic isolation of the CM—to which some authors consider M. bricenii endemic—might have promoted the notion that populations from this cordillera are differentiated enough from those in the Andes of Colombia, Ecuador, and Peru, to merit taxonomic recognition at the species level. In fact, the Táchira Depression has been postulated as a barrier that could have isolated and promoted the differentiation of members of the fauna and flora endemic to higher, cooler, and mesic environments of the CM, e.g., [23, 24, 25, 26, 27, 28].

Regardless of the taxonomic status of Mazama bricenii, information about its distribution is equivocal. Whereas some authors have assumed that M. bricenii is endemic to the Venezuelan Andes [15, 20, 29, 30, 31, 32, 33], others have assumed that the species’ range also includes the CO of Colombia, where it has not been definitively verified [16, 17, 19, 34, 35, 36]. Despite this uncertainty, M. bricenii has been included in the most recent list of mammals present in Colombia [37], whereas it was regarded as a Venezuelan endemic in the most recent list of Venezuelan mammals [33]. Assertions about the presence of M. bricenii in Colombia have not been supported by publicly available, verifiable evidence; however, it is plausible that this deer is present in the Colombian Andes. Dietrich [34] reported a vouchered record from Páramo del Tamá, which, although located within Venezuelan political borders, it is part of the CO of Colombia. Linares [38] recorded the species in the Serranía de Perijá, northwestern Venezuela, close to Colombian territory (i.e., the west versant of the sierra); however, it is unclear whether voucher specimens were deposited in a zoological collection. Critically, all of the aforementioned alleged records of M. bricenii for either the proximity to Colombian territory or for Colombia itself should have relied on information that allows unambiguous distinction between M. bricenii and M. rufina, however, such information has never been published.

Because proper assessments of taxonomic status of populations are essential for efficient conservation planning, e.g., [39, 40], and because firm documentation of a species’ presence within national borders is needed for countries to grant species protection, in the present study, we assess both of these aspects (taxonomy and extent of geographic distribution) for Mazama bricenii. This species is substantially different in pelage, skull morphology, overall body size, and geographic and ecological distribution to all other members of the genus Mazama (as currently understood) except M. rufina [5, 17, 41]. Consequently, we focused in comparisons and analyses of morphological and molecular data between M. bricenii and M. rufina. We first tested the validity of differences in qualitative cranial traits that we observed in preliminary side-by-side comparisons between topotypes of M. bricenii and M. rufina (Fig 1). We determined if these are consistent distinctions that hold taxonomic value across a larger number of specimens and a broader geographic sample (Fig 2). We then conducted both linear morphometric and molecular phylogenetic analyses (mtDNA) to assess the degree of distinctiveness of populations considered by various authors to correspond to M. bricenii. Finally, we employed Ecological Niche Modeling (ENM) analyses based on occurrence records and bioclimatic variables to test whether the Táchira Depression likely represented a barrier isolating Mazama populations in the CM. In addition to clarifying taxonomic boundaries within Mazama, our results prompted us to review the list of mammals currently considered endemic to the CM, and to propose criteria to evaluate putative cases of endemism in the CM.

[Figure omitted. See PDF.]

Fig 1. Possible qualitative cranial differences suggested from preliminary comparisons of topotypes of Mazama bricenii and M. rufina.

Top: topotype of Mazama bricenii from La Culata, Mérida, Venezuela (FMNH 20197); bottom: topotype of M. rufina from Volcán Pichincha, Pichincha, Ecuador (FMNH 44335). Both specimens are adult females. (A) Lacrimal fossa: narrower and substantially deeper in the specimen of M. bricenii than in that of M. rufina; (B) Frontal bones: slightly depressed anteriorly in the specimen of M. bricenii, but relatively straight in that of M. rufina. Illustrations by Megan Krol.

https://doi.org/10.1371/journal.pone.0129113.g001

[Figure omitted. See PDF.]

Fig 2. Linear measurements used in descriptive, univariate, and multivariate statistics.

See Materials and Methods for names and descriptions of measurements.

https://doi.org/10.1371/journal.pone.0129113.g002

Materials and Methods

Source of data

We examined and measured specimens housed in both North and South American collections (S1 File). Geographic data associated with these specimens were pooled together with information from the literature [34, 42] and used for Ecological Niche Modeling (ENM) analyses. The molecular data consisted of eleven sequences of the mitochondrial cytochrome-b gene (CYTB; 1140 bp). Six of them were obtained from degraded DNA extracted from museum specimens, from residual soft tissue attached to skeletons, or from maxilloturbinate bones [43]. For GenBank accession numbers and specimen’s information, see S1 File. The remaining sequences were downloaded from GenBank (accession numbers: JN632657, JN632671, NC020719, NC020721, NC020739). One of these sequences (NC020739) was mistakenly attributed to Pudu mephistophiles by Hassanin et al. [3]; comparisons of this sequence with those obtained from museum specimens of both Pudu mephistophiles and M. rufina allowed us to reidentify the species to which the sequence corresponds as M. rufina (page 40 of [3]) recognized the possibility that their identification was incorrect).

Laboratory methods

In order to avoid contamination from exogenous DNA, we conducted DNA extractions and preparation of reactions previous to PCR amplification in an isolated ancient DNA laboratory located in a separate building from the one containing the primary DNA laboratory and where PCR products of high molecular weight mammalian DNA have never been present. For DNA extractions, we used the method described by Wisely et al. [43], and subsequently employed various combinations of primers to amplify and to sequence short CYTB fragments (S2 File. Primer pairs used for amplification and sequencing of the CYTB gene). These reactions were performed in a six-stage touchdown protocol using a thermal cycler (MJ Research). After an incubation at 95°C for 10 min, the first stage consisted of 2 cycles of the following steps: denaturing at 95°C for 15 seconds, annealing at 60°C for 30 seconds, and extension at 72°C for 1 min. The second, third, fourth, and fifth stages were identical to the first except for lowered annealing temperatures of 58°C, 56°C, 54°C, and 52°C, respectively. The final stage consisted of 40 cycles with an annealing temperature of 50°C. We performed the PCR in 25 μl volumes containing 0.5 U AmpliTaq Polymerase (Applied Biosystems, Foster City, CA), 1X PCR AmpliTaq Buffer, 0.2 μM each dNTP, 0.4 μM of forward and 0.4 μM of reverse primers, 1.5 μM MgCl2, 10X BSA (New England Biolabs, Ipswich, MA), and 50–250 ng of genomic DNA template Successful amplifications were purified using ExoSAP (USB Corporation, Cleveland, OH) incubated at 37°C for 15 min followed by 80°C for 15 min. Both strands of each PCR product were cycle sequenced by subjecting them to a second amplification using a total of 10 μL sequencing reaction mixture, including 50–200 ng of PCR product, 10 pM of corresponding forward or reverse primer, 5X Big Dye Buffer (Applied Biosystems), 1/8 reaction of Big Dye version 3 (Applied Biosystems). The following conditions were used for the Dye Terminator Cycle Sequencing: 25 cycles consisting of denaturing at 96°C for 10 s, annealing at 50°C for 10 s and extension at 60°C for 4 min. These final products were cleaned using Sephadex filtration and then both the 3’ and 5’ strands were sequenced on a 50 cm array using the ABI PRISM 3130 Genetic Analyzer (Applied Biosystems). We employed Geneious v.7.1.5. (Biomatters; http://www.geneious.com/) to compile and edit the sequences that we generated.

Morphological comparisons and morphometric analyses

Specimens used in this study match characteristics that authors have employed to distinguish Mazama rufina and M. bricenii from other brockets, including the presence of a deep lacrimal fossa, black lower legs, a black mask on the head extending from the nose to the nape, and small skull and body size, among others [13, 16, 17, 34, 42]. However, because no information has been published (not even in their original taxonomic descriptions [6, 44] on how to distinguish M. rufina and M. bricenii, we had no option other than assigning species membership to each specimen based on geography (based on previous researchers’ assertions on the species distributions).

In a total of 22 specimens, we scored the state of two qualitative cranial characters observed in preliminary side-by-side comparisons between topotypes of Mazama bricenii (FMNH 20197) and M. rufina (FMNH 44335). These two traits were the relative depth of the lacrimal fossa—shallow in the topotype of M. rufina, but deeply depressed in the specimen of M. bricenii—and the shape of the frontals—relatively straight in the topotype of M. rufina, but depressed anteriorly in the specimen of M. bricenii (Fig 1).

Our morphometric analyses were based on 14 linear measurements of the cranium and mandible, described as follows (Fig 2): interorbital breadth (IB), maximum width across the anterior extremes of the orbits (this measurement is taken with the caliper tips right below the orifice of the lacrimal duct of each orbit); frontal length (FL), length of the frontal at the midline of the skull, taken from the suture formed by the frontal and the nasal (i.e. anteriormost point of the frontal) to the suture of the frontal and the parietal (i.e. anteriormost point of the parietal); intercondylar width (IW), maximum distance between the external borders of the occipital condyles; zygomatic breadth (ZB), greatest distance between the outer margins of the squamosal arms of the zygomatic arches; palatine-premaxillary length (PPL); length from the suture between the palatine and the maxilla (i.e., anteriormost point of the palatine) to the anteriormost point of the premaxilla; basal length (BL), length from the anteriormost point on the lower border of the foramen magnum to the anteriormost point of the premaxilla; condylobasal length (CBL), length from the posteriormost point of the occipital condyles to the anteriormost point of the premaxilla; greatest length of nasals (GLN), length of the nasals at the midline of the skull, taken from the suture formed by the frontal and the nasals (i.e. anteriormost point of the frontal) to the anteriormost point of the nasal bones; maxillary diastema-premaxillary length (MDPL), distance from the anteriormost point of the second upper premolar to the anteriormost point of the premaxilla; upper tooth row length (UTRL), distance from the anteriormost point of the second upper premolar to the posteriormost point of the third upper molar; occipital condyle-premolar length (COPL), length from the posteriormost point of the occipital condyles to the anteriormost point of the second upper premolar; lower tooth row length (LTRL), distance from the anteriormost point of the second lower premolar to the posteriormost point of the third lower molar; notch height (NH), maximum distance from the lowest point of the notch of the mandible to the most ventral point of the angle; and jaw length (JL), maximum distance of the mandible from the most lateral point of the alveolar margin of the canine socket to the posteriormost point of the angle. Measurements were taken to the nearest 0.01 mm with digital calipers.

We conducted univariate and multivariate analyses to assess possible secondary sexual dimorphism and differences among geographic groups. These analyses and descriptive statistics were based only on adult specimens (those with complete permanent dentition); two specimens (CVULA I-2657, I-8559) with marked toothwear, a sign of advanced age, were excluded from the analyses. To evaluate the extent of sexual size dimorphism, we compared the means of each measurement via a 2-tailed two-sample t tests performed in the R Language and Environment for Statistical Computing (hereafter referred to as “R”) [45] with the function t.test. This test assumes that the data are normally distributed; Shapiro-Wilks normality tests conducted for each measurement (with R´s native function shapiro.test) could not reject the null hypothesis of data being normally distributed. For multiple-comparison corrections, we used the simple and conservative Bonferroni correction (hereafter “BoC”) [46] and the more liberal sequential correction proposed by Benjamini and Hochberg [47] (hereafter “BHC”), both calculated with the R package “stats” (function p.adjust) [45]. We carried out two sets of comparisons; the first was based only on data taken from specimens from the CM (6 ♀♀, 2 ♂♂), whereas the second was based on data from all specimens (15 ♀♀, 6 ♂♂). Failure to detect significant differences between sexes (all BoC p values ≥ 0.447; all BHC p values ≥ 0.285) allowed us to use data from both sexes pooled together in subsequent analyses, thus increasing statistical power for these analyses. The latter consisted of three sets of comparisons between measurement means of geographically defined groups, as follows: (1) a comparison between the two putative species based on data from specimens from either the type locality or nearby sites (hereafter referred to as “regional topotypes”) of both species—i.e. specimens of M. bricenii from the CM (n = 7) vs. Ecuadorian specimens of M. rufina from the Pichincha and Cotopaxi mountains (n = 4); (2) a comparison between specimens to the east (putatively representing M. bricenii; n = 7) and west (M. rufina; n = 14) of the Táchira Depression (see Fig 3 for geographic references)—this comparison aimed to discover differences aligned with the hypothesis that the Táchira Depression promoted morphological differentiation in populations to the east of that depression via geographic isolation; (3) a comparison between a group formed by Venezuelan specimens plus Colombian specimens from the CO (n = 13) vs. a group formed by Colombian specimens from the Cordilleras Central and Occidental plus Ecuadorian specimens (n = 8)—this comparison aimed to uncover differences between samples from regions from which previous authors have assigned specimens to either M. bricenii or M. rufina (see Introduction). Descriptive statistics for geographic groups were calculated with the R package psych [48].

[Figure omitted. See PDF.]

Fig 3. Map of Mazama localities.

Circles represent localities; numbers correspond to entries in S1 File. Progressively darker shading indicates areas with elevations of 1000–1500 m (pale gray) and above 1500 m (dark gray). Localities numbers 12 and 25 represent the type localities of Mazama rufina and M. bricenii, respectively.

https://doi.org/10.1371/journal.pone.0129113.g003

We conducted a principal component analysis (PCA) to detect possible morphometric discontinuities among geographic groups. Because multivariate analyses require complete datasets, we estimated missing values (those lacking due to cranial or mandibular damage) with the Probabilistic Principal Component Analysis method as implemented in the R package pcaMethods [49]. After transforming measurement values into natural logarithms, we conducted the PCA based on the correlation matrix using the Paleontological Statistics Software Package for Education and Data Analysis (PAST ver. 3.02) [50].

Phylogenetic analyses and genetic distances

Before conducting phylogenetic analyses, we checked the quality of each DNA sequence, aligned them, and employed the resulting matrix to determine both the most suitable partition scheme and the best fitting model of nucleotide substitutions (for each data subset, if any). All sequences were translated to amino acid sequences using Geneious ver. 7.1.5 and examined to assure that no premature stop codons were present. Subsequently, sequences were aligned using default options in MAFFT ver. 7.017 [51] as implemented in Geneious. Partition schemes and model of nucleotide substitution for each data subsets were selected with PartitionFinder ver. 1.0.1 [52] using the corrected Bayesian Information Criterion (BIC) and considering models applicable in MrBayes [53].

We conducted phylogenetic analyses using maximum likelihood (ML) and Bayesian inference (BI) optimality criteria. For all analyses, we used one sequence of Odocoileus cariacou (name follows the provisional taxonomic arrangement for Odocoileus proposed by Molina and Molinari [4], and Molinari [25]) and two of Mazama americana as outgroups (GenBank accession numbers JN632671, NC020719, JN632657). According to a previous study, these two species are closely related to M. rufina [3]. The ML analysis was conducted with 20 independent searches in the Genetic Algorithm for Rapid Likelihood Inference (GARLI 2.0) [54] using the default settings. The Bayesian inference analysis was conducted in MrBayes v. 3. 2 [53]. The search started with a random tree, and the Markov chains were run for 100,000,000 generations; trees were sampled every 1000 generations. Default values were kept for the ‘‘relburnin” and ‘‘burninfrac” options in MrBayes (i.e. relburnin = yes; burninfrac = 0.25); therefore, the first 25,000,000 generations (25,000 trees) were discarded as burn-in. To assess nodal support, posterior probability estimates were obtained based on the remaining (75,000) trees. For the ML analysis, nodal support assessed with nonparametric bootstrapping [55]. This analysis was also conducted with GARLI and was based on 1000 searches (100 pseudoreplicated data matrices, and 10 searches for each of them). The degree of support received by individual nodes in the ML bootstrap analysis was categorized as follows: strong if bootstrap value ≥ 75%; moderate if bootstrap value > 50% and < 75%; negligible if value ≤ 50%. For the BI analysis, the categorization of nodal support is as follows: strong if posterior probability ≥ 0.95; negligible support if posterior probability values < 0.95.

A high degree of sequence divergence is neither necessary nor sufficient for species recognition [56, 57, 58]; however, pairwise genetic distances provide a heuristic basis for comparisons of genetic variation within and among lineages [59]. Thus, we calculated the average uncorrected (p) distance within each haplogroup and the average pairwise p distances among them. We also report the commonly used Kimura 2-parameter-corrected distances to facilitate comparisons with data from the literature for non-volant, terrestrial mammals. Genetic distances were calculated using MEGA ver. 5.2.1 [60].

Ecological Niche Modeling

We constructed Ecological Niche Models (ENMs) using the maximum entropy algorithm Maxent [61] to assess whether the Táchira Depression represents a current or past climatic barrier to dispersal. To accomplish this, we trained preliminary models using an appropriate study region (see below) and then projected them onto a larger study region twice, once while employing current climatic conditions and the second time using climatic conditions estimated for the Last Glacial Maximum (LGM).

Occurrence records were gathered from our own examination of museum specimens and from the literature [34, 42], and were georeferenced with maps, gazetteers, and digital resources (S1 File). Because we did not have a priori information regarding which occurrence records corresponded to Mazama bricenii and which to M. rufina, or, alternatively, if all records would actually belong to a single species, we first applied the methods described above to determine the taxonomic status of M. bricenii. To mitigate potential effects of sampling bias and spatial autocorrelation [62], we employed the R package spThin [63] to spatially filter occurrence records to have a minimum distance of at least 10 km among them (i.e. so that no record was closer than 10 km to any other record). After filtering, 22 occurrence records from unique localities remained for modeling.

The environmental data used to create the models consisted of 19 bioclimatic variables from WorldClim 1.4 [64] that previous studies found to be important in determining mammal species distributions, e.g., [65, 66, 67]. To define an appropriate study region for model calibration, we followed the criteria proposed by Anderson and Raza [65], and the operational strategy used by Gutiérrez et al. [67]. Thus, we employed ArcGIS 10.2 to first create a minimum convex hull that enclosed the filtered occurrence records, and then buffered this at a distance of 50 km. As mentioned above, the models were then projected to a larger study region (extend: 7.0° S–13.5° N and 60.0° W–81.5° W) under both current climatic conditions and also using estimated paleoclimatic data for the LGM from two general circulation models (GCMs): the Community Climate System Model (CCSM4) [68] and the Model for Interdisciplinary Research on Climate (MIROC-ESM) [69]. The environmental dataset corresponding to present-day climatic conditions was at a resolution of 30 seconds, whereas those corresponding to LGM climatic conditions had coarser resolution of 2.5 minutes.

Ecological Niche Models were created using Maxent ver. 3.3.3h [61]. Since several recent studies have demonstrated the importance of selecting model Maxent settings carefully to balance model’s performance and complexity, we tuned the value of the regularization multiplier and determined the optimal selection of feature classes [66, 70, 71, 72]. All tuning experiments were implemented using the R package ENMeval [73]. The regularization multiplier was varied from 0.5 to 6.0 in increments of 0.5, and the following four feature classes (or combinations thereof) were tested: (1) linear; (2) linear and quadratic; (3) hinge; and (4) linear, quadratic, and hinge. ENMeval allows for several data-partitioning schemes; for this study we employed the “checkerboard1” approach, a variation on the ‘masked geographically structured’ data-partitioning strategy described in Radosavljevic and Anderson [71]. Model performance was assessed using the Akaike Information Criterion corrected for small sample sizes (AICc) [70, 72]. The final model was constructed with the combination of regularization multiplier and feature classes that yielded the lowest value of AICc, and employing all filtered, georeferenced occurrence records. As AICc can only select the ‘best’ model from among a set of models, and does not directly assess model performance, we also inspected omission rate and test AUC of the models selected as optimal. Additionally, since we were projecting these models onto a different region than that used for calibrating the model, we inspected the multivariate environmental similarity surfaces and clamping maps produced by Maxent to determine whether any environmental variables on the larger region were outside the range of climatic conditions present in the calibration study regions.

Final models were converted into binary maps of “suitable” vs “unsuitable” habitat using the minimum training presence threshold calculated by Maxent. This threshold considers as suitable all pixels that have a suitability score equal to or greater than the lowest suitability predicted for a known occurrence. While this threshold value is likely to include some habitat that may be marginally suitable (due to real and artifactual causes; [71, 74]), for our purposes it is an objective and logical choice—our goal is to test whether a warm and dry area (i.e., the Táchira Depression) acts as a barrier to dispersal, thus even marginal habitat should be identified.

Results

Morphological comparisons and morphometric analyses

Two differences in qualitative cranial traits observed between Mazama bricenii and M. rufina in initial comparisons that included topotypes of both species—i.e. depth of lacrimal fossa and shape of frontal bones (Fig 1)—proved to be extremely variable, and do not hold taxonomic value. In a total of 22 specimens examined, we found (a) deep lacrimal fossae and frontal bones posteriorly depressed in eight females, five from the CM, two from the CO, and one from Ecuador; (b) shallow lacrimal fossae and straight frontal bones in two females from Ecuador; (c) deep lacrimal fossae but straight frontal bones in three specimens, one female from CM and one male and one female from the Cordillera Central of Colombia. In addition, depth of the lacrimal fossa in two specimens was scored as intermediate between “deep” and “shallow” (see Fig 1), and the shape of the frontal bones in nine specimens was intermediate between “depressed” and “straight”. These observations do not reveal any consistent pattern of geographic variation.

We did not detect significant differences in size in comparisons between topotypes of Mazama bricenii and M. rufina, or in comparisons between specimens separated by the Táchira Depression (in all of these analyses BoC p values ≥ 0.469 and BHC p values ≥ 0.219). Except for two measurements, comparisons between a group formed by specimens from Venezuela and the CO of Colombia—where previous authors have reported deer allegedly identified as M. bricenii—as compared with a group formed by specimens of M. rufina from Ecuador and from the Cordillera Central of Colombia failed to detect significant differences (p values ≥ 0.116 for both BoC and BHC). However, the former group has significantly narrower skulls (mean of zygomatic arches = 71.06 mm; p values for both BoC and BHC for t-tests ≤ 0.032) and marginally significantly shorter jaws (mean of jaw length = 120.99 mm; p values for the BHC for t-tests = 0.070) than specimens from Ecuador and the Cordillera Central of Colombia (75.90 mm and 129.75 mm, respectively). Descriptive statistics for geographic groups are shown in Tables 1 and 2.

[Figure omitted. See PDF.]

Table 1. Measurements of Mazama specimens from Venezuela and the Cordillera Oriental of Colombia.

https://doi.org/10.1371/journal.pone.0129113.t001

[Figure omitted. See PDF.]

Table 2. Measurements of Mazama specimens from the Cordillera Central of Colombia and Ecuador.

https://doi.org/10.1371/journal.pone.0129113.t002

The principal components analysis did not reveal morphometric discontinuities among geographic groups. The first and second components explained 49.99% and 15.90% of the variance, respectively. Because the sign of the loadings of the first component (PC1) are all the same, we interpret PC1 as an axis that captured primarily differences in size among specimens (Table 3). Seven measurements (ZB, PPL, BL, CBL, GLN, DLM, JL) had the largest loadings in PC1. Some loadings on the second component (PC2) are negative and others positive; hence, we interpret PC2 as an axis that primarily captured differences in cranial proportions among specimens. In this component, the measurements with the largest loadings were GLN and DLM (both loading negatively) and UTRL and LTRL (both loading positively). A scatter plot constructed with specimens’ scores on the first two components shows substantial overlap among geographic groups in both components (Fig 4). Specimens from the CM, including the type locality of Mazama bricenii, are widely scattered across PC1 and PC2, overlapping with specimens from the Colombian and Ecuadorian Andes.

[Figure omitted. See PDF.]

Table 3. Results of the principal component analysis based on measurements of specimens of Mazama from the northern Andes.

https://doi.org/10.1371/journal.pone.0129113.t003

[Figure omitted. See PDF.]

Fig 4. Plot of specimen scores on the first two axes of the principal components analysis of skull measurements of adult specimens of Mazama.

Solid and open symbols represent male and female specimens, respectively. Geographic provenance represented as follows: black circles: Cordillera de Mérida (Venezuela); red circles: Cordillera Oriental (Colombia); red triangles: Cordillera Central (Colombia); blue triangles: Ecuador. PC1 is a size axis in which larger specimens appear toward the right side of the axis, whereas PC2 represents differences in cranial proportions (Table 3).

https://doi.org/10.1371/journal.pone.0129113.g004

Molecular data and models of molecular evolution

Despite reliance on degraded DNA from museum specimens, our sequence matrix had only ca. 7% missing data (i.e., entries coded as unknown). PartitionFinder found the best partitioning scheme to be one formed by two subsets, one subset containing bases in the first and second codon positions, and another subset containing bases in the third codon position. The structure of the model for both subsets corresponds to the Hasegawa-Kishino-Yano model [75] with a proportion of invariant sites (i.e., HKY+I). The fact that PartitionFinder assessed that the best fit-model was the same for both subsets but still kept them separate means that the parameter values for these two subsets were sufficiently different that they were best modeled separately.

Phylogenetic analyses and genetic distances

Our analyses of CYTB sequences did not uncover phylogeographic patterns nor did they establish the phylogenetic distinction of Mazama bricenii. After examining the average standard deviation of split of our BI analysis and confirming that it reached stationarity, we compared its topology with that resulting from the ML analysis. Both trees recovered the same branching patterns with nearly the same degree of nodal support. The two samples from the CM (one topotype and one near-topotype of M. bricenii) were recovered as sister to each other with strong nodal support (haplogroup A in Fig 5), but embedded within the diversity of M. rufina. Sister to this haplogroup (A) was another (B) formed by samples from unknown localities within Colombia (reported by Hassanin [3]), one sample from the CO of Colombia, and the sequence of a topotype of M. rufina (from Pichincha, Ecuador). This latter haplogroup received moderate and strong support in the ML and BI analyses, respectively, whereas the sister relationship between haplogroups A and B received non-neglible support only in the ML analysis. One sequence from Ecuador and another from Colombia were recovered as successive sisters to the clade formed by haplogroups A and B (always with either moderate or strong support; Fig 5).

[Figure omitted. See PDF.]

Fig 5. Phylogenetic tree of cytochrome-b sequences of Mazama from the northern Andes.

This is the best topology resulting from the maximum-likelihood analysis. Nodal support is indicated at each node, except when the involved relationship received negligible support. Bootstrap values (from the maximum likelihood analysis) and posterior probabilities (from the Bayesian inference analysis) are indicated before and after the slash (“/”). Three topotypes (one M. “bricenii” and two M. rufina) are indicated with bold type (see detailed locality information in S1 File). The length of each sequence (number of base pairs, bp) is indicated at each terminal label. Asterisks denote sequences obtained from DNA extracted from museum specimens; all other sequences were downloaded from GenBank.

https://doi.org/10.1371/journal.pone.0129113.g005

None of the genetic distances calculated were particularly high. The mean p-distances within haplogroups A (from the CM; Fig 5) and B were 2.0% and 1.7%, respectively. The distances between haplogroups A and B were 3.0% and 3.1% for the p- and K2P-corrected metrics, respectively. With respect to a group containing all sequences from specimens from west of the Táchira Depression, haplogroup A showed distances of 2.6% and 2.7% for the p- and K2P-corrected metrics, respectively; the mean p-distance within the former group was 1.8%.

Ecological Niche Modeling and Climatic Suitability in the Táchira Depression

Ecological Niche Models (ENMs) indicate that the Táchira Depression does not represent a current climatic barrier to dispersal of Mazama rufina and did not during the cooler, drier climates of the Last Glacial Maximum (LGM; Fig 6; for results of tuning see S1 Table. Results from tuning experiments using ENMeval). At the Minimum Training Presence (MTP) threshold the model projected to current climatic conditions indicated suitable conditions in the entire Táchira Depression. Under past conditions, suitable conditions existed across on an even more extensive area, including sites of much lower elevation. Both model projections onto datasets of LGM climate conditions (CCSM in Fig 6; MIROC-ESM in S1 Fig Maxent models of abiotically suitable areas for Mazama rufina projected onto estimated climatic conditions of the Last Glacial Maximum.) indicated that the Táchira Depression contained suitable conditions for M. rufina at that time; however, the two projections differed in their overall predictions. While there is still much work to be done on assessing the effects of alternate paleo-climate reconstructions on ENMs, at least one study has shown that models based on CCSM tend to agree more strongly with independent data [76]. Inspections of clamping maps and MESS surfaces did not indicate that any climatic conditions in the projection region (and time period) were outside of the range of conditions used to calibrate the models.

[Figure omitted. See PDF.]

Fig 6. Abiotically suitable areas for Mazama rufina as predicted by ecological niche modeling analyses.

In green are areas predicted to be suitable under current climate conditions (A, B) and during the Last Glacial Maximum (C, D). Dotted line in A circumscribes the region shown in the close up panels (B, D) and contains the Táchira Depression. Both sets of models indicate extensive suitable conditions in the area of the Táchira Depression, suggesting long term and continuous habitat connectivity between the Cordillera Oriental and Cordillera de Mérida. The Minimum Training Presence threshold of Maxent was used to convert continuous values of predicted suitability into a binary prediction, which classifies each pixel of the image (map) as suitable or unsuitable.

https://doi.org/10.1371/journal.pone.0129113.g006

Discussion

Mazama bricenii is not a valid taxon

Results of our analyses based on morphological and molecular data do not support recognition of Mazama bricenii as a valid taxon at either the species or subspecies level, and should be regarded as a junior synonym of Mazama rufina. Specimens from the CM cannot be distinguished from the rest of the specimens based on cranial morphology. Principal component analysis (PCA) is an ordination technique that allows detection of morphometric discontinuities in size, proportion, or both among taxa or geographic groups in which such discontinuities exist, e.g., [77, 78, 79, 80, 81]; however, the wide overlap in specimen scores on the first two components evince lack of morphometric discontinuity. In addition, the fact that the qualitative cranial traits that we evaluated did not show consistent differences among geographic groups is congruent with our conclusion that specimens from regions where previous authors have identified specimens as “M. bricenii” are cranially indistinguishable from Mazama rufina. This lack of distinctiveness is in contrast to the cases of the only two other known medium to large-sized mammals considered endemic to the CM, the Mérida mountain coati, Nasuella meridensis [79], and the Mérida páramo deer, Odocoileus lasiotis (see page 59 of [25]), both of which show substantial cranial differences with respect to their congeners in the Colombian Andes.

Phylogenetic analyses of molecular data also failed to support the validity of Mazama bricenii. Because we used a fast-evolving marker, if samples from the CM—or from both the CM and the CO—represented a taxon distinct from M. rufina, then we might expect them to be recovered as a reciprocally monophyletic haplogroup with respect to M. rufina. However, the two sequences from the CM were embedded within the shallow diversity of M. rufina, rendering that taxon paraphyletic. It is likely that future studies using markers with even faster mutation rates (e.g., the mitochondrial control region) could recover the samples from the CM—or from the CM and the CO—as a reciprocally monophyletic haplogroup with respect to M. rufina; however, even considering that potential topology, recognition of “bricenii” at the species level would be compromised by its lack of morphological diagnosabilty [82].

The genetic distances between geographically defined groups attributable to Mazama bricenii and M. rufina (ca. 3%) scarcely exceeds within-group variation (ca. 2%). This pattern does not suggest the presence of more than a single species among analyzed populations. Very few studies have reported genetic distances calculated from CYTB sequences between putative sister species of deer. Several of these studies instead reported within-group distance ranges for all of their focal species rather than a matrix of pairwise comparisons. Despite scarcity of data to compare with, it seems clear that CYTB genetic distances between recovered haplogroups A and B (Fig 5) are not particularly high (ca. 3%) relative to reported ranges of interspecific distances in other cervids, e.g. among various Neotropical cervids (including taxa in multiple genera; 8.6–13.9% for K2P-corrected distances) [2], among various species of muntjac, Muntiacus (6.5–8.7% K2P-corrected distances) [83], or between the closely related species Cervus elaphus and C. canadensis (5.7% uncorrected p-distance) [84].

In summary, our results from analyses of morphological and molecular data demonstrate that populations previously referred by authors as “Mazama bricenii” or “Mazama rufina bricenii” do not merit taxonomic recognition, and populations previously referred to by these names should be regarded as Mazama rufina. Although sufficient evidence supporting the validity of M. bricenii has never been published, the long-lived notion that populations from the CM were a valid species endemic to that cordillera was widely accepted simply because the Táchira Depression has been regarded as an important barrier to dispersal between the CO of Colombia and the CM. Statements like the following, by influential mammalogist J. A. Allen (page 529 of [13]), illustrate this way of thinking:

“For example, Mazama rufina of Mount Pichincha in Ecuador and M. brincenii of the paramo of the Sierra de Merida in Venezuela so closely resemble each other in size, in coloration, and in the peculiar character of the pelage, that if their known ranges were contiguous they would naturally be regarded as local forms of a single species, but their wide separation by regions of much lower elevation and very different climatic conditions renders improbable any continuous distribution and consequent geographical intergradation.”

Biogeographic significance of the Táchira Depression and mammalian endemism in the CM

The warm and dry climatic conditions present in the Táchira Depression now separate cooler and more mesic habitat types (cloud forest and páramo) in the CM from similar habitats in the CO of Colombia. Because of its distinctive climate and large area, the Táchira Depression could indeed represent an important barrier to dispersal for species with lower vagility and strictly restricted to cloud forest, páramo, or both. The suspicion that the Táchira Depression represented a barrier for dispersal of Mazama “bricenii” (= M. rufina) very likely biased the view that regarded it as a valid species, implying the assumption that the “species” may had differentiated in isolation in the CM. Nevertheless, specimen-based research has shown that the Táchira Depression should not be assumed to be a barrier for species not restricted to cool and mesic habitats typical of high elevations of the northern Andes. Examples include a heteromyid rodent (Heteromys australis) and a didelphid marsupial (Marmosa waterhousei), which were recorded in the CM until recently [85, 86]. Both species were known to occur in the Colombian CO [87, 88], but the current climatic conditions of the Táchira Depression made their presence in the CM highly unexpected. Although H. australis inhabits very mesic evergreen forests up to ca. 2500 m in elevation, and M. waterhousei inhabits humid lowland and mountain forests from 50 m to 1100 m, neither of these species is restricted to habitat types with low temperatures [85, 88, 89]. As previously suggested [85, 86], finding these species in the CM indicates that the Táchira Depression was not always a barrier for them. During the LGM, about 26,500–20,000 years before present [90], the altitudinal zonation of mountains was affected with the descent of upper vegetation belts to lower elevations. This phenomenon could have led to a connection of cool and mesic habitat between the CO and the CM over the area currently occupied by the warm and dry Táchira Depression [85, 86]. Congruent with this possibility, the projection of the ENM of M. rufina onto estimated climatic conditions of the Táchira Depression during the LGM revealed that suitable climatic conditions for the species existed then. The ENM analyses show that, even at present, suitable climatic conditions (at least minimally) exist for M. rufina in the Táchira Depression (Fig 6). Similarly, results of ENM analyses in a previous study predicted suitable climatic conditions at present in the Táchira Depression for another species, Nephelomys meridensis [65], supposedly endemic to the CM (cf [91]).

The most recent isolation of cool and mesic habitat in the CM should not have initiated until the beginning of the current interglacial period ca. 15,000–12,000 years before present [92, 93], or perhaps even more recently. Thus, even for species restricted to cloud forest and páramo, the period of time of isolation may have been too short for populations to differentiate. In the case of Mazama rufina, suitable climatic conditions are still present in the Táchira Depression, therefore, populations in the CM may have never become fully isolated. This hypothesis may apply also to other members of the cordillera’s biota. To explore this possibility, we revisited the list of mammals that have been considered endemic to the CM [30]. We first corrected previous omissions of species’ records reported in the literature, and then updated the list of endemics based on recently reported records and taxonomic changes. Four species were implicated in this verification process: one didelphid marsupial (Gracilinanus dryas), one echimyid rodent (Olallamys edax), and two sigmodontine rodents (Nephelomys meridensis and Neusticomys mussoi). All of these species occur in Colombia as well as in the CM as documented in previously overlooked or unavailable literature [91, 94, 95, 96]. If the original colonization of the CM by these species, or their ancestors, was from the CO—which seems plausible given the pattern of nestedness observed in other mountain systems in the region [97]—then the fact that these species are not endemic to the CM signifies that Táchira Depression did not represent a barrier to their dispersal and distribution, as previously implied [30]. In addition, another sigmodontine rodent, Oecomys flavicans, previously considered to be endemic to the CM, has recently been reported for the non-Andean Serranía de San Luis in northwestern Venezuela [97]. One addition to the list of CM endemics is the procyonid Nasuella meridensis, which was not recognized as a valid species until recently [26]. Our overview of literature identified six species currently considered endemic to the CM. These include one soricid, Cryptotis meridensis [28, 98], one procyonid, N. meridensis [26], one cervid, Odocoileus lasiotis [4, 25], and three sigmodontine rodents, Aepeomys reigi [99], Thomasomys vestitus [100, 101, 102], and an undescribed sigmodontine rodent of the genus Nephelomys [103, 104].

Results from our ENM analyses, and the fact that several species once thought to be endemic to the CM are actually present in the CO, indicate that the Táchira Depression should not be invoked a priori as a barrier that could have promoted isolation and subsequent differentiation of taxa in the CM, unless the following two criteria are met: logically, that (1) evidence (e.g. data indicating morphological and/or molecular distinctiveness) exists suggesting that the putative taxon endemic to the CM is distinctive from congeneric populations from other regions, particularly from the nearby Andes of Colombia, including the Venezuelan portion of the CO (El Tamá); and that (2) the putative endemic taxon is currently restricted to either cloud forest, páramo, or both. The latter is proposed because mammals present in habitats with either warm or dry (or both) climatic conditions have been consistently found both to the east and west of the Táchira Depression (i.e., present in the CM as well as in the CO).

Of the six mammal species currently known only from the CM, most do not satisfy these criteria. The unnamed species of Nephelomys does not match the second criterion (i.e. it is not restricted to cloud forest or páramo). This species is only known from a single locality in seasonal forest at an elevation of 1100 m [104] nearby Pregonero, a town in southwestern CM. The remaining five species are only known from cloud forest or páramo habitat, or both; however, two of them, Thomasomys vestitus and Odocoileus lasiotis, remain to be compellingly compared to congeners in the Colombian Andes, especially in the CO. In the case of T. vestitus, karyological data have been gathered, but did not reveal substantial differences with respect to several other species of Thomasomys that would otherwise support its taxonomic status. In the case of O. lasiotis, its taxonomic status has been assessed by comparison of qualitative and quantitative morphological data with respect to other Venezuelan and North American populations of Odocoileus, but similar comparisons with respect to populations from the Colombian Andes remain to be conducted. It is noteworthy, however, that examination of a small number of specimens revealed substantial differences in pelage and cranial morphology between O. lasiotis and specimens from the Colombian Andes, and a plausible hypothesis for the isolation of O. lasiotis in the páramos of the CM has been postulated [25]. Finally, three species, Aepeomys reigi, Cryptotis meridensis, and Nasuella meridensis, match both of our proposed criteria for regarding a taxon as endemic to the CM with confidence. These species are currently restricted to cloud forest and páramo above ca. 1600 m in elevation—in the case of A. reigi between 1600 and 3230 m [99], C. meridensis 1670–3950 m, and N. meridensis 1980–4000 m [26, 28, 38, 98, 105, 106, 107]. Thus, our review suggests that the CM is likely home to true mammalian endemism, comprising species restricted to highest elevation mesic habitats, but that the level of endemism in these mountains is probably not as substantial as previous reviews have suggested. With the exception of three species, all other species currently considered endemic to the CM may be better regarded more tentatively as putative endemics until further studies are undertaken to assess their taxonomic status and distribution, both of which are essential components for planning for their effective of conservation.

Supporting Information

[Figure omitted. See PDF.]

S1 Fig. Maxent models of abiotically suitable areas for Mazama rufina projected onto estimated climatic conditions of the Last Glacial Maximum.

A: projection onto the Community Climate System Model dataset (CCSM4); B: projection onto the Model for Interdisciplinary Research on Climate dataset (MIROC). While in both cases suitable conditions are predicted in the Táchira Depression, the shape and strength of the prediction varies.

https://doi.org/10.1371/journal.pone.0129113.s001

(TIF)

S1 File. Gazetteer, specimens examined, and GenBank accession numbers.

https://doi.org/10.1371/journal.pone.0129113.s002

(DOCX)

S2 File. Primer pairs used for amplification and sequencing of the CYTB gene.

https://doi.org/10.1371/journal.pone.0129113.s003

(DOCX)

S1 Table. Results from tuning experiments using ENMeval.

https://doi.org/10.1371/journal.pone.0129113.s004

(DOCX)

Acknowledgments

This study results from a larger initiative by EEG, JEM, and KMH to investigate the systematics and biogeography of Neotropical deer, and for which we have received the support of a number of colleagues and institutions. We would like to thank the curators and staff of institutions that allowed access to voucher material, especially Nancy Simmons, Robert Voss, Eileen Westwig, and Neil Duncan (AMNH); Hugo López-Arevalo and Catalina Cárdenas (ICN); Pascual Soriano (CVULA); Francisco Bisbal-E. and Javier Sánchez-H. (EBRG); John Phelps and William Stanley (FMNH); Carmen Ferreira, Sandra Giner, and Mercedes Salazar (MBUCV); Darrin Lunde, Nicole Edmison, Esther Langan, and Suzanne Peurach (NMNH/USNM). In addition, Santiago Burneo and María A. Camacho provided access to specimens of Mazama deposited at the Museo de Zoología de la Pontificia Universidad Católica del Ecuador. We are thankful to Megan Kroll for the illustration of specimens shown in Fig 1, to Paul Velazco for sending us photographs of some of the specimens reported here and providing logistic support during EEG visits to the AMNH, to Marcus Litchfield for assisting EEG during part of the laboratory work, and to Franziska Bauer, who provided helpful comments to EEG with regard to generation of molecular data from museum specimens. Manuel Ruiz-García kindly provided information about the locality of origin of a sequence we obtained from Genbank. We thank Robert Fleischer and Nancy McInerney for assistance and support in the genetics laboratory at the CCEG.

Author Contributions

Conceived and designed the experiments: EEG JEM KMH. Performed the experiments: ARR BDP EEG FJG JEM JM JMM MTRH. Analyzed the data: AR EEG. Contributed reagents/materials/analysis tools: EEG JEM KMH. Wrote the paper: EEG. Revised, commented, and edited earlier versions of the manuscript and contributed intellectually: EEG JEM AR JM BDP JMM ARR MTRH FJG KMH.

Citation: Gutiérrez EE, Maldonado JE, Radosavljevic A, Molinari J, Patterson BD, Martínez-C. JM, et al. (2015) The Taxonomic Status of Mazama bricenii and the Significance of the Táchira Depression for Mammalian Endemism in the Cordillera de Mérida, Venezuela. PLoS ONE 10(6): e0129113. https://doi.org/10.1371/journal.pone.0129113

References

1. Gilbert C, Ropiquet A, Hassanin A. Mitochondrial and nuclear phylogenies of Cervidae (Mammalia, Ruminantia): systematics, morphology, and biogeography. Mol Phylogenet Evol. 2006;40: 101–117. pmid:16584894

2. Duarte JMB, González S, Maldonado JE. The surprising evolutionary history of South American deer. Mol Phylogenet Evol. 2008;49: 17–22. pmid:18675919

3. Hassanin A, Delsuc F, Ropiquet A, Hammere C, van Vuuren BJ, Matthee C, et al. Pattern and timing of diversification of Cetartiodactyla (Mammalia, Laurasiatheria), as revealed by a comprehensive analysis of mitochondrial genomes. C R Biol. 2012;335: 32–50. pmid:22226162

4. Molina M, Molinari J. Taxonomy of Venezuelan white-tailed deer (Odocoileus, Cervidae, Mammalia), based on cranial and mandibular traits. Canadian J Zool. 1999;77: 632–645.

5. Hershkovitz P. Neotropical deer (Cervidae) part I. Pudus, genus Pudu Gray. Fieldiana Zool. 1982;1330: 1–86.

6. Thomas O. A new deer of the brocket group from Venezuela. Ann Mag Nat Hist. (series 8). 1908; 1: 349–350.

7. Beolens B, Watkins M, Grayson M. The eponym dictionary of mammals. Baltimore: Johns Hopkins University Press; 2009.

8. van der Hammen TYA, Cleef M. Development of the high Andean paramo flora and vegetation. In: Vuilleumier F, Monasterio M, editors. High altitude tropical biogeography. New York: Oxford University Press; 1986. pp. 153–201.

9. Tate GHH. The mammals of the Guiana region. Bull Amer Mus Nat Hist. 1939;LXXVI: 151–229.

10. Jorge W, Benirschke K. Centromeric heterochromatin and G-banding of the red brocket deer Mazama americana temama (Cervoidea, Artiodactyla) with probable non-Robertsonian translocation. Cytologia 1977;42: 711–721. pmid:923292

11. Groves CP, Grubb P. Relationships of living deer. In: Wemmer CM, editor. Biology and management of the Cervidae. Washington DC: Smithsonian Institution Press; 1987. pp. 21–59.

12. Tirira D. Libro rojo del los mamíferos del Ecuador. Vol. 1. Quito: SIMBIOE, EcoCiencia, Ministerio del Ambiente, and UICN. 2001.

13. Allen JA. Notes on American deer of the genus Mazama. Bull Amer Mus Nat Hist. 1915;34: 521–553.

14. Lydekker R. Artiodactyla, Families Cervidae (Deer), Tragulidae (Chevrotains), Camelidae (Camels and Llamas), Suidae (Pigs and Peccaries), and Hippopotamidae (Hippopotamuses). Vol. 4 of Catalogue of the Ungulate Mammals in the British Museum. London: British Museum (Natural History). 1915.

15. Geist V. Deer of the World: their Evolution, Behavior, and Ecology. Mechanicsburg: Stackpole Books; 1998.

16. Lizcano J, Yerena E, Alvarez SJ, Dietrich JR. Mérida Brocket Mazama bricenii (Thomas 1908). In: Barbanti JMD, Gonzalez S, editors. Neotropical cervidology: biology and medicine of Latin Amercian deer. Jaboticabal: Funep-IUCN; 2010. pp. 181–184.

17. Groves C, Grubb P. Ungulate Taxonomy. Baltimore: Johns Hopkins University Press; 2011.

18. Cabrera A. Catálogo de los mamíferos de América del Sur [part 2]. Rev Mus Argent Cienc Nat “Bernardino Rivadavia”, Cienc Zool. 1961;4: 309–732.

19. Bisbal FJ. Distribución y taxonomía del venado Matacán (Mazama sp) en Venezuela. Acta Biol Venezuelica. 1991;13: 89–104.

20. Whitehead GK. The whitehead encyclopedia of deer. Shrewsbury: Swan Hill Press; 1993.

21. Wemmer C, Mccarthy A, Blouch R, Moore D, editors. Deer: status survey and conservation action plan. Cambridge: IUCN/SSC Deer specialist group; 1998.

22. Rumiz DI, Pardo E, Eulert CF, Arispe R, Wallace RB, Gómez H, Ríos-Uzeda B. New records and a status assessment of a rare dwarf brocket deer from the montane forests of Bolivia. J Zool. 2007;271: 428–436.

23. Vuilleumier F, Ewert DN. The distribution of birds in the Venezuelan paramos. Bull Amer Mus Nat Hist. 1978;162: 47–90.

24. Cracraft J. Historical biogeography and patterns of differentiation within the South American avifauna: areas of endemism. Ornithological Monograph. 1985;36: 49–84.

25. Molinari J. Variación geográfica en los venados de cola blanca (Cervidae, Odocoileus) de Venezuela, con énfasis en O. margaritae, la especie enana de la Isla de Margarita. Mem Fund La Salle Cienc Nat. 2007;167: 29–72.

26. Helgen KM, Kays R, Helgen LE, Tsuchiya-Jerep MTN, Pinto CM, Klaus-Peter Koepfli K-P, et al. Taxonomic boundaries and geographic distributions revealed by an integrative systematic overview of the mountain coatis, Nasuella (Carnivora: Procyonidae). Small Carnivore Conservat. 2009;41: 65–74.

27. Castaño JH, Soriano PJ, Ataroff M. Características físicas y ambientales del Ramal Calderas, Andes de Venezuela. In: Rial-B A, Señaris JC, Lasso CA, Flores AL, editors, Evaluación rápida de la biodiversidad y aspectos socioecosistémicos del Ramal Calderas, Andes de Venezuela. Boletín RAP de Evaluación Biológica N° 56, Conservation International; 2010. pp. 45–49.

28. Quiroga-Carmona M, Molinari J. Description of a new shrew of the genus Cryptotis (Mammalia: Soricomorpha: Soricidae) from the Sierra de Aroa, an isolated mountain range in northwestern Venezuela, with remarks on biogeography and conservation. Zootaxa. 2012;3441: 1–20.

29. Czernay S. Spiesshirsche und Pudus. Die Neue Brehm Bucherei. 1987;581: 1–84.

30. Soriano P, Díaz de Pascual A, Ochoa-G J, Aguilera M. Biogeographic analysis of the mammal communities in the Venezuelan Andes. Interciencia. 1999;24: 17–25.

31. Weber M, Gonzalez S. Latin American deer diversity and conservation: A review of status and distribution. Ecoscience. 2003;10(4): 443–454.

32. Grubb P. Order Artiodactyla. In: Wilson DE, Reeder DM, editors. Mammal species of the world. 3rd ed. Baltimore: Johns Hopkins University Press. 2005. pp 629–722.

33. Sánchez-H J, Lew D. Lista actualizada y comentada de los mamíferos de Venezuela. Mem Fund La Salle Cienc Nat. 2012 [“2010”];173–174: 173–238.

34. Dietrich JR. Biology of the brocket deer (genus Mazama) in northern Venezuela. Ph.D. dissertation, the University of Basel. 1993.

35. Trolle M, Emmons LH. A record of dwarf brocket from lowland Madre de Dios, Peru. Deer Specialist Group News 2004;19: 2–5.

36. Mattioli S. Family Cervidae (Deer). In: Wilson DE, Mittermeier RA, editors. Handbook of the Mammals of the World. Volume 2. Hoofed Mammals. Barcelona: Lynx Editions; 2011. pp. 350–443.

37. Solari S, Muñoz-Saba Y, Rodríguez-Mahecha JV, Defler TR, Ramírez-Chaves HE, et al. Riqueza, endemismo y conservación de los mamíferos de Colombia. Mastozool Neotrop. 2013;20: 301–365.

38. Linares OJ. Mamíferos de Venezuela. Caracas: Sociedad Conservacionista Audubon de Venezuela; 1998.

39. Brooks T, Helgen KM. Biodiversity: a standard for species. Nature. 2010;467: 540–541. pmid:20882005

40. Gutiérrez EE, Helgen KM. Outdated taxonomy blocks conservation. Nature. 2013;495: 314.

41. Duarte JMB, González S, editors. Neotropical Cervidology: Biology and Medicine of Latin American Deer. Jaboticabal: Funep-IUCN. 2010.

42. Barrio J. First records and conservation status of Mazama rufina (Cervidae, Artiodactyla) from Peru. Mastozool. Neotrop. 2010;17(1): 117–122. Available: http://www.scielo.org.ar/scielo.php?script=sci_arttext&pid=S0327-93832010000100012&lng=es&nrm=iso

43. Wisely SM, Maldonado JE, Fleischer RC. A technique for sampling ancient DNA that minimizes damage to museum specimens. Conserv Genet. 2004;5: 105–107.

44. Pucheran J. Note sur une espèce nouvelle de Cerf (Cervus rufinus Bourc. et Puch). Rev et Mag Zool. 1851;2(III): 561–565.

45. R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing. Vienna, Austria. 2014. Available: http://www.R-project.org

46. Dunn OJ. Multiple comparisons among means. J Amer Stat Asso. 1961;56(293): 52–64.

47. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Series B Stat Methodol. 1995;57: 289–300.

48. Revelle W. psych: Procedures for Personality and Psychological Research, Northwestern University, Evanston, Illinois, USA, 2014. Available: http://CRAN.R-project.org/package=psych Version 1.4.8.

49. Stacklies W, Redestig H, Scholz M, Walther D, Selbig J. pcaMethods—a bioconductor package providing PCA methods for incomplete data. Bioinformatics 2007;23: 1164–1167. pmid:17344241

50. Hammer Ø, Harper DAT, Ryan PD. PAST: Paleontological statistics software package for education and data analysis. Palaeontologia Electronica 2001;4(1): 9pp. Available: http://palaeo-electronica.org/2001_1/past/issue1_01.htm

51. Katoh K, Standley DM. MAFFT Multiple Sequence Alignment Software Version 7: Improvements in Performance and Usability. Mol Biol Evol. 2013;30: 772–780. pmid:23329690

52. Lanfear R, Calcott B, Ho SYW, Guindon S. PartitionFinder: combined selection of partitioning schemes and substitution models for phylogenetic analyses. Mol Biol Evol. 2012;29: 1695–1701. pmid:22319168

53. Ronquist F, Teslenko M, van der Mark P, Ayres DL, Darling A, Höhna S, et al. MrBayes 3.2: efficient Bayesian phylogenetic inference and model choice across a large model space. Syst Biol. 2012;61: 539–542. pmid:22357727

54. Zwickl DJ. Genetic algorithm approaches for the phylogenetic analysis of large biological sequence datasets under the maximum likelihood criterion. Ph.D. dissertation, the University of Texas at Austin. 2006.

55. Felsenstein J. Confidence limits on phylogenies: an approach using the bootstrap. Evolution 1985;39: 783–791.

56. Ferguson JWH. On the use of genetic divergence for identifying species. Biol J Linn Soc. 2002;75: 509–516.

57. Baker RJ, Bradley RD. Speciation in mammals and the genetic species concept. J Mammal. 2006;87: 643–662. pmid:19890476

58. Dávalos LM, Russell AL. Sex-biased dispersal produces high error rates in mitochondrial distance-based and tree-based species delimitation. J Mamm. 2014;95: 781–791.

59. Gutiérrez EE, Jansa SA, Voss RS. Molecular systematics of mouse opossums (Didelphidae: Marmosa): assessing species limits using mitochondrial DNA sequences, with comments on phylogenetic relationships and biogeography. Amer Mus Novitates 2010;3692: 1–22.

60. Tamura K, Peterson D, Peterson N, Stecher G, Nei M, Kumar S. MEGA5: Molecular evolutionary genetics analysis using maximum likelihood, evolutionary distance, and maximum parsimony methods. Mol Biol Evol. 2011;28: 2731–2739. pmid:21546353

61. Phillips SJ, Anderson RP, Schapire RE. Maximum entropy modeling of species geographic distributions. Ecol Modell. 2006;190: 231–259.

62. Boria RA, Olson LE, Goodman SM, Anderson RP. Spatial filtering to reduce sampling bias can improve the performance of ecological niche models. Ecol Modell. 2014;275: 73–77.

63. Aiello-Lammens MA, Boria RA, Radosavljevic A, Vilela B, Anderson RP. spThin: an R package for spatial thinning of species occurrence records for use in ecological niche models. Ecography 2015;

64. Hijmans RJ, Cameron SE, Parra JL, Jones PG, Jarvis A. Very high resolution interpolated climate surfaces for global land areas. Int J Climatol. 2005;25: 1965–1978.

65. Anderson RP, Raza A. The effect of the extent of the study region on GIS models of species geographic distributions and estimates of niche evolution: preliminary tests with montane rodents (genus Nephelomys) in Venezuela. J Biogeogr. 2010;37: 1378–1393.

66. Anderson RP, Gonzalez I. Species-specific tuning increases robustness to sampling bias in models of species distributions: An implementation with Maxent. Ecol Modell. 2011;222: 2796–2811.

67. Gutiérrez EE, Boria RA, Anderson RP. Can biotic interactions cause allopatry? Niche models, competition, and distributions of South American mouse opossums. Ecography Ecography 2014;37: 741–753.

68. Gent PR, Danabasoglu G, Donner LJ, Holland MM, Hunke EC, Jayne SR, et al. The Community Climate System Model Version 4. J Clim. 2011;24: 4973–4991.

69. Watanabe S, Hajima T, Sudo K, Nagashima T, Takemura T, Okajima H, et al. MIROC-ESM 2010: model description and basic results of CMIP5-20c3m experiments. Geosci Model Dev. 2011;4: 845–872.

70. Warren DL, Seifert SN. Ecological niche modeling in Maxent: the importance of model complexity and the performance of model selection criteria. Ecol Appl. 2011;21: 335–342. pmid:21563566

71. Radosavljevic A, Anderson RP. Making better Maxent models of species distributions: complexity, overfitting, and evaluation. J Biogeogr. 2014;41: 629–643.

72. Warren DL, Wright AN, Seifert SN, Shaffer HB. Incorporating model complexity and spatial sampling bias into ecological niche models of climate change risks faced by 90 California vertebrate species of concern. Divers Distrib. 2014;20: 334–343.

73. Muscarella R, Galante PJ, Soley-Guardia M, Boria RA, Kass JM, Uriarte M, et al. ENMeval: An R package for conducting spatially independent evaluations and estimating optimal model complexity for Maxent ecological niche models. Methods Ecol Evol. 2014;

74. Soley-Guardia M, Radosavljevic A, Rivera JL, Anderson RP. The effect of spatially marginal localities in modelling species niches and distributions. J Biogeogr. 2014;41: 1390–1401.

75. Hasegawa M, Kishino H, Yano T. Dating of the human-ape splitting by a molecular clock of mitochondrial DNA. J Mol Evol. 1985;22: 160–174. pmid:3934395

76. Fernández‐Mazuecos M, Vargas P. Congruence between distribution modelling and phylogeographical analyses reveals Quaternary survival of a toadflax species (Linaria elegans) in oceanic climate areas. New Phytol. 2013;198: 1274–1289. pmid:23496320

77. Gutiérrez EE, Molinari J. Morphometrics and taxonomy of bats of the genus Pteronotus (subgenus Phyllodia) in Venezuela. J Mammal. 2008;89(2): 292–305.

78. Anderson RP, Gutiérrez EE. Taxonomy, distribution, and natural history of the genus Heteromys (Rodentia: Heteromyidae) in central and eastern Venezuela, with the description of a new species from the Cordillera de la Costa. In Voss R. S. and Carleton M. D. (editors), Systematic mammalogy: contributions in honor of Guy G. Musser. Bull Amer Mus Nat Hist. 2009;331: 33–93.

79. Helgen KM, Lim NTL, Helgen LE. The hog-badger is not an edentate: systematics and evolution of Arctonyx (Mammalia: Mustelidae). Zool J Linn Soc. 2008;154: 353–385.

80. Helgen KM, Pinto CM, Kays R, Helgen LE, Tsuchiya MTN, Quinn A, et al. Taxonomic revision of the olingos (Bassaricyon), with description of a new species, the Olinguito. ZooKeys. 2013;324:1–83. pmid:24003317

81. Bornholdt R, Helgen KM, Koepfli KP, de Oliveira LR, Lucherini M, Eizirik E. Taxonomic revision of the genus Galictis (Carnivora: Mustelidae): species delimitation, morphological diagnosis, and refined mapping of geographic distribution. Zool J Linn Soc. 2013;167: 449–472.

82. Mousseau T, Sikes DS. Almost but not quite a subspecies: a case of genetic but not morphological diagnosability in Nicrophorus (Coleoptera: Silphidae). Biol J Linn Soc. 2011;102: 311–333.

83. James J, Ramakrishnan U, Datta A. Molecular evidence for the occurrence of the leaf deer Muntiacus putaoensis in Arunachal Pradesh, north-east India. Conserv Genet. 2008;9: 927–931.

84. Polziehn RO, Strobeck C. A phylogenetic comparison of red deer and wapiti using mitochondrial DNA. Mol Phylogenet Evol. 2002;22: 342–356. pmid:11884159

85. Anderson RP, Soriano PJ. The occurrence and biogeographic significance of the southern spiny pocket mouse Heteromys australis in Venezuela. Z. Säugetierkunde. 1999;64: 121–125.

86. Gutiérrez EE, Soriano PJ, Rossi RV, Murillo J, Ochoa-G J, Aguilera M. Occurrence of Marmosa waterhousei in the Venezuelan Andes, with comments on its biogeographic significance. Mammalia. 2011;75: 381–386.

87. Williams F, Genoways HH, Braun K. Taxonomy and Systematics In: Genoway H. H. and Brown J. H., editors. Biology of the Heteromyidae. Special Publ, American Soc Mammal. 1993;10: 38–196.

88. Rossi RV. Revisão taxonômica de Marmosa Gray, 1821 (Didelphimorphia, Didelphidae). Ph.D. dissertation, Universidade de São Paulo. 2005.

89. Anderson RP. Preliminary review of the systematics and biogeography of the spiny pocket mice (Heteromys) of Colombia. Revista de la Academia Colombiana de Ciencias Exactas, Físicas y Naturales, 23 (Suplemento Especial). 1999 [2000]; 613–630.

90. Clark PU, Dyke AS, Shakun JD, Carlson AE, Clark J, Wohlfarth B, et al. The Last Glacial Maximum. Science. 2009; 325: 710–714. pmid:19661421

91. Prado JR, Percequillo AR. Geographic distribution of the genera of the tribe Oryzomyini (Rodentia: Cricetidae: Sigmodontinae) in South America: patterns of distribution and diversity. Arq Zool (São Paulo). 2013;44: 1–120.

92. Severinghaus JP, Brook EJ. Abrupt climate change at the end of the last glacial period inferred from trapped air in polar ice. Science. 1999;286: 930–934. pmid:10542141

93. Augustin L, Barbante C, Barnes PRF, Barnola JM, Bigler M, Castellano E, et al. Eight glacial cycles from an Antarctic ice core. Nature. 2004;429: 623–628. pmid:15190344

94. Thomas O. Some notes on the Echimyinae. Annals Ann Mag Nat Hist. (series 8). 1916;18: 294–301.

95. Hershkovitz P. The South American gracile mouse opossums, genus Gracilinanus Gardner and Creighton, 1989 Marmosidae, Marsupialia): a taxonomic review with notes on general morphology and relationships. Fieldiana Zool. 1992;39: 1–56.

96. Rodríguez-Posada ME. Primer registro del ratón de agua del Táchira, Neusticomys mussoi (Rodentia, Cricetidae) en Colombia. Mastozool Neotrop. 2014;21(2): 367–372.

97. Anderson RP, Gutiérrez EE, Ochoa-G J, García FJ, Aguilera M. Faunal nestedness and species–area relationship for small non-volant mammals in “sky islands” of northern Venezuela. Stud Neotrop Fauna Environ. 2012;47: 157–170.

98. Woodman N. A new species of small-eared shrew from Colombia and Venezuela (Mammalia: Soricomorpha: Soricidae: Genus Cryptotis). Proc Biol Soc Wash. 2002;115: 249–272.

99. Ochoa-G J, Aguilera M, Pacheco V, Soriano PJ. A new species of Aepeomys Thomas, 1898 (Rodentia: Muridae) from the Andes of Venezuela. Mamm Biol. 2001;66: 228–237.

100. Thomas O. On seven new small mammals from Ecuador and Venezuela. Ann Mag Nat Hist. (series 7). 1898;1: 451–457.

101. Tate GHH. The taxonomic history of South and Central American oryzomyine genera of rodents (excluding Oryzomys): Nesoryzomys, Zygodontomys, Chilomys, Delomys, Phaenomys, Rhagomys, Rhipidomys, Nyctomys, Oecomys, Thomsomys, Inomys, Aepeomys, Neacomys and Scolomys. Amer Mus Novitates. 1932;581: 1–28.

102. Aguilera M, Pérez-Zapata A, Ochoa-G J, Pascual S. Karyology of Aepeomys and Thomasomys (Rodentia: Muridae) from the Venezuelan Andes. J Mammal. 2000;81: 52–58.

103. Aguilera M, Pérez-Zapata A, Martino A. Cytogenetics and karyosystematics of Oryzomys albigularis (Rodentia, Cricetidae) from Venezuela. Cytogenet Cell Genet. 1995;69: 44–49. pmid:7835085

104. Márquez EJ, Aguilera-M M, Corti M. Morphometric and chromosomal variation in populations of Oryzomys albigularis (Muridae: Sigmodontinae) from Venezuela: multivariate aspects. Z. Säugetierkunde. 2000;65: 84–99.

105. Thomas O. New Neotropical mammals, with a note on the species of Reithrodon. Ann Mag Nat Hist. (series 7). 1901;8: 246–255.

106. Handley CO Jr. Mammals of the Smithsonian Venezuelan Project. Brigham Young Univ Sci Bull, biol ser. 1976;20(5): 1–91.

107. Bisbal FJ. Distribution and habitat association of the carnivores in Venezuela. In: Redford K, Eisenberg J, editors. Advances in neotropical mammalogy. Gainesville: Sandhill Crane Press; 1989. pp. 339–362.