Materiales y Métodos

Consideraciones taxonómicas. Para la asignación de las poblaciones estudiadas a las distintas especies de Abies mesoamericanos seguimos a ^{Aguirre-Planter et al. (2012)}, que a su vez se apoyó en ^{Liu (1971)} y ^{Farjon y Rushforth (1989)}, quienes consideran que en México existen ocho especies de Abies, seis de las cuales son endémicas al país. Siguiendo este criterio, algunos datos de presencia se reasignaron según los parámetros marcados por estos autores. Por ejemplo, apoyados por datos morfológicos y genéticos se consideró a Abies mexicana Martínez como sinónimo de Abies vejarii Martínez, por lo que los datos de ocurrencia de ambos taxa fueron fusionados. De igual manera, los datos de Abies colimensis Rushforth & Narave se asignaron a los de A. religiosa por ser una especie sinónima de esta última (^{Farjon, 2013}). En el caso de los demás abetos del occidente de México, la situación es un poco más complicada ya que de acuerdo con ^{Rushfort (1989)} los registros de Abies guatemalensis Rehder var. jaliscana Martínez y de Abies religiosa var. emarginata Loock & Martínez (Martínez, 1948) estarían en sinonimia con A. flinckii (^{Rushforth, 1989}). Esta última especie, aunque no es reconocida por algunas autoridades, ha mostrado una diferencia genética significativa con respecto a sus congéneres (^{Aguirre-Planter et al., 2000}, ²⁰¹²; ^{Jaramillo-Correa et al., 2008}), por lo que es de interés ecológico y evolutivo ver si ésta se refleja al nivel de su nicho ecológico. Sin embargo, en un trabajo reciente ^{Vázquez-García et al. (2014)} exponen evidencia fenológica y morfológica suficiente para proponer que las poblaciones del oeste de Jalisco registradas antiguamente como A. guatemalensis var. jaliscana deben de ser elevadas al nivel de especie como A. jaliscana (Martínez) Mantilla, Shalisko & A.Vázquez, siendo diferentes del resto de las poblaciones asignadas a A. flinckii (^{Rushforth, 1989}). Ante esta evidencia se decidió modelar las poblaciones consideradas por ^{Rushforth (1989)} como A. flinckii (de aquí en adelante, A. flinckii sensu lato) y las poblaciones a A. jaliscana y A. flinckii (sensu Vázquez-García) por separado. Finalmente, se determinaron por separado los nichos de las dos variedades de A. duranguensis Martínez (A. durangensis var. coahuilensis y A. durangensis var. durangensis), que aunque a nivel genético son muy similares (^{Jaramillo-Correa et al., 2008}; ^{Aguirre-Planter et al., 2012}), su distribución disyunta en las Sierras Madre Oriental y Occidental, respectivamente, podría indicar nichos ecológicos muy distintos. Cabe resaltar que una de las especies consideradas válidas por ^{Farjon (2013)}, Abies hidalgensis Debreczy, I.Rácz & Guízar, no se incluyó en el presente estudio por contar sólo con dos localidades conocidas; esto se reflejaría en modelos con un soporte estadístico muy pobre.

Obtención de los datos de presencia. Los registros de presencia de cada una de las especies se obtuvieron de diversas fuentes. La primera de ellas fue la base de datos de las colectas de ^{Eguiarte y Furnier (1997)}, depositados en la CONABIO bajo la clave B138. Posteriormente se consultó el portal Global Biodiversity Information Facility (“^GBIF”<http://www.gbif.org> consultado entre enero y diciembre de 2011), que incluye datos de herbarios de todo el mundo, cuyos registros fueron complementados con los del Herbario Nacional (MEXU) del Instituto de Biología de la Universidad Nacional Autónoma de México (UNAM). Algunos de los datos antiguos de esta última base de datos tuvieron que ser geo-referenciados conforme a los estándares delineados por ^{Wieczorek et al. (2004)} y al manual de procedimientos de geo-referenciación de localidades de CONABIO (^{Nuñez-Merchand, 2008}) antes de ser incluidos en el presente estudio. Los datos obtenidos fueron depurados en varias fases antes de iniciar las modelaciones del nicho ecológico de cada taxón. En la primera de ellas se eliminaron los datos repetidos y aquellos cuya identidad taxonómica estaba en duda. Posteriormente, se excluyeron los datos mal georeferenciados y/o que representaban valores atípicos geográficos; es decir aquellos puntos en donde se reportó una especie de abeto, pero que se encuentran fuera de su distribución natural, como en costas o desiertos. Luego se excluyeron los valores atípicos climáticos; es decir, las localidades de colecta que se encuentran en zonas en las que es difícil determinar si constituyen o no parte de la distribución natural de la especie. Para ubicar estos valores atípicos se utilizaron dos métodos implementados en DIVA-GIS ver. 7.5 (^{Hijmans et al., 2012}): (1) el método de Reverse jackknife (^{Chapman, 2005}), que ha sido recomendado para conjuntos de datos con una distribución normal de valores (^{Scheldeman y van Zonneveld, 2011}), como aquellos de muchas observaciones para cada taxón (n > 20) y (2) el método de 1.5x Interquartile range (1.5 IQR; ^{Scheldeman y van Zonneveld, 2011}), que se recomienda para conjuntos de datos con poca cantidad de observaciones por taxón (por ejemplo, n < 20). Para todos los análisis se consideró que la cantidad mínima de variables en las que un punto de presencia debía tener valores extremos para ser considerado un punto atípico climático debía ser de tres (^{Scheldeman y van Zonneveld, 2011}). Por último, para evitar problemas de sobre-modelación, los puntos de presencia se depuraron para dejar un sólo registro por celda (cada celda tenía 30 segundos de arco), lo cual se realizó utilizando la biblioteca Raster (^{Hijmans y van Etten, 2012}) para el leguaje de programación R (^{R Core Team, 2014}).

Obtención de datos ambientales y determinación del área accesible (M). Se obtuvieron las 19 variables bioclimáticas de Worldclim (^{Hijmans et al., 2005}), así como las capas de orientación, pendiente e índice topográfico de humedad (Wetness Index) de la base de datos geográfica Hydro1k (^{HYDRO1k Elevation Derivative Database}) en formato Esri grid a una resolución espacial de 30 segundos de arco para el área correspondiente a México, Estados Unidos y Centroamérica. Estas capas fueron recortadas en áreas más pequeñas para ser utilizadas en los análisis de cada especie. El criterio para el recorte fue reducirlas al área accesible para cada especie o área M; es decir, aquellas áreas en donde cada especie está o se supone que podría estar dado el conocimiento biológico que se tiene de ellas, de sus capacidades de dispersión, y por no existir grandes barreras ni discontinuidades ambientales que pudieran limitar su establecimiento. El concepto de área accesible para una especie hace referencia al diagrama BAM propuesto por ^{Soberón y Peterson (2005)}, este esquematiza que el área de distribución de esta especie en determinado tiempo es el resultado de una combinación de variables ambientales o scenopoeticas (A), de variables bióticas (B) y del área que la especie puede explorar (M).

Para ayudar a determinar el área M de cada especie se utilizaron las ecorregiones terrestres del mundo reportadas en ^{Olson et al. (2004)}. Se escogieron aquellas que coincidían con la ubicación de los sitios de colecta para cada especie, evitando ecorregiones que aunque similares, no tuvieran registros. Esto se hizo siguiendo la lógica de que aunque la dispersión del polen y las semillas de Abies es por viento, la germinación de las semillas y el establecimiento de las plántulas se favorecen en condiciones más o menos sombreadas y al resguardo de otros miembros de la misma especie, sugiriendo que la posibilidad de colonizar zonas abiertas o exentas de bosque es baja (^{Williams, 2009}). Las ecorregiones escogidas coinciden en su mayoría con zonas montañosas y/o templadas donde prosperan bosques de coníferas apropiados para el establecimiento de los abetos (^{Enright y Hill, 1995}). Solamente en dos casos (A. durangensis var. durangensis y A. vejarii) las poblaciones coincidieron con una sola ecorregión, para los otros, las fronteras entre ecorregiones contiguas se eliminó para tener una única área continua. Las áreas resultantes, que representan el área probablemente accesible (M), se usaron como “clip” para recortar el conjunto de capas bioclimáticas a emplearse en el análisis de cada especie, luego de ser transformadas a formato ASCII. La selección y edición de las ecorregiones, así como el recorte y transformación de las capas bioclimáticas, se realizó utilizando el programa ArcMap10 (^{ESRI, 2011}).

Reducción del conjunto de datos ambientales y extracción de información climática. Es sabido que utilizar una gran cantidad de capas bioclimáticas puede conducir a errores en las predicciones por sobreajuste de los modelos (^{Peterson y Nakazawa, 2008}). Para evitar esto, se redujo el número de variables que mostraran multicolinealidad a partir de un análisis de correlación de Pearson realizado con los valores de cada pixel de cada una de las capas utilizadas, conservando sólo una de las variables altamente correlacionadas (correlación ≥ 0.8). Este análisis se condujo por medio del programa ENMtools ver. 2.1 (^{Warren et al., 2010}). En la medida de lo posible se escogieron aquellas variables que pudieran ser interpretadas más fácilmente, es decir aquellas como la “temperatura media anual” o la “precipitación anual”, que conllevan un menor tratamiento matemático si se les compara con variables como la “oscilación anual de la temperatura”, que es el resultado de la sustracción de las variables “máxima temperatura del mes más cálido” y “mínima temperatura del mes más frío”. Cuando estos criterios no se cumplieron, las variables bioclimáticas a conservar se escogieron de manera aleatoria.

Modelado de nicho ecológico. Los análisis de distribución potencial se realizaron con el algoritmo de máxima entropía del programa Maxent ver. 3.3.3k. Este utiliza una técnica de aprendizaje de máquina para escoger los modelos más consistentes a partir la información disponible (^{Phillips et al., 2006}). Este programa se escogió dado que se ha demostrado con distintas simulaciones que genera buenas predicciones, aun cuando se utilizan muestras pequeñas (< 10; ^{Phillips et al., 2006}; ^{Pearson et al., 2007}). Los datos de presencia finales y las variables climáticas seleccionadas recortadas con el área M correspondiente sirvieron de insumo para generar los modelos de nicho ecológico potencial para cada especie. Los parámetros de elaboración utilizados fueron aquellos que venían por omisión en el programa (^{Phillips y Dudik, 2008}), a excepción de las opciones Extrapolate y Do clamping que se desactivaron para evitar extrapolaciones artificiales en los valores extremos de las variables ecológicas (^{Elith et al., 2011}), dado que nuestras especies están restringidas a zonas montañosas.

Se obtuvo una salida de tipo logística usando un umbral de mínima presencia de entrenamiento (minimum training presence) que incluyó el 75 % de los puntos para el entrenamiento del modelo y 25 % para validar el mismo. Dichas validaciones se realizaron inicialmente con una prueba binomial, en la que se supone que los modelos son mejores que el azar cuando arrojan un valor de P < 0.01. Posteriormente se utilizó el valor AUC (area under the curve) del análisis de ROC (receiver operating characteristic) para determinar la validez de los modelos. Los modelos con valores de AUC entre 0.7-0.9 para los puntos de entrenamiento y prueba se consideraron razonablemente buenos y los modelos con valores por arriba de 0.9 se catalogaron como muy buenos (^{Peterson et al., 2011}). Sin embargo, es de notar que la utilidad de los análisis ROC ha sido cuestionada para los algoritmos que utilizan datos de solo presencia, ya que en teoría este análisis también requiere datos de ausencias verdaderas (y no solo de pseudoausencias), además de que pondera igual los errores de omisión que los de comisión (^{Lobo et al., 2007}; ^{Peterson et al., 2008}). Por tales motivos, los modelos finales se validaron también por medio del análisis ROC parcial, que ha sido diseñado para tratar de subsanar estas deficiencias (^{Peterson et al., 2008}). Este análisis se realizó con en el programa Tool for Partial-ROC (^{Narayani, 2008}), utilizando un 50 % de los puntos de evaluación independientes re-muestreados en 1,000 réplicas bootstrap y fijando un error de omisión no mayor al 5 % (1-omission threshold > 0.95), luego se realizó una prueba z para determinar si los valores de proporciones AUC parciales de los modelos de nicho fueron estadísticamente mejores que un modelo al azar (AUC = 1.0). Los insumos para el análisis, o sea el número total de pixeles similares con sus valores de idoneidad y el valor de idoneidad para el pixel de cada sitio de presencia para cada análisis, así como los histogramas de frecuencias de las proporciones AUC para los valores de las distribuciones nulas arrojadas por Tool for Partial-ROC, fueron obtenidos con el programa R y la librería Raster. Los archivo tipo ASCII generados por el programa Maxent de los modelos de nicho ecológico para cada especie fueron proyectados posteriormente sobre un mapa del área de distribución correspondiente por medio del lenguaje de programación R (^{R Core Team 2014}), utilizando las bibliotecas Raster (^{Hijmans y van Etten, 2012}), Maptools (^{Bivand y Lewin-Koh, 2015})y Rworldmap (^{South, 2013}).

Caracterización ambiental y diferencias entre nichos. Para determinar si existen diferencias entre los nichos de las especies estudiadas, particularmente aquellas que aún son fuente de controversias taxonómicas, se realizó un análisis de componentes principales (ACP) con los datos bioclimáticos y topográficos de todas las localidades para cada una de las especies incluidas en este estudio. Posteriormente los valores (scores) de los tres primeros ejes del ACP fueron utilizados para probar si existían diferencias significativas entre los nichos de las especies por medio de un análisis no paramétrico de varianza multivariado (MANOVA), ya que no cumplían con los requisitos de normalidad y homoscedasticidad. Los valores significativos de este MANOVA indican una separación de las especies en el espacio ambiental, por lo que posteriormente se realizó un análisis de Kruskal-Wallis para determinar en cuales de los tres ejes existía diferenciación entre las especies. Posteriormente se realizó una prueba post hoc de Gao, que es una prueba de comparaciones pareadas no paramétrica (^{Gao et al., 2008}), para identificar los pares de especies estadísticamente diferentes en cada uno de los ejes. La extracción de los valores bioclimáticos y topográficos, así como los análisis estadísticos, se realizaron con la ayuda del lenguaje de programación R (^{R Core Team, 2014}) utilizando las bibliotecas Raster (^{Hijmans y van Etten, 2012}), Vegan (^{Oksanen et al., 2013}) y Nparcomp (^{Konietschke, 2015}).

Nichos ecológicos y Áreas Naturales Protegidas. Con el ánimo de contribuir a la conservación de los Abies a partir de nuestros resultados y tener datos para su posible manejo, se identificó el porcentaje de registros para cada especie y el área de los nichos ecológicos construidos con estos que coinciden con un área natural protegida en México. Para hacerlo primero se calculó el porcentaje de puntos de colecta para cada especie que están dentro de las poligonales de las áreas naturales protegidas federales y estatales (^{Bezaury-Creel et al., 2009}). Después se calcularon para cada especie el número de pixeles por nivel de idoneidad ambiental (bajo, medio y alto) considerando un umbral de mínima presencia de entrenamiento de los modelos de nicho obtenidos, y se sobrepusieron a las poligonales de las áreas naturales protegidas federales y estatales, posteriormente se calculó el porcentaje de pixeles que quedaron dentro de estas áreas con respecto al total de pixeles por nivel de idoneidad ambiental.

La proyección de los nichos ecológicos en México, la sobreposición de las poligonales de las áreas naturales protegidas y el cálculo de los porcentajes de áreas de distribución potencial que están protegidos, totales y por nivel de idoneidad del hábitat, se realizó con ayuda del lenguaje de programación R (^{R Core Team, 2014}) y las bibliotecas Raster (^{Hijmans y van Etten, 2012}), Maptools (^{Bivand y Lewin-Koh, 2015}).

Resultados

El número final de registros conservados para cada especie después del proceso de depuración de las bases de datos, y que fueron utilizados para modelar sus nichos ecológicos, se presentan en la Tabla 1. La mayor cantidad de datos eliminados fueron de sitios repetidos (incluyendo las localidades en un mismo pixel) o mal georreferenciados, siendo las especies con más datos suprimidos Abies concolor (con 406 de 536 registros iniciales), A. religiosa (284 de 398) y A. guatemalensis (200 de 259). De las otras especies se obtuvieron registros más depurados, dada su distribución restringida, por lo que prácticamente no se eliminaron datos. Los análisis de Reverse jackknife y de 1.5x Interquartile range (1.5 IQR) permitieron detectar sólo 29 localidades con valores climáticos atípicos para todas las especies (Tabla 2), siendo A. hickelii Flous & Gaussen la especie que perdió más registros con este criterio (7) y A. durangensis var. coahuilensis ningún registro; las gráficas resultantes de dichos análisis se pueden solicitar a los autores. En promedio, se conservaron 52.25 registros por especie, siendo A. concolor el taxón con más localidades (127) y A. vejarii y A. durangensis var. coahuilensis las de menor representación (con 25 y 14, respectivamente). Estos últimos valores no son sorprendentes dada la rareza de estos taxa (^{Farjon, 1990}; ^{Farjon y Filer, 2013}). La base depurada con la relación de coordenadas y nombres de las localidades, colectores y fechas de colecta empleadas para este estudio puede ser solicitada a los autores o consultada en formato DarwinCore en el portal de CONABIO ligada al proyecto JM015 (<http://www.conabio.gob.mx/institucion/cgi-bin/datos.cgi?Letras=JM&Numero=15>).

Tabla 1 Procedencia de los registros de Abies Mesoamericanos usados en la modelación de su nicho ecológico potencial. 

CONABIO: Comisión Nacional para el Conocimiento y Uso de la Biodiversidad.

IBUNAM-MEXU: Herbario Nacional de México bajo custodia del Instituto de Biología de la UNAM.

GBIF: Global Biodiversity Information Facility.

*En GBIF y en la base del INECOL.

Tabla 2 Relación de puntos climáticos extremos excluidos de los análisis finales. Todos los puntos fueron determinados por el método de Reverse jackknife exceptuando el punto correspondiente a A. vejarii localizado por el método 1.5 IQR. 

Los registros retenidos se distribuyeron en un total de 98 ecorregiones terrestres, siendo Abies concolor la especie que se distribuye en un mayor número de regiones (43), mientras que A. durangensis var. durangensis y A. vejarii únicamente estuvieron distribuidas en una ecorregión cada una. En el Apéndice 1 se puede consultar una relación de las ecorregiones terrestres del mundo utilizadas para obtener el área M de cada especie.

En total se eligieron de cinco (Abies hickelii) a ocho capas (A. concolor y A. guatemalensis) de variables no-colineales para modelar el nicho de las diferentes especies de Abies mexicanos. La relación final de dichas capas se enlista en la Tabla 3. Las capas elegidas en más ocasiones fueron bio19 (precipitación del cuarto más frío), común para seis especies, seguida de elevación, bio2 (intervalo promedio de temperaturas diurnas), bio18 (precipitación del cuarto más cálido) y pendiente, que fueron compartidas por cinco especies. De este conjunto de variables, la altitud tuvo una alta correlación negativa (< -0.9) con alguna o todas las capas relacionadas con temperatura (bio1 a bio11) en todas las especies, exceptuando a A. concolor para la cual las correlaciones negativas más altas fueron solo de -0.696 con bio11 (temperatura promedio del cuarto más frío) y de -0.692 con bio1 (temperatura media anual). Esto indica que la disminución de temperatura con el incremento de la altitud es más marcado en las zonas montañosas de áreas tropicales e intertropicales, que en las zonas de latitudes extratropicales, como en el oeste de Estados Unidos, en donde se distribuye A. concolor. Las variables que se relacionaron con el aspecto del terreno (pendiente y orientación) no presentaron correlación alguna con ninguna de las demás variables y no parecieron ser tan importantes para caracterizar y marcar diferencias entre los nichos de las especies, con la excepción de A. durangensis var. coahuilensis (ver más adelante).

Tabla 3 Relación de las capas bioclimáticas y topográficas utilizadas para cada uno de los análisis finales de nicho ecológico. 

Los modelos de nicho ecológico para todas las especies analizadas (Figura 1) y para Abies flinckii sensu lato (Figura 2) fueron mejores que el azar según la prueba binomial, que en todos los casos mostró valores significativos (P < 0.01). Estos modelos fueron igualmente catalogados como “muy buenos” según los valores de las pruebas ROC estándar, que fueron mayores a 0.900 para todas las especies. Finalmente, las pruebas de ROC parcial arrojaron razones de AUC con valores significativos (P < 0.0001) mayores a 1 (Tabla 4), indicando nuevamente que los modelos obtenidos son estadísticamente mejores que el azar.

Figura 1 Proyección geográfica de los modelos de nicho ecológico para las especies de Abies de México y Centroamérica, incluyendo A. concolor. La escala de colores representa las áreas de idoneidad del hábitat, siendo mayor en los tonos más rojos. 

Figura 2 Proyección geográfica del modelo de nicho ecológico para Abies flinckii s. l. La escala de colores representa las áreas de idoneidad del hábitat, siendo mayor en los tonos más rojos. 

Tabla 4 Valores de las proporciones AUC del análisis de ROC parcial. 

*Desviación estándar.

En general, las especies presentaron zonas de alta idoneidad (> 0.666) discontinuas y restringidas a ciertas zonas montañosas (ver Figura 1). Sin embargo, hubo algunas diferencias entre taxa, pues para Abies concolor las poblaciones están en alturas que promedian los 1,670 m s.n.m., mientras que en las especies mesoamericanas las condiciones templadas adecuadas para los abetos (ver más adelante) sólo se alcanzan en alturas que promedian los 2,480 m s.n.m. No obstante, existen algunas poblaciones por debajo de los 1,000 m s.n.m., incluyendo 23 localidades de A. concolor, nueve de A. guatemalensis y dos de A. duranguensis var. coahuilensis, indicando que algunos sitios relativamente bajos también pueden tener condiciones climáticas adecuadas para los abetos.

El análisis de componentes principales (ACP) para todas las especies mostró que los tres primeros componentes principales explican el 70 % de la varianza (el 90 % solo se alcanza hasta el séptimo componente) siendo cada uno de estos tres componente responsable del 40, 21.5 y 9 % de la varianza, respectivamente. La gráfica de los tres primeros componentes (Figura 3) muestra que todas las especies están más o menos solapadas ambientalmente con excepción de Abies concolor, que es la especie que más se separa de las demás. Los resultados del MANOVA no paramétrico mostraron que al menos una especie es estadísticamente distinta (F = 56. 431, g.l. = 9, P = 0.0009) en su espacio ambiental. Las pruebas de Kruskal-Wallis muestran que hay diferencias estadísticas entre las especies para cada uno de los tres primeros componentes principales (Tabla 5). Asimismo, la prueba de Gao (Apéndice 2) revela que en general la mayoría de las especies analizadas son similares entre sí con excepción de A. concolor, que para el primer componente (CP1) es la única especie que es estadísticamente distinta en sus atributos ambientales con respecto a todas las demás, aunque con los otros componentes muestra similitudes con A. durangensis var. coahulensis (CP2 y CP3) y con A. flinckii (CP2).

Figura 3 Análisis de componentes principales (ACP) de las variables bioclimáticas y topográficas que muestran las diferencias ecológicas entre las especies de abetos Mesoamericanos. Los puntos correspondientes de Abies flinckii s.l. no se marcaron de color ya que se superponen con los de A. flinckii y con los de A. jaliscana. 

Tabla 5 Sumario del análisis de multivariado de varianza (Manova) no paramétrico y del análisis de Kruskal-Wallis entre los tres primeros componentes del ACP. 

Los datos crudos para todas las especies (datos no mostrados) indican que Abies concolor ocurre en las zonas con mayor estacionalidad de temperaturas (bio4) con un promedio para todos los sitos de 9.5 °C. Por otra parte, los demás abetos se distribuyen en zonas con una estacionalidad de temperatura menos marcada y un poco más cálidas, con temperaturas promedio anuales que van de los 12 °C en las localidades de A. religiosa a los 16.4 °C en las poblaciones de A. guatemalensis. La precipitación a lo largo del área de distribución de A. concolor también es muy distinta a la de los demás taxa mexicanos, con la mayor cantidad de humedad concentrada alrededor de los meses más fríos (bio16 y bio19) y con un promedio de 554 mm por año, mientras que para los otras especies la precipitación está concentrada en los meses de verano o repartida a lo largo del año y con promedios anuales que van de 598 mm para A. durangensis var. coahuilensis a 1,322 mm para A. hickelii.

Cuando se sobrepusieron las poligonales de las zonas naturales protegidas estatales y federales a los modelos de nicho ecológico de las especies de abeto presentes en México, se observó que el porcentaje de colectas y de las zonas de mayor idoneidad ambiental (> 0.6) que está dentro de un área natural protegida (ANP) federal o estatal varía mucho con cada especie, pero en ningún caso es del 100 % (Apéndice 3). En general, las especies con una distribución geográfica menor fueron las que mostraron mayor porcentaje de registros y de áreas de alta idoneidad protegidas, como Abies vejarii con un 56 y 76.16 % o Abies durangensis var. coahuilensis con un 78.57 y 63.65 % respectivamente. En el caso de A. concolor, ningún punto de colecta en México cae dentro de una ANP, aunque de la fracción de nicho ecológico potencial de esta especie presente en México (3.1 %) el porcentaje de áreas con alta idoneidad que pudieran estar protegidas dentro de una ANP es de un 63.87 %. En el caso de A. religiosa, una de las especies con mayor distribución en el Faja Volcánica Transmexicana, el porcentaje de puntos de colecta y de áreas con alta idoneidad protegidas es de un 41.67 % y un 60.9 % respectivamente, mientras que para otras especies en la Faja Volcánica Transmexicana como A. flinckii, A. flinckii s. l. y A. hickelii este porcentaje es más bajo y corresponde al 26.09 % y 33.2, 19.51 % y 31.94 y al 7.89 % y 24. 3 %, respectivamente; para A. jaliscana se registran 11.11 % de colectas dentro de una ANP, pero sólo el 1.96 % de su distribución potencial coincidió con una, constituyéndose en la especie menos protegida del país. La siguiente especie menos protegida es A. guatemalensis, con sólo 8.16 % de registros y un 8.9 % de áreas de alta idoneidad ambiental dentro de una ANP. Cabe anotar que paradójicamente, esta especie es la más abundante de Mesoamérica con un 83.7% de registros y un 78.6 % de su distribución potencial ubicada en México. Otra especie que también tiene una baja protección es A. durangensis var. durangensis, que presenta solo el 20.59 % de los registros dentro de una ANP con áreas de alta idoneidad muy fragmentadas y alejadas sobre la Sierra Madre Oriental y de las cuales solo el 17.45 % estaría protegida.

Discusión

Este estudio es el primero en determinar y comparar el nicho ecológico potencial de la mayoría de las especies del género Abies en México. Los modelos de nicho ecológico obtenidos presentaron un desempeño alto de acuerdo a los resultados de las pruebas binomial, de ROC y ROC parcial (^{Peterson et al., 2008}). Esto se vio corroborado en el hecho de que para todas las especies las zonas con mayor idoneidad ambiental fueron disyuntas y estuvieron limitadas a las partes altas de las sierras coincidiendo con los conocimientos básicos de la biología de las mismas (p.ej. ^{Farjon, 2013}). Los análisis multivariados derivados de los resultados de estos modelos mostraron que la especie ecológicamente más divergente fue A. concolor mientras que los taxa restantes comparten características ecológicas similares. Es interesante anotar que A. concolor ya había mostrado una divergencia significativa a nivel genético (^{Aguirre-Planter et al., 2012}), lo que sugiere que para esta especie, la diferenciación evolutiva y ecológica parecen estar correlacionadas (^{Van Valen, 1976}; ver también ^{Aguirre-Planter et al., 2012} para una discusión al respecto). Estudios similares con marcadores nucleares (Aguirre-Planter et al., 2000; ^{Múgica-Gallart, 2013}) y de cloroplasto (^{Jaramillo-Correa et al., 2008}) también habían indicado que A. flinckii s.l. también es una entidad genética independiente de las demás especies de Abies en México, pero los resultados indican que esta especie es poco distinguible ecológicamente del resto de sus congéneres. Esta misma conclusión se puede ampliar a A. jaliscana y A. flinckii s.s., las dos especies que Vázquez-García y colaboradores reconocen como diferentes a partir de datos morfológicos. Estas similitudes podrían sugerir una evolución por conservadurismo de nicho para las especies de este gran clado, aunque esta hipótesis aún deberá ser probada en estudios futuros y tal vez integrada a análisis filogenéticos con un mayor número de marcadores nucleares.

Volviendo a Abies concolor, esta especie se encuentra actualmente en pocas poblaciones en el norte de México, con la mayor parte de su distribución en grandes áreas del oeste y sudoeste de los Estados Unidos, en regiones que se caracterizan por tener inviernos fríos y precipitaciones concentradas durante el trimestre más frío del año. Las demás especies de abetos mexicanos están distribuidas en zonas con inviernos más templados y con las precipitaciones más altas concentradas en el trimestre más cálido del año. La divergencia genética que se ha observado entre A. concolor y los demás Abies mesoamericanos podría entonces estar relacionada con adaptaciones a estos factores ambientales, que ya han sido señalados como motores importantes en la evolución de otras coníferas (p. ej. ^{Eckert et al., 2010}; ^{Mosca et al., 2012}; ^{Prunier et al., 2012}). Sin embargo, esto también deberá ser probado con estudios genómicos de asociación similares a los efectuados en coníferas modelo como Pinus taeda o Picea glauca (^{Eckert et al., 2010}; ^{Prunier et al., 2012}).

Las similitudes entre los nichos ecológicos de las demás especies de abeto de México son probablemente un reflejo de su origen relativamente reciente y la falta de monofilia recíproca (ver ^{Aguirre-Planter et al., 2012}; ^{Xiang et al., 2015}). Ésta se explica por probables ciclos de aislamiento y contacto secundario entre poblaciones y especies después de la migración del género hacia Mesoamérica entre el Mioceno y Plioceno (^{Eguiarte y Furnier, 1997}; ^{Graham, 1999}; ^{Aguirre-Planter et al., 2000}, ²⁰¹²). De esta manera, los abetos mesoamericanos parecen encontrarse en una “zona gris de especiación” (^{de Queiroz, 1998}, ²⁰⁰⁵), en donde, debido al reciente tiempo de separación de las poblaciones, la monofilia recíproca es incompleta y las poblaciones pueden no formar grupos fenotípicamente y ecológicamente distinguibles o diagnosticables, por lo que varios conceptos de especie entran en conflicto y no hay un consenso sobre cuantas especies existen realmente. Los ciclos de aislamiento y contacto arriba mencionados probablemente estuvieron mediados por el “seguimiento ecológico” que las poblaciones de abetos hicieron de las condiciones templadas características de sus nichos. Este seguimiento ecológico parece haber influenciado los patrones genéticos y de distribución de poblaciones de muchas otras coníferas (^{Gugger et al., 2011}; ^{Moreno-Letelier et al., 2013}) y ha sido sustentado con datos palinológicos que indican que los bosques templados de Mesoamérica se extendieron por debajo de los 1,000 m s.n.m. durante el último máximo glacial hace unos 21 ± 2 ka (^{Bush et al., 2009}). De esta forma, durante las fases de expansión en el pasado, los bosques de coníferas aumentaron sus áreas, facilitando el intercambio de genes entre poblaciones y especies, mientras que en las fases más cálidas (como el presente interglacial), las poblaciones se fragmentaron y se desplazaron hacia las zonas más altas y templadas (^{Ramírez-Barahona y Eguiarte-Fruns, 2013}).

El empleo de nuevos algoritmos y aproximaciones estadísticas, como la utilización de kernels de suavización acoplados a análisis de solapamiento de nicho ecológico para evitar errores por sesgos en los muestreos (^{Broennimann et al., 2012}), además de la utilización de capas a mayor resolución y otras generadas a partir de trabajos de percepción remota, podrían mejorar la resolución de las diferencias ecológicas y climáticas en las especies estudiadas. Por lo que se espera que los resultados de dichos análisis puedan ser integrados en futuros estudios comparativos a nivel filogenético-ecológico para los Abies mexicanos. Por el momento, los resultados ecológicos, genéticos y morfométricos disponibles para este género en el país parecen indicar que la taxonomía tradicional basada en caracteres morfológicos ha conducido a confusiones debido a una elevada plasticidad fenotípica en muchos taxa (^{De Kroon et al., 2005}). Por ejemplo, en un estudio previo ya se había mostrado que existe una gran similitud morfométrica entre las especies de Abies aquí estudiadas, hasta el punto de ser virtualmente indistinguibles (^{Strandby et al., 2009}). Desde este punto de vista, parece indispensable incorporar nuevos métodos para afinar las descripciones taxonómicas disponibles, como las simulaciones de nicho ecológico, el uso de marcadores moleculares y el establecimiento de ensayos de procedencias en campo, en donde se pueda estimar plasticidad fenotípica de las especies de forma adecuada.

Los estudios sobre modelado de nicho ecológico de coníferas, y en especial de abetos en México, son pocos. Uno de los primeros estudios sobre este tema fue el de ^{Téllez-Valdés et al. (2004)} en el cual se exploraba el modelado bioclimático como una herramienta de manejo forestal en cuatro especies de pino de importancia comercial. Posteriormente, ^{Sáenz-Romero et al. (2012)} modelaron el nicho ecológico de Abies religiosa en relación con su importancia como refugio de la mariposa monarca y las afectaciones ante distintos escenarios de cambio climático. Los otros estudios existentes se enfocan en utilizar el modelado de nicho ecológico como herramienta auxiliar para explicar algunos patrones filogeográficos (^{Gugger et al., 2011}; ^{Moreno-Letelier et al., 2013}). Es así que la importancia de este trabajo radica en que los datos aquí generados son los primeros en comparar ecológicamente a casi todas las especies de abetos en México y nos ayudarán a poner a prueba distintas hipótesis filogeográficas y al mismo tiempo obtener datos de partida para estudios de conservación y cambio climático (ver más adelante).

Áreas Naturales Protegidas y conservación de abetos en México. Actualmente, la mayoría de las zonas de mayor idoneidad ambiental para los abetos mexicanos están fuera de las áreas naturales protegidas (ANP), con excepción de Abies religiosa y de especies con distribuciones geográficas restringidas, como A. vejarii o A. duranguensis var. coahuilensis. Lo anterior coincide con otros estudios que muestran que en México las ANP resguardan una proporción muy baja de especies (< 20 %), sobre todo las de aquellas en alguna categoría de riesgo (^{Londoño-Murcia y Sánchez-Cordero, 2011}), que sólo tienen un pequeño porcentaje de su distribución potencial dentro de un ANP (^{Contreras-Medina et al., 2010}). Es de resaltar que, dado que los nichos de los Abies mexicanos son muy estrechos y que los niveles de diferenciación genética nuclear entre poblaciones son muy altos (^{Eguiarte y Furnier, 1997}; ^{Aguirre-Planter et al., 2000}), se debe tratar de conservar la mayor parte de éstas como se propone para el género Agave en México (^{Eguiarte et al., 2013}). En este sentido, sería recomendable generar nuevas áreas naturales protegidas que incluyan las zonas de mayor idoneidad ambiental, especialmente para los taxa poco protegidos, como A. guatemalensis y A. durangensis var. durangensis.

Un caso especial es el de Abies flinckii y A. jaliscana, ya que según la NOM059-SEMARNAT-2010 (^{SEMARNAT, 2010}) A. flinckii es una especie en peligro. Sin embargo, la nueva propuesta de reconocer como una especie válida (A. jaliscana) a las poblaciones más occidentales de A. flinckii s.l. (^{Vázquez-García et al., 2014}), deja a A. flinckii en un estado aun más vulnerable de lo que se había presupuestado y A. jaliscana casi sin protección (al menos en un sentido estricto). Este último punto se ve agravado por el hecho de que esta nueva especie es la que tiene una menor área potencial dentro de una ANP, por lo que sería de especial interés ampliar la zonas de protección para esta especie y para A. flinckii (independientemente del estatus taxonómico de cada una). Esta sugerencia se ve reforzada por la clara diferenciación genética que estos dos taxa en su conjunto (A. flinckii s.l.) muestran con el resto de abetos de Mesoamérica. Por otro lado, sería urgente realizar estudios genéticos y ecológicos sobre A. hidalgensis para verificar su estatus taxonómico y, si es el caso, resguardar sus poblaciones. En este sentido serían necesarias nuevas exploraciones en la zona donde se distribuye (Hidalgo, cerca de Metepec) para aumentar el número de localidades disponibles en bases de datos públicas. Debe recordarse que este taxón fue excluido del presente estudio, pues apenas se contó con dos localidades (ver ^{Debreczy y Rácz, 1995}) en dichas bases, lo que es insuficiente para desarrollar modelos de nicho adecuados.

Debido al acelerado calentamiento global (^{Zeng et al., 2004}; ^{Peng et al., 2014}), es muy probable que las zonas de mayor idoneidad ambiental para los Abies comiencen a diferir de las encontradas en este estudio (ver ^{Sáenz-Romero et al., 2012}). Actualmente, ya se han registrado cambios en el tamaño y localización de las zonas de idoneidad ambiental para muchas plantas en todo el mundo, registrándose cambios más severos en zonas montañosas, donde muchas especies tienden a elevar su límite altitudinal dependiendo la especie y la zona geográfica (^{Wilson et al., 2005}; ^{Lenoir et al., 2008}; ^{Parolo y Rossi, 2008}; ^{Telwala et al., 2013}). De esta forma, es urgente realizar análisis con diferentes escenarios de cambio climático para determinar cuáles serán estás nuevas zonas de idoneidad ambiental y determinar cuáles de los poblaciones actuales están en riesgo de desaparecer. Con estos datos se podrían proponer nuevos límites para las áreas naturales protegidas y tratar de implementar programas de rescate para estas especies. A este respecto, ^{Sáenz-Romero et al. (2012)} propusieron para A. religiosa que la migración asistida a zonas de mayor altitud e idoneidad ambiental debe ser una estrategia a considerarse. Para las demás especies de abeto se podrían sugerir estrategias similares; sin embargo, estudios más detallados sobre la diversidad genética, la plasticidad fenotípica y el potencial evolutivo de cada población y especie son necesarios para poder mejorar los programas de conservación existentes (^{Benito-Garzón et al., 2011}).