INTRODUCCION

El primer paso para la conservación de cualquier organismo es poder situarlo en un taxón discreto. Una vez un grupo de organismos es reconocido como una especie, otro paso fundamental es determinar en cuantas unidades genéticamente diferentes se divide, determinar la estructura en y entre esas unidades, y la reconstrucción de la historia evolutiva de las mismas. Eso es fundamental para la conservación de la biodiversidad y, también, para una correcta definición de las políticas de conservación de muchas especies (^{Haig 1998}). Sin embargo, en muchas ocasiones, existen discrepancias entre los especialistas en la determinación de qué es una especie, o cual es el estatus de un taxón determinado, debido al concepto de especie que manejen, al igual que por el tipo de datos que se ha utilizado (genética, morfología, etología, ecología, paleontología, etc.) Ruiz-García et al. (²⁰¹²; ²⁰¹⁴; ^2016a;^b). Por ejemplo, (^{Wilting et al. 2015}) mostraron que los planes diseñados para la conservación del tigre (Panthera tigris) han sido fallidos por la falta de consenso en el número propuesto de subespecies.

Los felinos de la América Neotropical no son una excepción. Actualmente, la comunidad científica identifica 11 especies de felinos en Latino América. Ocho de ellos, forman un grupo monofilético denominado el linaje del ocelote (^{Wurster-Hill & Centerwall 1982}; ^{Collier & O’Brien 1985}; ^{Pecon Slattery et al. 1994}; ^{Johnson et al. 1996}; 2006; Johnson & O’Brien 1997; ^{Slattery & O’Brien 1998}). Esas ocho especies son el ocelote (Leopardus pardalis), el margay (Leopardus wiedii), el gato andino (Leopardus jacobita), el gato de las pampas, o colocolo, (Leopardus colocolo) (sin embargo, ^{Garcia-Perea 1994} morfológicamente distinguió tres especies en el seno de este taxón), la güiña o huiña (Leopardus guigna), el gato de Geoffroy (Leopardus geoffroyi), el tigrillo (Leopardus tigrinus), y el tigrillo del sur de Brasil (Leopardus guttulus). Esta última especie ha sido diferenciada recientemente del L. tigrinus por medio de análisis moleculares exhaustivos (^{Trigo et al. 2008}; ²⁰¹³; ²⁰¹⁴). En los últimos dos años, otras dos especies del género Leopardus han sido propuestas. ^{Nascimento & Feijó (2017)} mostraron que la población de tigrillos del noreste de Brasil podría ser considerada una nueva especie utilizando caracteres morfológicos, Leopardus emiliae. Igualmente, Ruiz-García et al. (^2018b; ^2019a), utilizando datos moleculares y morfológicos, detectaron la presencia de una presunta nueva especie de felino en el sur de Colombia que ha sido denominada Leopardus narinensis. Sin embargo, todavía no se ha dado una ratificación total por parte de los especialistas de esas dos últimas presuntas nuevas especies del género Leopardus.

Una de esas especies de este género es el margay (L. wieddi), cuya distribución se da desde México hasta el sur de Brasil, Uruguay y norte de Argentina (Figura 1). Es una especie fuertemente asociada a bosques primarios (^{Bisbal 1989}) y parece ser poco tolerante a disturbios en esos bosques, aunque ha sido reportado fuera de áreas boscosas (^{Vaughan 1983}; ^{Tello 1986}) e, incluso, ha sido encontrada en plantaciones de cocoa y de café en Venezuela (^{Mondolfi 1986}). Se encuentra en raras ocasiones por encima de los 1.500 metros sobre el nivel del mar (^{de Oliveira 1994}). Esta especie está categorizada como casi amenazada por IUCN y en el apéndice I de CITES (^{de Oliveira et al. 2015}).

Morfológicamente, se ha considerado la existencia de 10 subespecies de margay (^{Cabrera 1957}; ^{Hall 1981}; ^{de Oliveira 1998}): (1) L. w. amazonica (^{Cabrera 1917}) con localidad tipo en Tabatinga, Amazonas, Brasil; (2) L. w. boliviae (^{Pocock 1941}) con localidad tipo en Buena Vista, Santa Cruz, Bolivia; (3) L. w. glaucula (^{Thomas 1903}) con localidad tipo en Beltrán, Jalisco, México; (4) L. w. nicaraguae (^{Allen 1919}) con localidad tipo en Chinandega, Nicaragua; L. w. oaxacensis (^{Nelson & Goldman 1931}) con localidad tipo en Cerro San Felipe, Oaxaca, México; L. w. pirrensis (^{Goldman 1914}) con localidad tipo en Cana, Darién, Panamá; (7) L. w. salvinia (^{Pocock 1941}) con localidad tipo en Vera Paz, Guatemala; (8) L. w. vigens (^{Thomas 1904}) con localidad tipo en Igarapé-Assu, Pará, Brasil; (9) L. w. wiedii (^{Schinz 1821}) con localidad tipo en Morro de Azará, rio Mucurí, Bahía, Brasil; y (10) L. w. yucatanica (^{Nelson & Goldman 1931}) con localidad tipo en Mérida, Yucatán, México.

Recientemente, ^{Kitchener et al. (2017)} reconocieron tres subespecies de margay, haciéndolas coincidir con los tres grupos moleculares encontrados por ^{Eizirik et al. (1998)}. L. w. wiedii correspondiente al grupo del noreste de Sudamérica, L. w. vigens correspondiente al grupo del sur de Sudamérica, y L. w. glauculus correspondiente al grupo centro americano. Sin embargo, ^{Nascimento (2010)}, basado en caracteres morfológicos, no encontró variaciones geográficas significativas en el margay.

Los estudios genéticos con esta especie son escasos y de alcance limitado. ^{Eizirik et al. (1998)} analizaron 24 especímenes de L. wieddi procedentes de cautiverio en diversas instituciones de México, Guatemala, Nicaragua, Costa Rica, Panamá, Guyana Francesa, Brasil, Bolivia, y Paraguay. Se analizó un fragmento de la región control del ADN mitocondrial (HVS-I) estimándose niveles de diversidad genética elevadísimos (diversidad haplotípica, H_d=0.985; diversidad nucleotídica, π=0.183). De hecho, la última estimación ha sido una de las más elevadas encontradas para cualquier especie de vertebrado. También detectaron tres agrupaciones geográficas diferentes. Una en Centro América (con dos subgrupos: México-Guatemala; Nicaragua-Costa Rica), otra en el noreste de Sudamérica (norte del Amazonas brasileño y Guyana Francesa), y la tercera en el sur de Sudamérica (este y sur de Brasil, Bolivia, y Paraguay).

Concluyeron que el Río Amazonas y el estrecho del Darién fueron los causantes de la fragmentación de esta especie en esas tres unidades genéticas. ^{Ruiz-García (2001)} analizó cinco marcadores microsatélites en 14 especímenes de esta especie. Doce de ellos procedieron de la Amazonía colombiana (presuntamente de la subespecie L. w. amazonica) y dos procedieron del Departamento de Santa Cruz en Bolivia (presuntamente de la subespecie L. w. boliviae). Los niveles de diversidad genética encontrados para esos marcadores nucleares fueron también muy elevados (H=0.846), al igual que se había determinado para el ADN mitocondrial. ^{Grisolia et al. (2007)} analizaron 25 especímenes de L. wieddi procedentes de diferentes zoológicos de Brasil mediante cuatro marcadores microsatélites. Estimaron un número promedio de alelos por microsatélite muy elevado (11), además de un Contenido de Información Polimórfica (PIC=0.85) y una heterocigosis esperada (0.85), también muy elevados coincidiendo con lo encontrado por ^{Ruiz-García (2001)}.

Sin embargo, hasta la fecha, no se ha realizado un estudio filogeográfico profundo de esta especie con una muestra lo suficientemente amplia para determinar cuál es el número de unidades genéticas diferenciadas y, presuntamente, el número de subespecies que pudiera contener esta especie. Para ello, analizamos 118 especímenes de L. wieddi procedentes de 12 países diferentes de Latino América (México, Guatemala, Belice, Nicaragua, Costa Rica, Panamá, Colombia, Guyana Francesa, Ecuador, Perú, Bolivia, y Brasil) mediante tres genes mitocondriales (mtATP8, mt16S rRNA, y mtND5; 859 pares de bases, pb). Los genes mitocondriales son una herramienta de gran utilidad para tareas filogenéticas y filogeográficas debido a la rápida acumulación de mutaciones, tiempos de coalescencia rápidos, inexistencia de intrones, un alto número de copias por célula (lo cual permite obtener ADN mitocondrial con mayor facilidad en muestras degradadas, o de poca calidad, que ADN nuclear), e inexistencia de recombinación genética debido a una herencia haploide materna (^{Avise 1994}). A pesar de que todos los genes mitocondriales pueden considerarse un único locus ligado, las presiones selectivas y las tasas de evolución pueden ser muy heterogéneas a lo largo de todo el ADN mitocondrial (^{Galtier et al. 2006}; ^{Nabholz et al. 2008}). Debido a todo ello, los genes mitocondriales son más precisos que la mayor parte de los genes nucleares para la reconstrucción de la historia evolutiva de especies muy relacionadas o de subespecies en el interior de especies dadas (^{Moore 1995}).

En el presente estudio se utilizaron fragmentos de tres genes mitocondriales (mtATP8, mt16S rRNA, and mtND5). Estos genes han sido elegidos porque tienen elevada variabilidad genética en felinos debido a su rápida evolución, son de fácil alineación, y son los que han sido secuenciados más frecuentemente en felinos neotropicales, lo cual permite potenciales com paraciones entre diferentes poblaciones o especies de estos felinos (^{Johnson et al. 1996}; ¹⁹⁹⁸; ¹⁹⁹⁶; ^{Johnson & O’Brien 1997}; ^{Slattery & O’Brien 1998}; ^{Culver et al. 2000}). Aunque es objeto de un estudio mostrado en otro sitio, donde se compara la estructura genética del ocelote y del margay (^{Ruiz-García et al. 2019b}), una fracción de los 118 margays analizados fueron secuenciados para sus mitogenomas completos (55 ejemplares; aproximadamente 16.300 pb) y mostraron resultados perfectamente compatibles con los obtenidos con los tres genes mitocondriales citados. Preferimos no utilizar el fragmento de la región control (HVS-I) empleada por ^{Eizirik et al. (1998)} ya que está repleto de inserciones y deleciones que dificulta enormemente la alineación correcta de las respectivas secuencias, lo cual puede ser la causa de la diversidad genética excepcionalmente elevada para el margay encontrada por esos autores.

Los objetivos principales del presente trabajo son los siguientes: 1) Estimar los niveles de diversidad genética para marcadores mitocondriales en el margay; 2) Determinar cuántos grupos genéticamente diferenciados pueden detectarse en el margay y en qué áreas geográficas se encuentran estos; 3) Determinar la posible existencia de estructura espacial a lo largo de la distribución geográfica de esta especie, y 4) Determinar la historia demográfica de la misma a lo largo de su trayectoria evolutiva.

MATERIALES Y METODOS

Muestras colectadas

Se analizaron muestras de 118 margays (pelos, dientes, huesos, y trozos de pieles) procedentes de México (tres muestras representando putativamente L. w. glaucula), Guatemala (ocho muestras representando putativamente L. w. salvinia), Belice (una muestra representando putativamente L. w. salvinia), Nicaragua (una muestra representando putativamente L. w. nicaraguae), Costa Rica (una muestra representando la subespecie L. w. pirrensis), Panamá (dos muestras representando putativamente L. w. pirrensis), Colombia (32 muestras representando putativamente L. w. pirrensis y L. w. amazonica), Guyana Francesa (tres muestras representando putativamente L. w. vigens), Ecuador (dos muestras representando putativamente L. w. amazonica), Perú (32 muestras representando putativamente L. w. amazonica), Bolivia (26 muestras representando putativamente L. w. boliviae) y Brasil (siete muestras representando la subespecie L. w. amazonica) (Fig. 1 y Tabla 1). La totalidad de las muestras representa siete de las 10 subespecies reconocidas morfológicamente. Estas muestras proceden básicamente de animales cazados por comunidades indígenas o colonas a través del Neotrópico (trozos de piel, dientes, y huesos). Algunos ejemplares estaban vivos y se conservaban como mascotas en esas comunidades (pelos con bulbo). Para ello, se solicitó permiso a las autoridades de las comunidades para colectar ese material biológico que ya estaba presente en las comunidades.Las comunidades solo se visitaron una vez, siendo las donaciones voluntarias y no se indujo, en ningún instante,la donación de muestras mediante ninguna transacción económica, o de ningún tipo. Durante el proceso de muestreo (tanto de tejidos de animales previamente cazados,como de ejemplares vivos presentes en las comunidades),se entrevistó a los cazadores de esas comunidades que mencionaron que los especímenes procedían, como mucho,de áreas alejadas entre 5-15 km de las mismas. Las muestrasse obtuvieron en un periodo de 10 años (2006-2016).

Extracción, amplificación y secuenciación de ADN

La extracción de ADN total se llevó a cabo mediante la resina Cheelex 100 TM (^{Walsh et al. 1991}), con la que se procesaron las muestras de pelo con bulbo y de diente (pulpa), y mediante el método de Fenol-Cloroformo (^{Sambrook et al. 1989}) para las muestras de piel. Las muestras de hueso se procesaron con DNeasy Blood and Tissue Kit (Qiagen, Hilden, Germany). Para la amplificación por medio de la Reacción en Cadena de la Polimerasa (PCR), se utilizaron los cebadores descritos por ^{Johnson et al. (1998)} para fragmentos de los genes ATP8 (ATP8.1F: GCATTAACCTTTTAAGTTAAAG y ATP8.2R: GGCGAATAGATTTTCGTTCA, longitud del amplicón: 180 pb) y ND5 (ND5-1F: GTGCAACTCCAAATAAAAG Y ND5-2R: GGGTCTGAGTTTATATATC; longitud del amplicón: 315 pb) y los cebadores descritos por ^{Hoelzel & Green (1992)} (1992) para 16S rRNA (16S.1F: GTGCAAAGGTAGCATAATCA y 16S.4R: TGTCCTGATCCAACATCGAG, longitud del amplicón: 364 pb).

Las reacciones de amplificación se realizaron en un volumen final de 25 µl, conteniendo MgCl2 1 mM (CorpoGen), dNTPs 0.2 mM (BioLabs), 0.1 µM de cada cebador y una unidad de Taq Polimerasa (CorpoGen) con 20-200 ng de ADN molde (dependiendo del método de extracción). Las temperaturas de la PCR se llevaron a cabo a 95°C por 10 minutos seguidos de 35 ciclos de 35 segundos a 95°C, 30 segundos a 52-59°C y 45 segundos a 72°C, con una extensión final de 10 minutos a 72°C. Los productos de la amplificación se corrieron en un gel de agarosa al 2% teñido con bromuro de etidio y fueron visualizados con rayos ultravioleta en un transiluminador. Ambas hebras fueron secuenciadas usando BigDye Terminator v3.1 (Applied Biosystems, Inc.), cuyos productos fueron analizados en un ABI 3730 (Applied Biosystems, Inc.). Las secuencias fueron ensambladas y editadas con el programa Sequencher 4.7 (Gene Codes, Corp., Ann Arbor, MI).

Los números de acceso de GenBank de los 118 ejemplares de margay analizados son los siguientes: de MN005818 a MN005935.

Tabla 1 Localización geográfica y coordenadas de los 118 margays (Leopardus wiedii) secuenciados para tres genes mitocondriales (mtATP8, mt16S rRNA, y mtND5), clasificados por subespecies putativas y localidad geográfica. 

Fig. 1 Mapa con la distribución geográfica del margay (Leopardus wieddi), países muestreados con las subespecies morfológicas putativas muestreadas y tamaños muestrales respectivos. 

Diversidad genética

Diversos estadísticos se utilizaron para determinar la diversidad genética de la muestra total estudiada de margays: sitios polimóricos (S), número de haplotipos (NH), la diversidad haplotípica (H_d), y la diversidad nucleotídica (π). Esos estadísticos fueron calculados con el programa DNAsp 5.1 (^{Librado & Rozas 2009}).

Análisis filogenéticos

Las secuencias para los genes mtATP-8, mt16S rRNA y mtND5 fueron concatenadas, ya que mostraron tener la misma señal filogenética cada una de ellas. Para ello se utilizó la prueba ILD (“Incongruence-length difference test”; ^{Farris et al. 1995}) mediante el programa mILD (^{Planet & Sarkar 2005}). Mediante los programas MrModeltest v2.3 (^{Nylander 2004}) y MEGA 6.05 (^{Tamura et al. 2013}) se determinó el mejor modelo de substitución nucleotídica. Para ello, se utilizó el criterio de información de Akaike (AIC) (^{Akaike 1974}).

Se obtuvieron dos tipos de árboles filogenéticos diferentes. El primero fue un árbol Bayesiano (AB) para determinar las relaciones de los 118 especímenes de margay estudiados y poder detectar posibles grupos monofiléticos que pudieran corresponderse con las subespecies putativas de esta especie. Para ello, se utilizó el programa BEAST v. 1.8.1 (^{Drummond et al. 2012}). Se corrieron cuatro iteraciones independientes usando tres particiones de datos (codón 1, codón 2, codón 3) con seis Cadenas Markov Monte Carlo (MCMC) muestreadas cada 10.000 generaciones por 40 millones de generaciones después de un periodo de “burn-in” de 4 millones de generaciones. La convergencia se verificó utilizando el programa Tracer v1.6 (^{Rambaut et al. 2013}). Se contrastó la verosimilitud para cada generación y se estimó el tamaño efectivo de la muestra (ESS>200) de todos los parámetros a través de cuatro análisis independientes para determinar la convergencia y resultados óptimos. Los resultados de las diferentes corridas se combinaron mediante el programa LogCombiner v1.8.0 (^{Rambaut & Drummond 2013}) y mediante el programa TreeAnnotator v1.8.0 (^{Rambaut & Drummond 2013}). Se utilizó un modelo de especiación de Yule y un reloj molecular relajado con una tasa de distribución log-normal no correlacionada (^{Drummond et al. 2006}). Los valores de probabilidad posterior proporcionaron una evaluación del grado de apoyo de cada nodo en el árbol. Los árboles fueron visualizados ± Importar tabla mediante el programa FigTree v. 1.4 (^{Rambaut 2012}). A partir del AB se obtuvieron los Tiempos de Divergencia del Ancestro Común más Reciente (TMRCA) para los diferentes nodos con probabilidades posteriores elevadas. Debido a que la información sobre el registro fósil para el género Leopardus es escasa, el tiempo prior para la raíz del árbol que se utilizó fue de 3.05 1.1 millones de años (MA) (intervalo de confidencia del 97,5%: 0.89 – 5.21 MA). Este valor se estimó obteniendo el promedio de las estimaciones propuestas para la diversificación del género Leopardus por ^{Johnson et al. (2006)} (valor medio=2.91 MA, valor mínimo= 2.02 MA, valor máximo=4.25 MA) y ^{Li et al. (2016)} (valor medio=3.14 MA, valor mínimo= 1.64 MA, valor máximo=5.03 MA).

Para intentar ratificar los posibles clados detectados en el AB, se obtuvo un árbol de Máxima Verosimilitud (AMV). Para ello se utilizó el programa RAxML v.7.2.6 (^{Stamatakis 2006}). Para seleccionar el mejor modelo de ajuste se utilizaron 50 iteraciones independientes usando tres particiones de datos (codón 1, codón 2, codón 3). Para cada análisis se utilizó el modelo GTR + G + I (^{Tavaré 1986}) para la obtención del AMV y el soporte topológico se estimó mediante 500 réplicas de bootstrap (^{Stamatakis 2006}).

Adicionalmente, para determinar las relaciones entre los haplotipos analizados se utilizó el algoritmo Median Joining Network (MJN) (^{Bandelt et al. 1999}) utilizando el programa Network 4.6.0.1 (Fluxus Technology Ltd). Una ventaja de MJN comparado con los árboles filogenéticos tradicionales es que permite explícitamente la coexistencia de haplotipos ancestrales y descendientes, mientras que los árboles tradicionales tratan todas las secuencias como taxones terminales (^{Posada & Crandall 2001}). Esto permite observar cuáles haplotipos fueron los primeros en originarse y cuáles son los haplotipos derivados más recientes.

Análisis genético espacial

Se realizaron dos tipos de análisis espaciales. El primero de ellos fue una prueba de Mantel (^{Mantel 1967}). Se comparó la matriz de distancias genéticas de 2 parámetros de Kimura (^{Kimura 1980}) entre individuos con la matriz de distancias geográficas entre los mismos. Ambas distancias fueron logaritmizadas y la significación estadística se estimó con 10.000 permutaciones. Para este análisis se utilizó el programa Alleles In Space (AIS) (^{Miller 2005}).

El segundo fue un análisis de autocorrelación espacial mediante el estadístico Ay (^{Miller 2005}). Ay se puede interpretar como la distancia genética promedio entre pares de individuos que caen dentro de una clase de distancia (CD) especificada. Ay toma un valor de 0 cuando todos los individuos dentro de una CD son genéticamente idénticos y toma un valor de 1 cuando todos los individuos dentro de una CD son completamente diferentes. La probabilidad para cada CD se obtuvo utilizando 10.000 permutaciones. Para llevar a cabo este análisis, se definieron cinco CD: 1CD (0-968 km), 2 CD (968-1.937 km), 3 DC (1.937-2.905 km), 4 DC (2.905-3.874 km), y 5 CD (3.874-5.243 km). Los especímenes fueron conectados con la red de ^{Gabriel & Sokal (1969)} y con otras variedades de redes (^{Ruiz-García 1993}; ¹⁹⁹⁴; ¹⁹⁹⁷; ¹⁹⁹⁹; ^{Ruiz-García & Alvarez 2000}).

Cambios demográficos

Para detectar posibles cambios en la historia demográfica del margay se utilizaron tres procedimientos diferentes.

El primero fue la distribución de las diferencias de nucleótidos entre haplotipos/individuos o “mismatch distri bution” (^{Rogers & Harpending 1992}; ^{Rogers et al. 1996}) mediante el programa DNAsp v5.1 (^{Librado & Rozas 2009}). Para ello, se generó la curva esperada de acuerdo a un modelo demográfico de población constante. Se utilizó el estadístico rg (^{Harpending et al. 1993}; ^{Harpending 1994}) para determinar la similitud entre las curvas teóricas y observadas generadas en este análisis. Para aceptar, o rechazar, el modelo de demográfico de población constante se utilizó 10,000 simulaciones de coalescencia. Este análisis permite estimar temporalmente cuando comenzó el cambio demográfico detectado a partir de la expresión τ = 2 µ t, dónde τ es la tasa de mutación por millón de años y t el número de generaciones desde que ocurrió el cambio demográfico. ^{de Oliveira (1998)} describió que la madurez sexual en el margay se alcanza entre los 1 y 2 años de vida. Por ello, se utilizó la estima de 1 o 2 años por generación.

El segundo procedimiento consistió con varios estadísticos, siendo éstos las pruebas F_S de ^{Fu (1997)}, D de ^{Tajima (1989)} y R2 de ^{Ramos-Onsins & Rozas (2002)}. Se obtuvo un intervalo de confianza del 95% y las probabilidades se obtuvieron con 10.000 permutaciones de coalescencia. El programa DNAsp 5.1 (^{Librado & Rozas 2009}) fue utilizado para estimar esos estadísticos.

El tercer procedimiento fue un “Bayesian Skyline Plot” (BSP) mediante los programas BEAST v. 1.8.1. y TRACER v1.6. La opción “Coalescence-Bayesian skyline” fue seleccionada con cuatro pasos con 40 000 000 de generaciones (los primeros cuatro millones se descartaron mediante “burn-in”), con una función kappa con log Normal (1, 1.25) y el tamaño poblacional “Skyline” siguiendo una distribución uniforme (0, infinito; valor inicial 80). Con Tracer v1.5 (^{Rambaut et al. 2013}), se analizaron las densidades marginales de las divisiones temporales y se seleccionó la opción de reconstrucción “Bayesian Skyline”. Se utilizó una variante constante (de un paso) en el “Bayesian Skyline” con el tiempo máximo en el gráfico como la densidad posterior del 95%. La calibración del AB se utilizó para darle un valor al rootHeight. Por ello se analizó el posible cambio demográfico del margay en los últimos dos millones de años.

RESULTADOS

Diversidad genética

Para la muestra de 118 especímenes analizados, se encontró S=240, NH=85, H_d=0.976 0.009, y π=0.035 0.0032. El margay mostró tener una elevada diversidad genética mitocondrial, tal como se había encontrado en los escasos trabajos publicados anteriormente con diferentes tipos de marcadores, tanto mitocondriales como nucleares (^{Eizirik et al. 1998}; ^{Ruiz-García 2001}). Por lo tanto, no existe discrepancia con el uso de los tres marcadores mitocondriales empleados en este estudio con los otros marcadores utilizados por otros autores.

Análisis filogenético

El AB (Fig. 2) mostró que las secuencias encontrados para el margay se diferencian y son monofiléticas con respecto a los de otra especie del género Leopardus, en este caso el gato andino (L. jacobita) (Probabilidad Posterior, PP=1). El inicio de la diversificación de los haplotipos del margay fue estimada en 2.73 MA, coincidiendo con la fase final del Plioceno. En el AB se detectaron 14 clados principales con PP superiores a 0.6.

El clado 1 (PP=1) se diversificó hace unos 0.85 MA y contuvo especímenes de cuatro subespecies putativas morfológicas [L. w. boliviae (1 espécimen de Bolivia), L. w. amazonica (3 especímenes de Colombia), L. w. pirrensis (4 especímenes de la costa norte y pacífica de Colombia), y L. w. glaucula (1 espécimen de México)]. El clado 2 (PP=0.99) comenzó su diversificación alrededor de 1.43 MA. Estuvo integrado tan solo por tres especímenes representando dos de las subespecies morfológicas [L. w. amazonica (2 especímenes de la Amazonía peruana y colombiana) y L. w. nicaraguae (1 espécimen de Nicaragua)]. El clado 3 (PP=0.97) empezó su diversificación alrededor de 0.55 MA y solo estuvo integrado por dos especímenes de L. w. amazonica, ambos de los Llanos Orientales colombianos. El clado 4 (PP=0.92) empezó su diversificación hace 1.21 MA e incluyó a L. w. vigens (1 espécimen de Guyana Francesa), L. w. boliviae (1 espécimen de Bolivia), y, especialmente, L. w. amazonica (12 especímenes procedentes de la Amazonía peruana, colombiana y brasileña, y también de los Llanos Orientales colombianos). El clado 5 (PP=1) es el que presentó una diversificación inicial más antigua estimada en torno a 1.9 MA incluyendo a L. w. boliviae (4 especímenes de Bolivia), L. w. glaucula (1 espécimen de México), L. w. amazonica (2 especímenes de la Amazonía peruana) y L. w. pirrensis (2 especímenes del Pacífico colombiano). El clado 6 (PP=0.71) empezó a diferenciarse alrededor de 0.34 MA conformado únicamente por cinco especímenes de L. w. amazonica procedentes de la Amazonía occidental (Perú, Colombia, y Brasil). El clado 7 (PP=0.86), que comenzó su diversificación alrededor de 0.34 MA, está constituida exclusivamente por cuatro especímenes bolivianos (L. w. boliviae). El clado 8 (PP=0.93), que se originó alrededor de 0.53 MA, también estuvo exclusivamente conformado por siete especímenes bolivianos (L. w. boliviae). El clado 9 (PP=0.99), siendo el segundo más antiguo en originarse hace 1.43 MA, junto con el segundo clado, estuvo compuesto por cinco especímenes de L. w. amazonica procedentes de la Amazonía peruana y brasileña. El clado 10 (PP=0.8) se empezó a diversificar hace 0.62 MA, y estuvo compuesto únicamente por tres especímenes de dos presuntas subespecies morfológicas, L. w. pirrensis (1 espécimen del noroeste de Colombia) y L. w. amazonica (2 especímenes del centro de Perú). El clado 11 (PP=0.92) se diversificó hace 0.27 MA, y estuvo compuesto únicamente por dos especímenes de L. w. vigens procedentes de Guyana Francesa. El clado 12 (PP=0.67) empezó su diferenciación hace 0.87 MA, y estuvo conformado por 12 especímenes de L. w. amazonica, procedentes de la Amazonía peruana, colombiana y ecuatoriana, además de los Llanos Orientales colombianos, y dos especímenes de L. w. pirrensis del Pacífico colombiano. El clado 13 (PP=0.96), con una diferenciación inicial hace 0.59 MA, es de origen centroamericano y está constituido por dos sub-clados también con elevado soporte estadístico. El primero de ellos (PP=0.78), con una diversificación a partir de 0.34 MA, estuvo conformado por los nueve especímenes estudiados de L. w. salvinia procedentes de Guatemala y Belice. El segundo de ellos (PP=0.77), que comenzó su conformación hace 0.47 MA, estuvo compuesto por los únicos tres especímenes centroamericanos muestreados de L. w. pirrensis (dos procedentes de Panamá y uno de Costa Rica). Entre esos dos sub-clados se encontró también un espécimen de L. w. glaucula procedente de México. Finalmente, el clado 14 (PP=0.88), el de más reciente conformación (0.12 MA), estuvo compuesto únicamente por dos especímenes, uno de L. w. amazonica (Brasil) y otro de L. w. pirrensis (norte de Colombia).

Los 14 clados significativos no coincidieron en general con las subespecies morfológicas putativas y todos los procesos de diferenciación se dieron durante el Pleistoceno. El AMV (Fig. 3) ratificó 13 de los 14 grupos de margays encontrados en el AB con porcentajes de bootstraps elevados. El único grupo que no se conservó fue el clado 14 encontrado en el análisis anterior. Sin embargo, las posibles relaciones entre los 13 clados encontrados en ambos árboles filogenéticos fueron diferentes en función de cada uno de los árboles, lo cual pone de manifiesto que la estructura genética del margay es débil.

El procedimiento MJN (Fig. 4) mostró que los haplotipos de las supuestas subespecies putativas de margay se encuentran básicamente entremezclados. Por lo tanto, las conclusiones filogenéticas obtenidas a partir de los análisis con especímenes o con haplotipos fueron similares. El haplotipo 7 es compartido por especímenes clasificados a priori en las subespecies L. w. amazonica, L. w. pirrensis, y L. w. boliviae. El haplotipo 23 es compartido por especímenes de L. w. amazonica y L. w. pirrensis. Desparramados por esta red de haplotipos se encuentra algunas congregaciones de haplotipos que se corresponden con supuestas subespecies putativas. Este es el caso de varias agrupaciones polifiléticas de L. w. amazonica y de L. w. boliviae. La única agrupación monofilética de haplotipos que se correspondería a una supuesta subespecie morfológica es la que se estableció para L. w. salvinia.

Análisis genético espacial

La prueba de Mantel mostró una correlación no significativa entre las distancias genéticas y las distancias geográficas (r=0.073; p=0.06), ya que las distancias geográficas solo explicaron el 0.54 % de las distancias genéticas.

El análisis de autocorrelación espacial con 5 CD también mostró la inexistencia de estructura espacial tanto para el correlograma global (p=0.689) como para las CD individuales (Fig. 5). Así pues, en una distancia de más de 5.000 km, no se encontró evidencia de estructura espacial para el margay.

Cambios demográficos

El análisis de “Mismatch distribution” (Fig. 6) mostró una curva unimodal significativamente diferente a la distribución que se esperaría si la población hubiera tenido un tamaño constante (rg=0.0028, p=0.0074). Este resultado se ajusta bien a una expansión poblacional para la muestra global analizada para el margay. El valor de τ fue de 5.394. Aceptando que una generación en el margay fuera de un año, el inicio de esta expansión poblacional se habría dado hace 405 000 años. Si asumimos que una generación en el margay es de dos años, la misma se habría generado hace unos 202 000 años.

Los estadísticos utilizados para determinar presuntos cambios demográficos durante la historia evolutiva del margay también mostraron, todos ellos, una expansión poblacional significativa (Tajima D=- 2.215, p=0.0013; Fu FS=-48.485, p=0.000001; R2=0.038, p=0.0051).

Por último, el análisis BSP detectó un inicio de expansión poblacional hace 350 000-300 000 años que se vio fuertemente incrementado en los últimos 150 000 años (Fig. 7). Por lo tanto, los tres métodos empleados detectaron una importante expansión poblacional para el margay y, tanto el primer método como el tercero, coinciden en estimar esa expansión poblacional en los últimos 400 000-150 000 años, durante la última fase del Pleistoceno.

Fig. 2 Árbol Bayesiano (AB) con los 118 especímenes de margay (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16S rRNA, mtND5). En los nodos, probabilidades posteriores (PP) y en paréntesis el tiempo de divergencia de los clados respectivos en millones de años. Se señalan los 14 clados significativos encontrados en este análisis. 

DISCUSION

Diversidad genética

Como ya se ha mencionado, la muestra global de margay analizada para tres genes mitocondriales presentó una elevada diversidad genética (π=0.035). La misma es comparable a la de otros carnívoros neo-tropicales que también presentaron elevadas estimas de diversidad genética con los mismos, o similares, genes mitocondriales. Por ejemplo, estos son los casos de la taira (Eira barbara: π=0.042, ^{Ruiz-García et al. 2017}), del jaguar (Panthera onca: π=0.034, ^{Ruiz- García et al. 2013c}), y de la especie hermana del margay, el ocelote (Leopardus pardalis: π=0.031, ^{Ruiz-García et al. 2019b}). Sin embargo, la estima del margay es considerablemente superior a las encontradas en tres especies de zorros sudamericanos (Cerdocyon thous: π=0.019, ^{Tchaicka et al. 2007}; Lycalopex culpaeus: π=0.008, ^{Ruiz-García et al. 2013b}; Lycalopex sechurae: π=0.015, ^{Ruiz-García et al. 2013b}), en tres especies de nutrias (Lontra felina: π=0.005, ^{Vianna et al. 2010}; Lontra longicaudis: π=0.011, ^{Ruiz-García et al. 2018a}; Pteronura brasiliensis: π=0.008, ^{Ruiz- García et al. 2018a}), o en el puma (Puma concolor: π=0.003, ^{Culver et al. 2000}). Incluso, la diversidad genética estimada para el margay con los tres genes empleados (ATP8, 16S rRNA, y ND5) fue muy superior a la encontrada con la región control (D-loop) para especies de felinos neotropicales que tienen distribuciones geográficas más restringidas que el margay y que están en caso de mayor vulnerabilidad en términos de conservación. Estos son los casos de la guigna (Leopardus guigna: π=0.004, Napolitano et al. 2013) y del gato andino (Leopardus jacobita: π=0.004, ^{Cossíos et al. 2009}; ^{Ruiz-García et al. 2013a}). Por lo tanto, la diversidad genética encontrada con los tres marcadores genéticos empleados es lo suficientemente elevada para encontrar estructura genética y espacial en el seno de esta especie en caso que existiera. La ausencia de la misma no es atribuible al tipo de marcadores mitocondriales utilizados.

Fig. 3 Árbol de Máxima Verosimilitud (AMV) con los 118 especímenes de margay (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16S rRNA, mtND5). En los nodos, los porcentajes de bootstraps. Se señalan los 13 clados significativos encontrados en este análisis. 

Fig. 4 Red de haplotipos (“Median Joining Network”) con los haplotipos detectados en este análisis con 118 especímenes de margay (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16SrRNA, mtND5). Círculos rojos (mv) indican haplotipos no encontrados, o desaparecidos, que conectan los haplotipos detectados. 

Fig. 5 Correlograma obtenido con el estadístico de autocorrelación espacial Ay con cinco clases de distancias para los 118 especímenes de margay (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16SrRNA, mtND5). El correlograma no fue significativo (p = 0.689) indicando ausencia de estructura espacial para esta especie. 

Evolución temporal, espacial, y demográfica del margay

Hasta la fecha presente, solo se ha publicado un estudio relacionado con la filogeografía del margay (^{Eizirik et al. 1998}). Sin embargo, ese estudio se realizó únicamente con 24 especímenes procedentes básicamente de animales en cautiverio y dejando sin muestrear grandes zonas geográficas donde habita el margay. El presente estudio intentó cubrir esas carencias ya que el tamaño muestral se incrementó hasta 118 especímenes y se cubrieron muchas zonas geográficas que no fueron abarcadas por ^{Eizirik et al. (1998)}.

El inicio de la diversificación mitocondrial del margay hace 2.7 MA se correlaciona bien con loscambios climáticos que se originaron por el cierre del Istmo de Panamá durante el final del Plioceno3.1-2.8 MA; (^{Marshall et al. 1979}; ¹⁹⁸²; ^{Marshall 1985}; ¹⁹⁸⁸; ^{Webb 1985}; ^{Coates & Obando 1996}; ^{Webb 1997}). Existe discusión entre diversos autores en determinar cuando apareció el Pleistoceno. La mayor parte de los trabajos parecen apuntar el inicio del Pleistoceno alrededor de 1.8 MA (extinción de los discoastéridos, 1.8-1.6 MA, ^{Berggren 1987}; paleomagnetimo, 2.0-1.8 MA, ^{McDougall 1979}; Potasio/Argón, 1.86 MYA, ^{Jordá 1995}). Sin embargo, recientemente, otras fuentes reivindican el origen del Pleistoceno hace 2.59 MYA, incorporando el periodo Gelasiano del Plioceno al Pleistoceno (^{International Stratigraphic Chart 2008}). Por lo tanto, el primer indicio de diversificación mitocondrial para el margay (2.7 MA) se habría dado en la última fase del Plioceno.

Fig. 6 Análisis de “Mismatch distribution” mostrando una expansión poblacional significativa para la muestra de 118 especímenes de margays (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16SrRNA, mtND5). 

Fig. 7 Análisis “Bayesian Skyline Plot” mostrando una expansión poblacional significativa desde hace 300.000 años para la muestra de 118 especímenes de margays (Leopardus wieddi) analizados para tres genes mitocondriales (mtATP8, mt16S rRNA, mtND5). 

En el AB se observó que los clados 5, 2, 9 y 4 (solo se discuten los 13 clados que fueron comunes entre el AB y el AMV) fueron los que se diversificaron más antiguamente (1.9, 1.43, 1.43 y 1.21 MA, respectiva- mente) durante la primera fase del Pleistoceno. El primero incluyó un número limitado de especímenes (nueve) pero que ocuparon la mayor parte del rango geográfico analizado (desde México a Bolivia). Su diferenciación se dio dentro del periodo Gelasiano (transición Plio-Pleistocénica). El segundo es también un clado muy limitado (únicamente tres especímenes) pero que engloba especímenes amazónicos y centroamericanos (de Nicaragua). El tercero fue un clado exclusivamente conformado por ejemplares amazónicos (Perú y Brasil). El cuarto fue un clado que engloba un número mayor de ejemplares, básicamente cubriendo una buena parte de la Amazonía occidental y de los Llanos Orientales colombianos, pero, también con un ejemplar boliviano y otro de la Guyana Francesa, siendo exclusivamente un clado sudamericano.

Estos tres últimos clados (diferenciación 1.4-1.2 MA) parecen asociarse al periodo pre-Pastoniano (1.30-0.8 MA), el cual constituyó la etapa más fría y seca de la primera gran glaciación detectada en Norte América y en Europa Central (Nebraska-Günz). Fue un periodo extremadamente seco con importante fragmentación de bosques. En Sudamérica se relaciona con el Ensenadense tal como se comentará unas líneas más adelante. En este periodo se detectó también fragmentación genética en otras especies de Carnivora, tales como en el gato de pajonal (Leopardus colocolo; ^{Cossíos et al. 2009}), en ciertas formas de tigrillos (Leopardus tigrinus; ^{Ruiz-García et al. 2018b}), en el jaguarundí (^{Ruiz-García et al. 2018c}), en zorros del género Lycalopex (^{Ruiz-García et al. 2013b}) o en la taira (Eira barbara; ^{Ruiz-García et al. 2017}).

Otros clados significativos fueron, no obstante, de diversificación temporal mucho más reciente siendo el caso de los clados 13b, 7, 6, 13a, y 11 (0.47, 0.36, 0.34, 0.34, 0.27 MA, respectivamente). Este periodo coincide con el Bonaerense (^{Ameghino 1889}) y se extendió desde hace 0.5 a 0.13 MA. Se caracterizó por la biozona del perezoso gigante (Megatherium americanum), por el incremento de mamíferos holárticos en Sudamérica (por ejemplo, roedores y cérvidos), y por la existencia de suelos con una flora muy diversa. Fue un periodo de una cierta estabilidad climática con un ambiente húmedo (extensión de bosques) y cálido lo cual se correlaciona perfectamente con la expansión poblacional detectada para el margay con los tres procedimientos empleados. La aparición de algunos clados, que se diferenciaron en esa franja temporal, pudo haberse dado en un periodo más frío y seco anterior al Bonaerense, el Ensenadense (2.0 a 0.5 MA). Al inicio del Ensenadense se caracterizó por la biozona de Tolypeutes pampaeus (un armadillo). Al comienzo de este periodo el clima fue substancialmente cálido y la fauna de la provincia de Buenos Aires era típicamente tropical o chaqueña (Tapirus, Catagonus, Cyonasua, Clyomys). Sin embargo, desde hace un millón de años hasta hace 0.5 MA, el clima se enfrió y se secó considerablemente. El registro fósil de mamíferos de la provincia de Buenos Aires se caracterizó entonces por taxones de zonas esteparias secas, como la fauna actual de la Patagonia (Lama, Lestodelphys, Lyncodon) (^{Forasiepi et al. 2007}). Esto pudo producir una fragmentación de los bosques y de las antiguas poblaciones de margay en Latino América. Posteriormente algunas de esas poblaciones se expandieron de forma notable durante el Bonaerense y, a nivel genético, fueron detectadas sus expansiones durante ese periodo. No se puede descartar, no obstante, que la fragmentación de los grupos más modernos de margay se diera en algunos periodos, relativamente, breves de sequía y enfriamiento que se ocurrieron dentro del contexto de clima más cálido, húmedo y estable del Bonaerense. En Europa, la glaciación Riss (hace 0.2 MA) estaría intercalada temporalmente durante el Bonaerense sudamericano. Por ejemplo, ^{Hooghiemstra (1984)} detectó, tan solo para el Pleistoceno de Colombia, hasta 27 posibles cambios climáticos importantes con una periodicidad de 100 000 A. Sin embargo, solo hay clara evidencia de la última gran glaciación que comenzó hace 116 000 A (^{van der Hammen & González 1963}; ^{van der Hammen et al. 1981}; ^{Helmens 1988}; ^{Flórez 1992}). Sin embargo, hay una cierta evidencia de que se pudieron registrar glaciaciones en el noroeste de Sudamérica anteriores a los 128 000 años. ^{Herd (1982)} detectó una morrena de más de 100 000 años en el Departamento del Tolima (Colombia), lo cual significa que era anterior a la última gran glaciación que se inició hace 116 000 años. También, ^{van der Hammen et al. (1981)} (Sierra del Cocuy, Departamento de Boyacá, Colombia), ^{Helmens (1988)} (Sumapaz, Departamento de Cundinamarca, Colombia), o ^{Flórez (1992)} (Departamento de Caldas, Colombia) encontraron evidencias de glaciaciones anteriores a unos 200 000 años. Esto significa que la fuerte expansión detectada para el margay durante el Bonaerense pudo darse en poblaciones que se fragmentaron anteriormente durante el Ensenadense o, también, que se fragmentaron durante el propio periodo Bonaerense, en ciertos momentos más secos y que, posteriormente, habrían generado múltiples contactos secundarios.

Todo lo comentado pone de manifiesto que han existido múltiples procesos de diversificación en el interior del margay, que esos procesos de diversificación han ocurrido en diferentes épocas, y muchos de ellos se han sobrepuesto geográficamente unos a otros debido a la considerable capacidad de expansión poblacional de esta especie. Especialmente, las zonas colonizadas originalmente (distribución de clados más antiguos) fueron posteriormente colonizadas, especialmente, en los últimos 400.000 A, por grupos de margays originados más modernamente. Por ello, encontramos, en unas mismas zonas geográficas, varios clados (AB y AMV) viviendo en sim- patría. Algunos ejemplos son suficientes. Para el AB, en la Amazonía occidental (especialmente peruana y colombiana) detectamos hasta siete clados diferentes (clados 1, 2, 4, 5, 6, 9, y 12) coexistiendo en simpatría. En México, únicamente muestreamos tres ejemplares y cada uno de ellos pertenece a un clado diferente. Uno, en el clado 5 más antiguo y conjuntamente con especímenes sudamericanos (1.9 MA), otro en el clado 1, también con todos los otros especímenes sudamericanos (0.85 MA), y el tercero en el clado 13 (0.59 MA) que es exclusivamente centroamericano. En Bolivia (incluyendo las zonas amazónicas y centro-andinas del país), convergen cinco clados diferentes.

Tres de ellos (clados 1, 4 y 5) se encuentran ampliamente extendidos por la distribución geográfica del margay, mientras que otros dos son exclusivamente de origen boliviano. Los tres ejemplares muestreados en Guyana Francesa también ejemplarizan este concepto. Dos de los especímenes analizados conforman un clado (11) únicamente detectado en esa zona geográfica, pero el tercer espécimen pertenece al clado 4 ampliamente expandido por la cuenca del Amazonas. Estos resultados ponen de manifiesto que ni la Cordillera de los Andes, ni el tapón del Darién fueron obstáculos a múltiples conexiones entre las poblaciones de margay. Todo ello muestra una sobreimposición de grupos colonizadores que llegaron a las respectivas áreas geográficas en diferentes épocas temporales del Pleistoceno debido al elevado potencial de flujo génico para esta especie. Este fenómeno de sobreimposición de haplogrupos generados en diferentes épocas pero que actualmente viven en simpatría también ha sido detectado en otras especies de mamíferos neotropicales (Cebus albifrons, ^{Ruiz-García et al. 2010}); Tapirus terrestris, (^{Ruiz-García et al. 2016c}); Eira barbara, (^{Ruiz-García et al. 2017}); Puma yagouaroundi, (^{Ruiz-García et al. 2018c}).

Esta sobreimposición geográfica de diferentes haplogrupos es el fenómeno que explica la no detección de estructura espacial significativa, tanto para el test de Mantel como para la autocorrelación espacial. Utilizando los mismos marcadores mitocondriales, otros estudios detectaron una muy significativa estructura espacial en el tigrillo (Leopardus tigrina; ^{Ruiz-García et al. 2018b}) y en el ocelote (^{Ruiz-García et al. 2019b}). También, para el jaguar (^{Ruiz-García et al. 2019c}) y el jaguarundí (Puma yagouroundi; ^{Ruiz-García et al. 2016c}), la estructura espacial a través del Neotrópico fue significativa aunque no tan marcadamente como en los dos taxones anteriormente comentados. Por lo tanto, la inexistencia de estructura espacial en el margay no es atribuible a las características de los genes mitocondriales empleados porque esos mismos marcadores sí detectaron una estructura espacial significativa en otros felinos neotropicales.

Un análisis mitogenómico, con una fracción de los 118 especímenes de margay estudiados aquí (55 especímenes), mostraron la existencia significativa de los clados 2, 4, 5, 7, 8, 12, y 13b (^{Ruiz-García et al. 2019b}). Los otros clados detectados en el presente estudio no aparecieron porque no se pudo obtener el mitogenoma completo de especímenes representativos de esos clados.

Sistemática del margay

Todo lo comentado, conjuntamente con la poca es- tructura filogeográfica y la ausencia de estructura espacial significativa, entra en conflicto con las 10 subespecies morfológicas putativas y de localización geográfica específica establecida por ^{Cabrera (1957)}, ^{Hall (1981)}y ^{de Oliveira (1998)}, pero, también, entra en conflicto con las tres subespecies sugeridas por ^{Kitchener et al. (2017)} a partir del estudio molecular de ^{Eizirik et al. (1998)}.

Los especímenes supuestamente pertenecientes a la subespecie L. w. amazonica se encontraron en 10 grupos significativamente diferentes. Cuatro de ellos incluyeron especímenes exclusivamente clasificados a priori como L. w. amazonica. Por lo tanto, esos cuatro grupos son monofiléticos pero polifiléticos entre sí. Por lo tanto, en el área geográfica considerada tradicionalmente para L. w. amazonica, encontramos 10 grupos diferentes, cuatro de ellos exclusivos de la Amazonía occidental. El estudio de ^{Eizirik et al. (1998)} no incluyó ningún espécimen de esta amplia zona de Sudamérica por lo que no pudo aportar nada al conocimiento filogeográfico de esta especie en dicha área.

La distribución geográfica tradicional de L. w. pirrensis estaba establecido en el sur de Centroamérica (Costa Rica y Panamá) y costa norte de Colombia y zona este-andina (Pacífica) de Colombia y Ecuador. El presente estudio detectó ejemplares del norte de Colombia y de la zona este-andina (Pacífica) colombiana en cinco grupos diferentes y, en todos ellos, junto a especímenes de otras supuestas subespecies de origen oeste-andino. Si aceptamos que el origen del actual margay se dio en la Amazonía (^{Ruiz-García et al. 2019b}), al menos, en cinco ocasiones esta especie ha pasado de la zona oeste-andina (Llanos y Amazonía) a la este-andina (Pacífica) por lo que la Cordillera de los Andes no pareciera haber sido una barrera geográfica consistente para el margay. Por lo tanto, esta supuesta subespecie (L. w. pirrensis) no parece ser una entidad real en el norte y pacífico colombiano. ^{Eizirik et al. (1998)} tampoco muestrearon especímenes de esa zona y ese estudio no pudo contribuir a resolver la situación sistemática de esa población de margay. Sin embargo, sí detectamos un grupo monofilético y significativo para los especímenes muestreados en Panamá y Costa Rica que se corresponden al área geográfica del holotipo de esta presunta subespecie. Este grupo exclusivamente del sur de Centroamérica sí se correspondería, entonces, con el tradicional L. w. pirrensis.

En territorio del Departamento de Santa Cruz(Bolivia), se definió el holotipo de L. w. boliviae. A priori los especímenes muestreados en Bolivia, los clasificamos en esa subespecie. Sin embargo, encontramos animales bolivianos en cinco grupos significativamente diferentes. Dos de ellos son monofiléticos(todos los especímenes de esos grupos son originariamente bolivianos) pero polifiléticos entre sí. Si el taxón L. w. boliviae fuera real, sin embargo, no podemos determinar a cuál de los dos grupos monofiléticos correspondería el mismo. Adicionalmente,dos de los otros tres grupos donde se encontraron especímenes bolivianos estuvieron integrados básicamente por ejemplares a priori clasificados como L.w. amazonica. Una posible interpretación es que en el territorio boliviano existen dos linajes endémicos de esa zona, pero también se pueden localizar ejemplares que se originaron en otras áreas de la cuenca amazónica (incluyendo al norte del Río Amazonas) y que se sobre-impusieron geográficamente a los dos linajes endémicos de ese país. Por lo tanto, la validez de L. w. boliviae es muy dudosa. Igualmente, este resultado pone en evidencia que el Río Amazonas no fue una barrera geográfica para la dispersión efectiva de esta especie a un lado, y al otro, del mencionado Río. ^{Eizirik et al. (1998)} analizó un solo ejemplar procedente de Bolivia que quedó asociado al grupo del sur de Sudamérica detectado por ellos. ^{Kitchener et al. (2017)} clasificaron ese ejemplar como L. w. wiedii. Nosotros no estudiamos ningún espécimen del rango de distribución de esa presunta subespecie por lo que no podemos hacer inferencias acerca de ella. Sin embargo, resulta evidente que al haber estudiado 26 especímenes bolivianos (versus el único ejemplar de ^{Eizirik et al. 1998}), podemos afirmar que la situación taxonómica de esta especie en Bolivia es mucho más compleja que lo que se determinó con el estudio de un único individuo.

Dos de los tres especímenes analizados procedentes de Guyana Francesa formaron un grupo monofilético que se correspondería bien a la zona geográfica donde se definió el holotipo de L. w. vigens. ^{Eizirik et al. (1998)} analizaron unos cuantos ejemplares de esa zona geográfica y también detectaron un grupo diferenciado en la misma. ^{Kitchener et al. (2017) c}oncluyó que, en esta zona al norte de la desembocadura del Río Amazonas en Brasil y de las Guyanas, está presente L. w. vigens. Nosotros concordamos con esa idea. Sin embargo, no hay que olvidar que en la misma área geográfica también detectamos un ejemplar que hizo parte de un grupo mayoritariamente integrado por especímenes clasificados como L. w. amazonica. Por lo tanto, en esa área geográfica pueden estar viviendo en simpatría especímenes asignables a L. w. vigens con otros procedentes de otras áreas de la cuenca amazónica.

El caso de los taxones del margay en Centroamérica es muy interesante. ^{Kitchener et al. (2017)}, a partir de ^{Eizirik et al. (1998)}, concluyó que únicamente existía una subespecie, L. w. glaucula. Nuestros resultados no confirmaron esa afirmación. Ya comentamos más arriba que se pudo distinguir L. w. pirrensis en Panamá y sur de Costa Rica. Todos los especímenes muestreados en Guatemala y Belice también conformaron un grupo monofilético que se corresponde perfectamente con el origen del holotipo de L. w. salvinia. Sin embargo, el único ejemplar que analizamos de L. w. nicaraguae poseyó un haplotipo asociado a algunos haplotipos de origen sudamericano. El caso de los tres ejemplares mexicanos también es llamativo. Dos de ellos, se asociaron más intensamente a haplotipos sudamericanos que a los haplotipos de pirrensis o salvinia, mientras que el tercero se asoció de manera significativa a los de los dos taxones centroamericanos. Este último ejemplar podría potencialmente representar L. w. glaucula. Sin embargo, más especímenes de esa zona de México deberían estudiarse para poder determinar con exactitud esta aseveración. Al no haber muestreado, en México, especímenes correspondientes a las subespecies putativas L. w. oaxacensis y L. w. yucatanica no podemos hacer inferencias respecto a ellas. Por lo tanto, la situación en Centroamérica parece mucho más compleja que lo especificado por ^{Eizirik et al. (1998)} y ^{Kitchener et al. (2017)} con la existencia de un único grupo o una única subespecie (L. w. glaucula). Sin haber analizado ejemplares de dos presuntas subespecies de la zona, detectamos dos grupos bien diferenciados en Centroamérica (pirrensis y salvinia) y, al menos, cuatro penetraciones desde Sudamérica a Centroamérica. Las tres primeras, más antiguas, fuertemente relacionadas con actuales acervos sudamericanos, y una cuarta más moderna que dio lugar a los haplotipos exclusivamente encontrados en Centroamérica. Estos resultados ponen en evidencia que el tapón del Darién tampoco fue una barrera geográfica importante para la dispersión del margay.

Por lo tanto, el presente estudio molecular únicamente detectó tres grupos que pueden corresponderse a las subespecies morfológicas tradicionalmente consideradas: L. w. salvinia (Guatemala, Belice), L. w. pirrensis (pero restringida únicamente al sur de Costa Rica y Panamá), y L. w. vigens (noreste de Brasil al norte del Río Amazonas y Guyanas), aunque en ese territorio también pueden encontrase individuos y haplotipos característicos de otras zonas de la cuenca Amazónica. Igualmente, este estudio pone de manifiesto que el análisis de un número reducido de muestras y que no cubra ampliamente la distribución geográfica de la especie objeto de estudio puede mostrar una visión simplificada de la historia evolutiva y sistemática de la misma. Esto es lo que ocurre con el estudio de ^{Eizirik et al. (1998)}. Esas circunstancias favorecen encontrar un número más reducido de grupos diferenciables y que estas agrupaciones parezcan monofiléticas cuando no lo son necesariamente. ^{Eizirik et al. (1998)} concluyeron que el Río Amazonas y el tapón del Darién eran barreras geográficas que separaban grupos genéticos diferenciables de margays. Nuestro estudio mostró que, al multiplicar por cinco el tamaño muestral, ambos procesos geográficos no fueron barreras impermeables a los eventos colonizadores del margay. Nuestro análisis también mostró que la Cordillera de los Andes tampoco imposibilitó la colonización de esta especie a un lado y al otro de la misma. En concordancia con nuestros hallazgos, ^{Nascimento (2010)}, desde una perspectiva morfológica, tampoco encontró barreras geográficas que diferenciaran las poblaciones de margay.

Es necesario ampliar el número de especímenes a analizar, especialmente, en el rango geográfico de las tres presuntas subespecies no representadas en este análisis. También es necesario incrementar el análisis de muestras en ciertas áreas de la zona de distribución de L. w. glaucula (México) y de L. w. nicaraguae (norte de Costa Rica y Nicaragua). Igualmente, importante es el análisis de marcadores nucleares que permitan constatar el grado de mezcla entre los individuos que viven en simpatría pero que pertenecen a haplogrupos que han llegado en diferentes épocas a una misma zona geográfica. Si se hubiera dado introgresión genética podríamos encontrar que el ADN mitocondrial de un ejemplar se correspondiera con un agrupamiento diferente a la clasificación morfológica del mismo, mientras que el ADN nuclear podría coincidir con la clasificación morfológica preliminar. Igualmente, también podría ocurrir en el caso de que se diera hibridación que el ADN mitocondrial de ejemplares de una misma zona geográfica presente considerables diferencias por haber evolucionado en regiones geográficas diferentes, mientras que el ADN nuclear de los mismos ejemplares se haya homogeneizado por flujo genético efectivo entre los ancestros de esos especímenes.

Por último, no es descartable la posible hibridación entre el margay y otras especies del género Leopardus, especialmente con el ocelote. La importancia de la hibridación en la historia evolutiva de ciertas especies de ese género (por ejemplo, para L. tigrinus) ha sido constatado por Trigo et al. (²⁰⁰⁸; ²⁰¹³; ²⁰¹⁴) y ^{Ruiz-García et al. (2018b)}. Por lo tanto, es importante analizar comparativamente los resultados obtenidos para el margay con los obtenidos para otras especies del mismo género.