Área de estudio y trabajo de campo

El área de estudio se encuentra en la provincia fisiográfica Sierra Madre Oriental donde convergen dos subprovincias: Pliegue Saltillo Parras y Gran Sierra Plegada (INEGI, 2009a, b). El estudio incluye siete poblaciones representativas de la distribución natural de P. remota en México, las cuales fueron georreferenciadas con el uso de un GPS (Garmin eTrex 20x Ltd, Olathe, Kansas, EUA) mediante coordenadas UTM NAD 27.

La población de La Gavía (LG) está ubicada en el municipio Castaños, Coahuila, dentro de la subprovincia Pliegue Saltillo Parras (Fig. 1 A, B). Taxonómicamente es caracterizada como la localidad tipo de P. remota (^{Favela, 2004}). En el municipio Santa Catarina, subprovincia Gran Sierra Plegada, entre las montañas que limitan a Monterrey, Nuevo León, con Saltillo, Coahuila (^{Perry, 1991}; ^{Favela, 2004}) se ubica el resto de las poblaciones: Molino de Viento (MV), El Jonuco (JN), Casa Blanca 1 (CB1), Casa Blanca 2 (CB2) localidad tipo del sinónimo P. catarinae (^{Favela, 2004}), La Cantera (LC) y Cañón de Pinos (CP) (Cuadro 1, Figs. 1A, C).

Figura 1: A. Mapa con la ubicación de los siete sitios de estudio de P. remota (Little) D.K. Bailey & Hawksw. en el Noreste de México y con el rango de distribución de la especie (^{GBIF, 2021}), delimitado por el desierto Chihuahuense; B. Fisionomía de los árboles de la población tipo La Gavia (LG) de P. remota (^{Favela, 2004}); C. Fisionomía típica de los árboles del resto de las poblaciones catalogadas como P. catarinae Passini. 

Las poblaciones se distribuyen entre 1250-1496 m s.n.m. y, aunque las condiciones del hábitat son en general áridas, presentan importantes variaciones. Para la población LG distribuida al norte y alejada del resto de las poblaciones (entre 93-94 km de distancia lineal), el clima es de tipo muy seco semicálido con suelo Leptosol, temperatura entre 14-22 °C y precipitación de 200-500 mm (^{INEGI, 2009a}). Para las poblaciones del municipio Santa Catarina la temperatura oscila entre 10-22 °C, con precipitación de 200-900 mm (^{INEGI, 2009b}). En la población MV el clima es seco semicálido; en CP, semiseco templado y en el resto, seco templado. Los suelos son de tipo Litosol a excepción de MV donde es Xerosol Haplico (^{SEDUE, 2015}). Las poblaciones se asocian con especies de matorral desértico micrófilo (Acacia berlandieri Benth., Berberis trifoliolata Moric., Prosopis glandulosa Torr. y Gymnosperma glutinosum Less.) y matorral desértico rosetófilo (Agave lechuguilla Torr., A. berlandieri, Porlieria angustifolia A. Gray., Celtis pallida Torr. y Berberis trifoliata Torr.) (^{Mata et al., 2014}).

En cada población se recolectaron acículas de 15 a 22 árboles. Se obtuvo un total de 112 muestras, las cuales fueron guardadas en bolsas de plástico y etiquetadas con los siguientes datos de campo: número de individuo, especie, fecha y lugar de colecta. Así mismo, para cada árbol se tomaron datos de altura con la ayuda de un clinómetro digital (Forestry Suppliers, modelo Suunto PM5/66PC, Misisipi, EUA) y se midió el diámetro basal del tronco (DBA) con cinta diamétrica (Cuadro 1). Las muestras fueron transportadas en un contenedor con hielo al laboratorio de Genética de la Facultad de Agrobiología, Universidad Michoacana de San Nicolás de Hidalgo, y almacenadas en un ultracongelador (Thermo Fisher Scientific, modelo 931, Ohio, EUA) a -80 °C, hasta su procesamiento para la extracción de ADN.

Obtención de microsatélites de núcleo (SSRn)

El ADN se extrajo con el protocolo Miniprep CTAB 2× modificado para pinos (^{Vázquez-Lobo, 1996}). Se ensayaron 12 microsatélites de núcleo (SSRn): nueve reportados para P. strobus L. (^{Echt et al., 1996}) y tres para P. pinea L. (^{Pinzauti et al., 2012}), de los cuales cinco fueron reproducibles y polimórficos (Cuadro 2). La reacción de amplificación se hizo acorde con ^{Delgado Valerio et al. (2013)} y fue la siguiente: 2.5 μl buffer (1M Tris, pH 9.0), 2.0 μl dNTPs (4×, 100 mM), 0.5 μL de cada oligo (10 μM). El oligo derecho fue marcado con fluorescencia 6-FAM, VIC o NED, 11.3 μL de agua destilada, 1.5 μL, 50 mM MgCl₂ (4 mM), 0.2 μL de Taq polimerasa (1U) y 3.0 μL de DNA (50 ng/μl), con un volumen final de 23 μl.

Para las amplificaciones se usó un Termociclador Bio-Rad C1000™ (California, EUA) con un gradiente de temperatura dependiente de cada microsatélite (52-60 °C). El programa consistió en los siguientes pasos: 94 °C (5 min), gradiente de 18 ciclos a 95 °C (45 s), 52-60 °C (45 s) disminuyendo 0.5 °C cada ciclo y 72 °C (45 s), seguido de 20 ciclos similares al gradiente, pero con la menor temperatura de alineamiento alcanzada para cada microsatélite (42-50 °C), y una extensión final a 72 °C (5 min). Las amplificaciones fueron enviadas a un servicio externo a la Universidad de Illinois, EUA, para la obtención de los fragmentos con un secuenciador automático (Appied Biosystems, modelo 3730xl, Applied Biosystems, California, EUA), mientras que los genotipos fueron determinados con el programa Peack Scanner^TM v. 1 (Applied Biosystems, 2006).

Análisis de datos

Los estimadores de la variación genética fueron los siguientes: número de alelos por locus (n), promedio de alelos por locus (A), número efectivo de alelos por locus (A_e), heterocigosis observada (H_o) y esperada (H_e), e índice de endogamia (F). La desviación del equilibrio Hardy-Weinberg se estimó con el método no paramétrico de cadenas de Marcov tipo Monte Carlo (MCMC) con 100,000 permutaciones (^{Guo y Thompson, 1992}). Los análisis fueron obtenidos con el programa Arlequin v. 3.5.1.2 (^{Excoffier y Lischer, 2010}). La riqueza alélica (A_r) se estimó con el programa ADZ v. 1 (^{Szpiech et al., 2008}), que utiliza un enfoque de rarefacción para estandarizar las estimaciones al tamaño de muestra de la población más pequeña del conjunto de datos (^{El Mousadik y Petit, 1996}).

Para los cinco loci utilizados se hicieron pruebas exactas de desequilibrio de ligamiento a nivel de las poblaciones y una global con el método exacto de Fisher, obtenidas con el programa GENEPOP v. 4.4 (^{Rousset, 2008}). Así mismo, se determinó la frecuencia de alelos nulos con el programa FreeNa (^{Chapuis y Estoup, 2007}). El tamaño efectivo de la población (N_e) se obtuvo a partir del parámetro θ=4N_eµ, que es igual a cuatro veces el tamaño efectivo por la tasa de mutación (^{Kimura y Ohta, 1978}). Este parámetro (θ) se determinó a través del modelo de mutaciones de un paso SMM (^{Slatkin, 1995}), con el programa Arlequin v. 3.5.1.2 (^{Excoffier y Lischer, 2010}). La tasa de mutación que se asumió fue de 10^-3 por generación (^{Boys et al., 2005}; ^{Delgado Valerio et al., 2013}; ^{Rebolledo Camacho et al., 2018}).

Para inferir la consecuencia de procesos demográficos sobre el nivel de diversidad genética de las poblaciones, se probó la hipótesis de la reducción del tamaño de las poblaciones por el efecto de los cuellos de botella con dos métodos. El primero fue el método comparativo de la heterocigosidad obtenida (H_e) en el sentido de ^{Nei (1987)}, y la heterocigosidad esperada (H_eq) inferida del número de alelos observados para cualquier población y locus, bajo el supuesto de equilibrio entre deriva genética y mutación (^{Cornuet y Luikart, 1996}; ^{Piry et al., 1999}). Si se presenta un cuello de botella reciente, esta relación puede perderse y H_e podría ser temporalmente mayor que la H_eq, debido a que la pérdida de alelos es más rápida que la reducción de la heterocigosis (^{Luikart y Cornuet, 1998}; ^{Piry et al., 1999}). Para el análisis se usaron los modelos de mutación de un paso (SMM) y de mutación de dos fases (TPM), que asume mutaciones múltiples (30%) y simples (70%) (^{Di Rienzo et al., 1994}). La significancia fue obtenida con la prueba de Wilcoxon con 10,000 réplicas. Para este análisis se usó el programa Bottleneck v. 1.2.02 (^{Piry et al., 1999}).

En la segunda aproximación se usó el parámetro M radio (^{Garza y Williamson, 2001}); fue calculado como el radio medio del número de alelos (k) dividido por el intervalo del tamaño de los alelos (r) de cada loci, con el uso del programa Arlequin v. 3.5.1.2 (^{Excoffier y Lischer, 2010}). Un valor menor al Mc crítico=0.68 indica un evento de cuello de botella, debido a que durante el declinamiento de la población el número de alelos decrece más rápidamente que el intervalo en el tamaño de los alelos y disminuye por tanto el valor de M (^{Garza y Williamson, 2001}). Este método, a diferencia del anterior, es muy útil para detectar cuellos de botella ancestrales y mayores a 100 generaciones (^{Williamson-Natesan, 2005}).

Para la determinación de grupos genéticos se usó una aproximación Bayesiana con el programa STRUCTURE v. 2.3.4. (^{Pritchard et al., 2000}). Este programa infiere la estructura poblacional a partir de datos de genotipo de marcadores no ligados mediante un método de agrupamiento. Los modelos establecen que hay poblaciones (K grupos) que tienen frecuencias alélicas diferentes en cada locus analizado, donde los individuos se agrupan según la cercanía o grado de similitud genética. En este estudio se usó el modelo de mezcla genética, asumiendo información a priori de las poblaciones y frecuencias alélicas correlacionadas (^{Falush et al., 2003}). El número de K grupos analizados fue de 1 a 7 con 10 corridas independientes, 10,000 periodos de burnin y 100,000 iteraciones (MCMC). Subsecuentemente, el número de K grupos con la mayor probabilidad fue evaluado con el método complementario de ^{Evanno (Evanno et al., 2005}) e implementado en el programa STRUCTURE HARVESTER Web v. 0.6.94 (^{Earl y vonHoldt, 2012}). La asignación de la ancestría (q) de los individuos pertenecientes a cada grupo genético se obtuvo con el programa CLUMPP v. 1.1.2 (^{Jakobsson y Rosenberg, 2007}).

La estructura genética se obtuvo con un análisis jerárquico de varianza molecular (AMOVA) en dos niveles: entre poblaciones y entre individuos dentro de las poblaciones, con el uso del estadístico R _ST basado en el modelo de mutación de un paso (SMM) (^{Slatkin, 1995}). Este análisis también se hizo con los grupos genéticos obtenidos con el programa STRUCTURE (^{Pritchard et al., 2000}) en dos niveles: entre los K grupos genéticos y entre los individuos dentro de los K grupos. La significancia estadística se evaluó con 1000 permutaciones no paramétricas con el programa Arlequín v. 3.5.1.2 (^{Excoffier y Lischer, 2010}). El flujo genético entre pares de poblaciones se estimó con la formula Nm=(1/F_ST-1)/4 (^{Slatkin, 1993}).

Finalmente, para conocer si las poblaciones se ajustan a un modelo de aislamiento por distancia (IBD) (^{Slatkin, 1993}) se hicieron análisis pareados de Mantel con 10,000 permutaciones con el programa IBD v. 1.52 (^{Bohonak, 2002}). Las distancias genéticas obtenidas con los valores del estadístico R _ST entre pares de poblaciones (^{Slatkin, 1995}) se compararon con las distancias geográficas absolutas en kilómetros obtenidas a través de la transformación de Mercator con el programa QGIS v. 2.16 (^{Quantum GIS Development Team, 2016}).

Variación genética

Se presentó un total de 57 alelos en un intervalo de 8 (RPS2) a 17 (RPS20) por locus y un promedio de 5 alelos por locus. Los alelos más frecuentes fueron cuatro: 123pb (RPS34b), 220pb (RPS118), 143pb (RPS2) y 182pb (RPS120). El número de alelos de cada población fue relativamente diferente; el menor número se presentó en las poblaciones JN y CB2 (17 alelos) y el mayor número en las poblaciones MV (32) y LC (29). El promedio de alelos efectivos fue menor que la riqueza alélica entre las poblaciones: A_e de 1.962 (CB2) a 3.543 (MV) y A_r de 2.771 (CB2) a 4.178 (MV) (Cuadro 3). Se obtuvo un total de 22 alelos únicos, de dos a cinco por población y frecuencias entre 0.025 a 0.167. Las poblaciones de distribución más extrema presentaron un mayor número de alelos únicos: LG y MV ubicadas al norte de su distribución, cinco alelos únicos y LC sureña, cuatro. En la población JN no se detectó ningún alelo único. El promedio de la heterocigosis observada (H_o=0.375) fue menor a la esperada (H_e=0.557). El valor de H_o varió de 0.129 (CB2) a 0.550 (CP), y el de la heterocigosis esperada (H_e) de 0.427 (CB2) a 0.653 (MV y CP) (Cuadro 3). Todos los valores se desvían del equilibrio Hardy-Weinberg, pero con diferente grado de significancia estadística, lo que sugiere que existen diferencias en el déficit de la H_o respecto a la H_e entre las poblaciones y, por tanto, diferente intensidad de endogamia.

No se observó desequilibrio de ligamiento significativo para los pares de loci analizados a nivel de las poblaciones, ni a nivel global (P>0.05), lo cual indica que los loci segregan de manera independiente. Tampoco se presentaron alelos nulos con valores altos; solo en un locus (RPS20) con una proporción ligeramente mayor a 0.2 (0.228).

El promedio del índice de fijación fue alto y significativo (F_IS=0.259; P=0.0001). Las poblaciones CB2 (F_IS=0.554; P=0.001) y LG (F_IS=0.569; P=0.000) presentaron los valores más altos; las poblaciones CB1 (F_IS=0.198; P=0.011) y CP (F_IS=0.168; P=0.031), los más bajos y marginalmente significativos (Cuadro 3). El tamaño efectivo fue similar entre las poblaciones: valores de 375 (JN y CB1) a 427 (MV), incluso para las poblaciones con niveles importantes de endogamia (CB2, N_e=386 y LG, N_e=376). Por tanto, se puede considerar que todas las poblaciones contribuyen de manera similar al mantenimiento de la diversidad genética de la especie.

Demografía histórica

El análisis demográfico de la reducción reciente del tamaño de las poblaciones, basado en el modelo de exceso de heterocigosidad, muestra que en todas las poblaciones los valores de H_e y H_eq no son diferentes significativamente (P>0.05), lo que indica que no ha ocurrido una reducción reciente de su tamaño. Por el contrario, en tres de las poblaciones se presenta un déficit de H_e significativo (LG, MV y CB2), sugiriendo cuellos de botella ancestrales (Cuadro 4). El modelo de M radio muestra en todas las poblaciones valores menores al Mc crítico (0.68), los cuales se encuentran entre 0.186 (CB2) y 0.283 (MV) sugiriendo la existencia de declinamiento poblacional ancestral.

Grupos genéticos y aislamiento por distancia

El análisis Bayesiano de agrupamiento genético implementado en el programa STRUCTURE (^{Pritchard et al., 2000}), resultó en tres grupos genéticos con el mayor soporte estadístico (DeltaK=44.137) (Fig. 2). En el primer grupo genético K1 se incluyen las poblaciones LG y MV, donde LG tiene la distribución más al norte y alejada del resto. El segundo grupo K2 está conformado por las poblaciones JN, CB1 y CB2, de distribución más central, y LC y CP en el grupo genético K3, distribuidas en la zona más sureña. Sin embargo, todas las poblaciones muestran diferentes grados de ancestría de los tres grupos (valores de q entre 0.02 y 0.97). De los 112 individuos analizados, 64% (30 individuos) y 77% (45) de los genotipos pertenecen a los grupos genéticos K1 y K3 respectivamente, mientras que en el grupo K2 se presentó la mayor mezcla de genotipos de los tres grupos con 54% (29 individuos) de ancestría para este grupo. Debido a que no se presentaron diferencias genéticas entre la población de LG, considerada como tipo de P. remota, con las otras poblaciones catalogadas como P. catarinae y particularmente con la población MV, con la que conforma un mismo grupo genético (K1), se corrobora que P. catarinae es un sinónimo de P. remota.

Figura 2: Agrupamiento de las poblaciones de Pinus remota (Little) D.K. Bailey & Hawksw., en México, obtenido con el análisis Bayesiano de STRUCTURE (^{Pritchard et al., 2000}). A. Gráfico de barras que muestra el número de grupos K=3, con la mayor probabilidad obtenido con el método de ^{Evanno et al. (2005)}. Las barras verticales representan a cada individuo analizado en el estudio (112) y las barras con diferentes colores representa la cantidad de mezcla de los K grupos genéticos para cada individuo; B. Grupos genéticos (K) de mayor probabilidad con base al algoritmo DeltaK (^{Evanno et al., 2005}); C. Gráficos de pastel con la distribución de la proporción de ancestría (q) para cada uno de los K grupos genéticos obtenidos. 

El análisis jerárquico de varianza molecular (AMOVA) indica una diferenciación alta y significativa entre poblaciones (R _ST =0.158; P=0.000) y menor, pero significativa entre los grupos genéticos obtenidos con STRUCTURE (R _ST =0.082; P=0.000). En ambos casos la mayor variación se observó dentro de las poblaciones; 84.19% a nivel de las poblaciones y 91.80% con los grupos genéticos (Cuadro 5).

El promedio del Nm fue de 1.33 unidades, lo que sugiere un intercambio genético moderado. Los mayores valores de Nm se presentaron entre las poblaciones de JN y CB1, y LG y MV; los menores entre CB2 con LC y LG. El análisis de IBD no fue significativo (r²=0.084, P>0.05) e indica que no existe una asociación de la composición genética de las poblaciones con el área geográfica donde se distribuyen.

Variación genética y demografía histórica

Este estudio muestra que la variación genética de P. remota obtenida con el uso de SSRn es moderada y con una riqueza alélica similar entre las poblaciones. La H_e es comparable con otros estudios realizados con SSRn, como en P. oocarpa Schiede ex Schltdl. (H_e=0.642), P. patula Schltdl. & Cham y P. tecunumanii F. Schwerdtf. ex Eguiluz & J.P. Perry (H_e=0.567 y 0.646 respectivamente; ^{Dvorak et al., 2009}), P. strobus L. (H_e=0.740, ^{Zinck y Rajora, 2016}), P. mugo Turra (H_e=0.440; ^{Zukowska y Wachowiak, 2017}), P. sylvestris L. (H_e=0.772; ^{Bilgen y Kaya, 2017}), P. caribaea Morelet var. hondurensis (Sénécl.) W.H. Barret & Golfari (H_e=0.575; ^{Rebolledo Camacho et al., 2018}), P. pseudostrobus Lindl (H_e=0.515) y P. douglasiana Martínez (H_e=0.406) (^{Ramírez Enríquez et al., 2019}), lo que sugiere que contiene una diversidad genética que puede ser suficiente para el mantenimiento y posible adaptación de la especie a cambios ambientales.

Sin embargo, todas las poblaciones presentan niveles importantes de endogamia con un déficit significativo de heterócigos, por lo que no se encuentran en equilibrio Hardy-Weinberg. Se ha documentado que los valores altos de endogamia, y por tanto el déficit de heterócigos observados en las poblaciones, se pueden deber a la presencia de alelos nulos (alelos no expresados) que dan lugar a valores erróneos que desvían a las poblaciones de la panmixia (^{Shinde et al., 2003}; ^{Ramírez Enríquez et al., 2019}). No obstante, en el análisis para la detección de alelos nulos solo se encontró un locus (RPS20) con frecuencias superiores a 0.2 (0.228), que es un valor admisible para no subestimar significativamente la endogamia o la estructura genética de las poblaciones (^{Chapuis y Estoup, 2007}). Además, los loci han sido usados en varios estudios con P. strobus y no se reporta la presencia de alelos nulos (^{Marquardt y Epperson, 2004}; ^{Chhatre y Rajora, 2014}; ^{Zinck y Rajora, 2016}), ni presentan desequilibrio de ligamiento genético (^{Chhatre y Rajora, 2014}). Este comportamiento concuerda con los resultados obtenidos en el presente trabajo; por ello se puede considerar que la desviación del equilibrio Hardy-Weinberg obtenida está fundamentalmente determinada por el efecto de la endogamia y/o reproducción entre parientes cercanos. Reproductivamente para esta especie se ha reportado una capacidad germinativa muy baja (0.1 a 15.7%), así como una escasa producción de semillas viables (11 semillas/cono), lo cual se asocia al efecto de la endogamia (^{Lemus, 1999}) y da soporte a los resultados obtenidos en este estudio.

Así mismo, los tamaños efectivos son reducidos y menores a los reportados en otras especies de pinos como P. caribaea var. hondurensis (N_e entre 402-868) (^{Delgado et al., 2011}), con solo dos poblaciones fragmentadas y restringidas al sureste de México. Lo mismo sucede con algunas poblaciones de esta variedad y de P. caribaea Morelet var. caribaea y P. caribaea Morelet var. bahamensis (Griseb.) W.H. Barret, distribuidas en la cuenca del Caribe (N_e entre 161-537) (^{Rebolledo Camacho et al., 2018}) y de P. pseudostrobus, P. oocarpa y P. douglasiana, distribuidas en la Cuenca del Cupatitzio en Michoacán, México (N_e entre 183-2104, ^{Ramírez Enríquez et al., 2019}). Es muy probable que la fragmentación de las poblaciones de P. remota no sea reciente, pues los efectos de la deriva genética y/o de la endogamia han influido a través de su historia evolutiva y se reflejan en la disminución de la variabilidad genética, y por ende en los estimados de sus tamaños efectivos.

^{Garza y Williamson (2001)} mencionan que la permanencia de las poblaciones relativamente pequeñas y aisladas durante un largo tiempo, puede deberse a cuellos de botella antiguos (>100 generaciones). Lo anterior implica que los valores bajos de M radio que se observan en las poblaciones de P. remota y que son menores al Mc crítico (0.68; ^{Garza y Williamson, 2001}) pudieron ser causados por cuellos de botella de manera independiente y prolongada, manteniendo a las poblaciones relativamente aisladas y con la consecuente pérdida de variantes alélicas. Además, el análisis de heterocigosidad muestra a cuatro de seis poblaciones con valores de H_eq significativos, pero para déficit de heterocigosidad (H_eq>H_e), lo que sugiere procesos ancestrales de reducción poblacional (^{Cornuet y Luikart, 1996}; ^{Luikart y Cornuet, 1998}), como es el caso de P. resinosa Aiton en el Noreste de EUA (^{Boys et al., 2005}) y de las dos poblaciones de P. caribaea var. hondurensis distribuidas en México (^{Delgado et al., 2011}).

Estructura genética e historia evolutiva

Los tres grupos genéticos obtenidos con el análisis Bayesiano STRUCTURE pueden ser explicados por el movimiento de variantes alélicas que comparten entre las poblaciones de cada grupo. El grupo K2 resultó ser el de mayor mezclada; es decir, con un mayor Nm entre las poblaciones que lo conforman (JN, CBA y CB2) y que están ubicadas en la zona central del área de estudio. En contraste, se observa una tendencia a la diferenciación de los grupos extremos (K1 y K3), que muestran una ancestría superior a 64% de pertenencia a cada grupo genético. Curiosamente la población LG considerada como tipo de P. remota y que tiene la distribución más al norte y se encuentra alejada del resto de las poblaciones y aislada por una barrera geográfica importante (montaña rocosa llamada La Muralla, Coahuila), se incluye en el grupo K1 junto con una de las poblaciones consideradas como P. catarinae (MV). Lo anterior permite corroborar que es un sinónimo de P. remota, pues comparte variantes alélicas con esta población y diferentes grados de ancestría con el resto de los grupos genéticos obtenidos (ver Fig. 2).

Acorde con el estadístico R _ST , la estructura genética tanto a nivel de las poblaciones como a nivel de los grupos genéticos fue alta y significativa. Esto indica que el efecto del proceso de mutación es similar o mayor con respecto a la migración y puede ser el recurso primario más probable de la estructura genética observada (^{Balloux y Goudet, 2002}; ^{Hardy et al., 2003}) y con una menor variación genética entre las poblaciones que la observada dentro de las poblaciones, como es común que suceda en especies de coníferas (^{Ledig et al., 2001}; ^{Delgado y Piñero, 2008}). Aunque existen pocos estudios que han estimado la estructura genética con SSRn en pinos, el resultado es similar a los obtenidos en P. strobus (F_ST=0.104; ^{Zinck y Rajora, 2016}), P. pseudostrobus (F_ST =0.152), P. douglasiana y P. oocarpa (F_ST=0.131) para ambas especies (^{Dvorak et al., 2009}; ^{Ramírez Enríquez et al., 2019}); también a los reportados con el uso de otros marcadores genéticos como las isoenzimas (P. pinceana Gordon, F_ST=0.152, ^{Ledig et al., 2001}; F_ST=0.240, ^{Molina-Freaner et al., 2001}), o los ISSR (P. cembroides Zucc., F_ST=0.260, ^{Fuentes-Amaro et al., 2019}). Este comportamiento se ha observado para la mayoría de las especies de coníferas mexicanas, pues su distribución es fragmentada y aislada, lo que reduce los niveles de flujo genético y en consecuencia se incrementa la diferenciación genética con la fijación de distintos alelos en las poblaciones (^{Ledig et al., 2001}; ^{Delgado et al., 2002}).

En efecto, el promedio del Nm (1.33 unidades) sugiere un intercambio genético moderado y relativamente menor al obtenido en otras especies de pino (p. ej., P. pseudostrobus y P. douglasiana con valores de Nm de 1.40 y 1.66 respectivamente; ^{Ramírez Enríquez et al., 2019}; P. johannis M.-F. Robert, Nm=4.0, ^{García-Gómez et al., 2014}; P. cembroides, Nm de 1.52-1.75, ^{Fuentes-Amaro et al., 2019}). Las coníferas, en general, tienen mecanismos de dispersión por viento tanto de polen como de semillas, por lo que el estimado del flujo genético generalmente es superior a 1.5 unidades (^{Hamrick et al., 1995}). En el estudio de ^{Ledig et al. (2001}) se reporta un intervalo de 4 a 17 unidades para este grupo de plantas. Con base en ello, consideramos que el valor obtenido en P. remota es moderado y similar al observado en otras especies de coníferas de distribución restringida y endémicas de México como Picea martinezii T.F. Patterson (Nm=0.99-1.1; ^{Ledig et al., 2000}), P. greggii Engelm. ex Parl. (Nm=1.16; ^{Parraguirre Lezama et al., 2002}) y P. pinceana (Nm=1.39, Ledig et al., 2001). Además, como se sabe, las semillas de las especies de pinos piñoneros (grupo al cual pertenece P. remota) no tienen alas para su dispersión por medio del viento y la realizan fundamentalmente con ayuda de dispersores, como son las aves de la familia Corvidae (^{Gernandt et al., 2003}). Por ello, existe una mayor probabilidad de conformar fragmentos o parches asociados con las rutas de movimiento o migración de los dispersores, lo cual también influye en el incremento de la estructuración genética y disminución del Nm entre las poblaciones (^{Delgado et al., 2002}; ^{Favela, 2004}; ^{Ellstrand, 2014}). Aunado a lo anterior, las especies de pinos tienen capacidad para autofecundarse, pues son en su mayoría monoicas y auto-compatibles (^{Williams, 2008}). Por tanto, se puede considerar que debido a las características reproductivas y de historia de vida de P. remota, las poblaciones muestran niveles importantes de endogamia, con poca posibilidad de intercambiar variantes alélicas, con presencia de alelos únicos y de alelos compartidos con frecuencias bajas que acentúan la diferenciación genética entre ellas, representado en la conformación de los tres grupos genéticos obtenidos con el análisis Bayesiano.

No se presentó un patrón de IBD, lo cual indica que la distribución de la variación genética no está asociada con el área geográfica donde se distribuyen las poblaciones. Esto sugiere, nuevamente, que esta diferenciación genética podría estar en función de factores ambientales locales o microambientales, como pueden ser climáticos (grado de aridez), edáficos (tipo calizo), altitudinales (^{Favela, 2004}), o de barreras geográficas naturales (^{Garrido-Garduño y Vázquez-Domínguez, 2013}), y que probablemente también contribuyen a las diferencias morfológicas observadas en la fisionomía de los árboles (^{Snajberk y Zavarin, 1986}). La no existencia de un patrón de IBD se ha presentado en otras especies de la subsección Cembroides, como P. johannis y P. cembroides (^{Cuenca, 2003}; ^{Fuentes-Amaro et al., 2019}). Por ejemplo, P. cembroides tiene una distribución amplia; sin embargo, al tener poblaciones aisladas y fragmentadas se presenta ese comportamiento: valores altos de estructura genética, tamaños efectivos intermedios y una asociación geográfica nula, lo que es atribuible en gran medida a la deriva génica, a los sistemas reproductivos o a la presión de selección en respuesta a diferentes componentes ambientales (^{Ledig et al., 2001}; ^{Delgado et al., 2002}).

Acorde con los resultados expuestos se puede apoyar un escenario evolutivo para P. remota: que estuvo ampliamente distribuida y que, debido a su fragmentación por algún evento geológico o ambiental, las poblaciones quedaron aisladas y diferenciadas genéticamente. Los trabajos realizados en el desierto de Chihuahua sobre la vegetación del Pleistoceno tardío, conocido en América como glaciaciones del Wisconsin (100 mil a 10 mil años atrás), establecen de manera concluyente que P. remota estaba distribuida ampliamente en esta región con una extensión que abarcaba hasta el sur de Durango (^{Lanner y Van Devender, 1981}). Una aproximación del tiempo en el cual las poblaciones pudieron haber declinado se puede obtener al multiplicar el N_e promedio obtenido para la especie, por el intervalo de generaciones propuesto por ^{Cornuet y Luikart (1996)} de 2N_e a 4N_e para cuellos de botella recientes (menores a 80 generaciones). En P. remota las estimaciones son aproximadamente de 974 a 1588 generaciones atrás (mucho más de 100 generaciones), que al asumir un tiempo generacional a la primera reproducción (15 años; ^{Romero Manzanares et al., 1996}) sugieren que el declinamiento poblacional pudo haber ocurrido entre 11,910 y 23,820 años atrás, tiempo asociado claramente con las últimas glaciaciones en América. Este proceso pudo tener efectos importantes en la dinámica evolutiva de la especie. Por un lado, al quedar las poblaciones fragmentadas y aisladas por eventos de contracción de las poblaciones (lo cual se sustenta con los resultados del análisis de cuellos de botella obtenidos en este estudio) se promovió la diferenciación genética dando lugar a la disminución de la variabilidad genética y al incremento de los niveles de endogamia, y que debido a los mecanismos reproductivos (monoica y auto-compatible) y de dispersión limitada de la especie tuvo efectos negativos en su éxito reproductivo (^{Lemus, 1999}). Por otro lado, se promovió la divergencia adaptativa probablemente por presión de selección asociada a condiciones ambientales locales, lo que les ha permitido subsistir hasta la actualidad. Por tanto, las poblaciones de P. remota podrían considerarse como relictos remanentes del Wisconsin (^{Snajberk y Zavarin, 1986}). No obstante, aunque esta investigación hace un aporte importante al conocimiento genético y evolutivo de la especie, es necesario hacer estudios complementarios que incluyan otras poblaciones reportadas más al norte de México y en el sur de Texas, EUA (^{Farjon y Styles, 1997}), con la finalidad de comparar y verificar si han evolucionado por procesos demográficos similares a los obtenidos en este trabajo y tener un panorama más completo de su historia evolutiva. Así mismo, sería interesante explorar el uso de marcadores selectivos para probar hipótesis relacionadas con la adaptación local de los individuos y/o poblaciones a las variaciones microambientales en las que habita la especie.

Recomendaciones para la conservación

Los resultados obtenidos en este trabajo también son de utilidad para delinear algunas acciones puntuales para la conservación de la especie. Primeramente, se sugiere iniciar con la protección formal de toda el área de su distribución, con la finalidad de incrementar su tamaño censal de manera natural y que a mediano plazo se logre disminuir el índice de endogamia. De manera prioritaria se debe de proteger mínimamente una de las poblaciones incluidas en cada grupo genético obtenido con el análisis Bayesiano: del grupo K1 a la población MV, ubicada en la parte norte del área de estudio; del grupo K2 a la población CB1, ubicada en la zona central, y del grupo K3 a la población CP, de la zona sur. Estas poblaciones contienen el mayor número de alelos y de alelos únicos, niveles altos de H_o y H_e y una endogamia moderada. Así mismo, se requiere producir planta con semilla proveniente de estas poblaciones más diversas genéticamente para dirigir los esfuerzos de reforestación a las poblaciones con valores bajos de variación genética y tratar de incrementarla. Por ejemplo, la población CB2, que presentó los valores más bajos de H_o y más altos de endogamia, así como la población LG (población tipo de P. remota) con la distribución más extrema, al presentar un valor alto de endogamia y el segundo valor de H_o más bajo.

En general, consideramos que el rescate de la especie debe basarse en la conservación in situ dirigida a la conservación mínima de las tres poblaciones mencionadas (MV, CB1 y CP), que son clave para mantener el acervo genético de la especie. Si bien P. remota no tiene una utilidad económica inmediata, sus poblaciones proporcionan una gran cantidad de bienes y servicios ambientales (madera para construcción, oxígeno, producción de leña, regulación del clima, retención de agua, resguardo para la fauna silvestre, protección de la biodiversidad y espacio escénico y de recreación, entre otros). Particularmente las semillas o piñones son importantes como alimento para aves y roedores, y para el hombre, aunque su consumo es muy local. Los árboles tienen gran potencial para la restauración ecológica de zonas degradadas por su capacidad para crecer en bajas elevaciones y en condiciones locales heterogéneas y adversas, por lo que se podrían plantar en muchas áreas que ahora son completamente improductivas (^{Perry, 1991}). Además, la gran plasticidad genética que se ha observado en las características fisonómicas de los árboles que conforman las poblaciones (^{Snajberk y Zavarin, 1986}; ^{Favela, 2004}) podría ayudar a incrementar el éxito de las reforestaciones y plantaciones con esta especie. Las actividades propuestas para la conservación de P. remota son un gran reto, que solo se logrará con la participación coordinada de los ejidos y comunidades poseedoras de los bosques, así como de las autoridades locales, estatales y federales encargadas de salvaguardar los recursos forestales del país.