Área de estudio

Se exploraron cinco poblaciones del centro-occidente del Cinturón Volcánico Transmexicano, cuatro en Michoacán (Arocutín, San Andrés Coru, Nuevo San Juan Parangaricutiro y Zacapu) y una en la Sierra de Manantlán, Jalisco (Cuadro 1; Fig. 2). Estas localidades se encuentran en diferentes formaciones volcánicas. Las de Michoacán pertenecen al campo Michoacán-Guanajuato del Cinturón Volcánico Transmexicano. El malpaís de Zacapu está compuesto por varios flujos de lava agrupados en dos periodos con una edad de 27,000 a 21,300 a. C. y 1479 a 900 a. C (^{Reyes-Guzmán, 2018}). Al flujo de lava (El Capaxtiro) donde se realizó la colecta en este estudio, se le ha asignado una edad de entre 1500 a 2500 años de antigüedad (^{Reyes-Guzmán, 2018}). Por otra parte, en el pedregal Arocutín, la zona de muestreo fue en el cerro La Taza, que se encuentra al suroeste del Lago de Pátzcuaro, y se le ha asignado una edad aproximada de 8430 años (^{García-Quintana et al., 2016}; ^{Osorio-Ocampo, 2018}). Los pedregales cercanos a Nuevo San Juan Parangaricutiro fueron originados por varios eventos de actividad volcánica durante el Plioceno, Pleistoceno y Holoceno (^{Bocco, 1998}), aunque se sabe que 50% de todo ese territorio fue cubierto por cenizas o flujos de lava provenientes del volcán Paricutín durante su erupción de 1943 a 1952 (^{Medina et al., 2000}). En la localidad de San Andrés Coru, el volcán El Metate es considerado el más joven del campo volcánico Michoacán-Guanajuato, iniciando su actividad en el año 1250 d. C. (^{Chevrel et al., 2016}). Los individuos de B. stormiae se colectaron en el flujo de mayores dimensiones (flujo 6). En el estado de Jalisco, la formación de la Sierra de Manantlán se encuentra en la confluencia de la Sierra Madre del Sur y el Cinturón Volcánico Transmexicano (^{INE, 2000}). La zona de colecta de B. stormiae se ubica en el occidente de la Sierra de Manantlán en la Unidad Fisiográfica Las Joyas, la cual es la parte más alta y conservada, formada por roca volcánica intrusiva y extrusiva, que es producto del vulcanismo del Mioceno y Pleistoceno (^{Martínez y Ramírez, 1998}).

Figura 2: Ubicación geográfica de las poblaciones muestreadas de Balmea stormiae Martínez en Michoacán y Jalisco, México. 1. Arocutín; 2. Manantlán; 3. Nuevo San Juan; 4. San Andrés Coru; 5. Zacapu. 

Muestreo de material vegetal

Se incluyeron muestras de la mayor cantidad de individuos que estuvieran separados entre sí por al menos 15 m, pero dada la dificultad de desplazamiento dentro de los pedregales, solo se muestrearon entre 8 a 15 individuos por localidad (Cuadro 1). Las muestras por individuo consistieron en una rama de 35 cm (permiso de colecta SEMARNAT SGPA/DGVS/003505/18), de las cuales una hoja sin rastros de herbívoros o de alguna afectación fúngica se preservó en gel de sílice (J. T. Baker, Estado de México, México). Una vez deshidratadas se mantuvieron en un congelador a -20 °C en el laboratorio de Sistemática Molecular de Plantas, Facultad de Biología, Universidad Michoacana de San Nicolás de Hidalgo, donde se realizó la extracción del ADN con el kit de Invisorb ® Spin Plant Mini Kit (STRATEC, Birkenfeld, Alemania).

Amplificación y secuenciación de ADN

Se amplificó el intrón de cloroplasto rps16 siguiendo a ^{Oxelman et al. (1997)} con los cebadores rpsF y rpsr2 (Cuadro 2). Además, se incluyó una secuencia de este marcador obtenida de GenBank (^{GenBank, 2022}; Cuadro 3) de un individuo reportado para Colima por ^{Manns y Bremer (2010)}.

Las amplificaciones del ETS se realizaron siguiendo a ^{Becerra (2003)} utilizando su cebador BUR1 y el ETS-iR1 (Cuadro 2). Las PCRs se realizaron en un termociclador TC-412 (Techne, Staffordshire, UK). Los fragmentos de ADN resultantes se enviaron a secuenciar a Macrogen Inc., Corea del Sur. Los electroferogramas “forward” y “reverse” se revisaron visualmente y se ensamblaron con Sequencher v. 4.8 (Gene Codes Corp. Ann Arbor, Michigan, EUA). Las secuencias obtenidas se alinearon manualmente en el programa BioEdit v. 7.1.9 (^{Hall, 1999}). Los sitios variables detectados fueron revisados cuidadosamente para descartar errores.

En el caso de los ISSR se realizaron pruebas preliminares con 28 cebadores de los desarrollados por el Laboratorio de Biotecnología de Universidad de Columbia Británica, en Vancouver, British Columbia, Canadá. De estos se eligieron con base en la obtención de bandas claras y su reproducibilidad en todos los ensayos, a los oligonucleótidos denominados UBC807, UBC811, UBC814, UBC825, UBC841, UBC842 y UBC844 (Cuadro 2), anteriormente utilizados en varios estudios (^{Lau et al., 2005}; ^{Li y Jin, 2008}; ^{Muthusamy et al., 2008}; ^{Datta et al., 2010}; ^{Gama-Maia y Torres, 2016}; ^{George et al., 2017}). Las condiciones de amplificación siguieron el procedimiento realizado por ^{Santillán-Mendoza et al. (2018)}. Los productos de los amplicones se corrieron en geles de agarosa al 1.5%, con TBE 0.5×, a cinco voltios por centímetro (5 V/cm) y un marcador de peso molecular de 100 pb (New England Biolabs, UK). Los geles se sometieron a postinción en una dilución al 3× del colorante GelRed (Biotum, Inc, EUA) por 30 minutos.

Análisis de secuencias de ADN

Las ambigüedades en las secuencias de rps16 fueron tratadas como datos faltantes y los haplotipos se establecieron en DNAsp v. 5 (^{Librado y Rozas, 2009}). En el caso de los datos nucleares ETS, se evidenciaron los individuos heterocigotos ubicando dobles picos en los electroferogramas de las secuencias. Se detectaron y eliminaron los bloques con recombinación mediante IMgc (^{Woerner et al., 2007}) y las ambigüedades encontradas se resolvieron con el método Phase (^{Stephens et al., 2001}) del programa DNAsp v. 5 (^{Librado y Rozas, 2009}). Con este mismo programa se calcularon los índices de diversidad nucleotídica de ^{Nei (π, Nei, 1987}), diversidad haplotípica (Hd,^{Nei, 1987}), diferenciación genética (G_ST, ^{Nei, 1973}) y flujo genético (Nm, ^{Wright, 1951}), tanto para rps16 como para ETS. Las relaciones genealógicas entre haplotipos se infirieron con redes de “median-joining” en PopART v. 1.7 (^{Bandelt et al., 1999}).

La baja variabilidad genética obtenida en rps16 de este muestreo impidió realizar analisis adicionales con estas secuencias de cloroplasto. En el caso de ETS, para tratar de comparar variabilidad genética, compensando el diferente tamaño del muestreo entre las poblaciones, se obtuvo la riqueza de haplotipos (RH) con el programa SPAGeDi v. 1.5a (^{Hardy y Vekemans, 2002}). Los análisis que se describen a continuación se realizaron con las secuencias de ETS, considerando dos regiones (Jalisco y Michoacán) y las cinco poblaciones, una en Jalisco y cuatro dentro de Michoacán. Se realizaron análisis “mismatch distribution” y de partición de la variabilidad genética (AMOVA), así como la obtención de las F_ST pareadas en Arlequin v. 3.5.2 (^{Excoffier y Lischer, 2010}). Se evaluó la estructura filogeográfica comparando G_ST y N_ST en PERMUT-CPSSR v. 2.0 (^{Pons y Petit, 1996}). Las pruebas de neutralidad con Fs de ^{Fu (Fu, 1997}) y D de ^{Tajima (Tajima, 1989}) se realizaron en DNAsp v. 5 (^{Librado y Rozas, 2009}) y Arlequin v. 3.5.2 (^{Excoffier y Lischer, 2010}). Se evaluó el aislamiento por distancia con una prueba de ^{Mantel (1967)} en GenAlEx v. 6.5 (^{Peakall y Smouse, 2006}), usando los valores de las F_ST pareadas y las distancias geográficas lineales, realizando dos analisis, uno incluyendo todas las poblaciones y otro excluyendo Manantlán, la localidad más alejada.

Análisis de datos ISSR

Las bandas resultantes de las amplificaciones se codificaron como caracteres binarios (^{Wolfe, 2005}). Utilizando Popgene v. 1.31 (^{Yeh et al., 1999}), se estimó el número de alelos observados (Na), el número de alelos efectivos (Ne) (^{Kimura y Crow, 1964}), la diversidad genética de ^{Nei (h) (Nei, 1973}), el porcentaje de loci polimórficos (%P), el índice de ^{Shannon (I) (Shannon, 1948}) y la diferenciación genética (G_ST) (^{Nei, 1973}). El flujo genético se calculó como: Nm=0.5 (1-G_ST) / G_ST, con el supuesto de que todas las poblaciones tienen el mismo flujo genético. Con GenAlEx v. 6.5 (^{Peakall y Smouse, 2006}), se obtuvieron las distancias genéticas entre las poblaciones con el método de ^{Nei (1972)}. Se realizó el AMOVA con dos regiones (Jalisco y Michoacán) y poblaciones dentro de Michoacán. La prueba de ^{Mantel (1967)} se efectuó entre las distancias genéticas de Nei y las distancias geográficas lineales entre las poblaciones, realizando dos análisis, uno incluyendo todas las poblaciones y otro excluyendo Manantlán, la localidad más alejada.

La estructura genética se determinó con STRUCTURE v. 2.3 (^{Pritchard et al., 2000}) considerando los datos como marcadores dominantes (^{Falush et al., 2007}) usando el modelo de ancestría mezclada (Admixture model) y frecuencias alélicas correlacionadas (Allele Frequencies Correlated). El análisis se hizo con un periodo burn-in de 100,000, longitud de corrida de 1,000,000 generaciones de MCMC (cadenas Markovianas de Monte Carlo) y 10 réplicas independientes para cada valor de K (de 1 a 6). Se determinó el número adecuado de grupos con el valor de delta K de acuerdo con ^{Evanno et al. (2005)} en Structure Harvester v. 0.6.94 (^{Earl y VonHoldt, 2012}).

Diversidad genética y redes de haplotipos con rps16 y ETS

El alineamiento de rps16 consistió en 73 individuos y 716 posiciones, presentando baja diversidad haplotípica (Hd=0.335±0.056) y baja diversidad nucleotídica (π=0.0014±0.00024). Se encontraron dos haplotipos separados por tres sustituciones. El haplotipo con mayor frecuencia (H1) se presenta en 57 muestras de las cuatro poblaciones de Michoacán (Cuadro 3). El haplotipo menos frecuente se encontró en 16 muestras, 15 de Jalisco y una de Colima (Cuadro 3). Consecuentemente la red de haplotipos separa las poblaciones de Michoacán de las de Jalisco/Colima (Fig. 3C).

Figura 3: A. distribución geográfica de los ribotipos del ETS obtenidos de cinco poblaciones de Balmea stormiae Martínez del centro-occidente de México, donde los colores de la derecha diferencian los 16 ribotipos encontrados; B. red de los ribotipos referidos; C. red de haplotipos con rps16, el haplotipo H1 se encuentra en las poblaciones de Michoacán y el H2 en Manantlán (Jalisco) y Colima. 

En los datos de ETS se incluyó cada elemento del par de secuencias de los heterocigotos detectados, por lo que el alineamiento resultó de 78 secuencias, 560 posiciones y 17 sitios polimórficos (S). La diversidad haplotípica fue alta (Hd=0.754 ± 0.041) y la diversidad nucleotídica baja (π=0.0021). Las poblaciones con mayor diversidad haplotípica y nucleotídica son Arocutín y Zacapu; en cambio en Manantlán, Jalisco se presentó la menor diversidad (Cuadro 1). En el caso de la estimación de la riqueza de haplotipos (RH) se obtuvo que Zacapu presenta la mayor riqueza de haplotipos (RH=6), seguido de Arocutín (RH=4.85), mientras que Manantlán presentó la menor (RH=1.92). Se obtuvieron 16 ribotipos (H) dispuestos en una red de haplotipos en forma de estrella: al centro el ribotipo más común (H2), presente en 55 secuencias de Michoacán, seguido por un ribotipo (H15) presente en 20 secuencias de Jalisco (Cuadro 3). Al menos un ribotipo único se encontró en todas las poblaciones estudiadas. No se presentaron ribotipos compartidos entre Michoacán y Jalisco (Figs. 3A, B).

Estructura genética

El valor de N_ST no fue significativamente mayor que el G_ST, lo que indica ausencia de estructura filogeográfica (N_ST=0.349, G_ST=0.182, P>0.05). El AMOVA con los datos ETS indicó que 51.76% de la variación fue explicada por las diferencias dentro de las poblaciones, mientras que 7.57% de la diversidad genética se encontró entre poblaciones y aunque la variación entre regiones es de 40.67%, este valor no fue significativo (Cuadro 4). Los valores más altos de las F_ST pareadas con ETS (0.35 a 0.71, P<0.05) se presentaron entre la población de Manantlán (Jalisco) y las cuatro poblaciones de Michoacán (Cuadro 5). En Michoacán, Arocutín presentó los mayores niveles de diferenciación al compararla con las demás poblaciones, con valores de F_ST moderados y bajos (entre 0.15 y 0.099, P<0.05). Por otro lado, San Andrés Coru y Nuevo San Juan, localidades cercanas geográficamente, mostraron una F_ST baja (F_ST=0.073, P<0.05), mientras que Zacapu presentó F_ST bajos respecto a San Andrés Coru (0.024, P>0.05) y Nuevo San Juan (0.0085, P>0.05). El flujo genético global, incluyendo todas las poblaciones, fue relativamente bajo (Nm=0.81), mientras que el flujo genético sin incluir la población de Jalisco fue moderado (Nm=8.95).

Historia demográfica

Los valores globales de D de Tajima (-1.88, P<0.05), de Fs de Fu (-11.45, P<0.05) y la distribución unimodal en el “mismatch distribution” sugirieron que hay señal de una expansión demográfica reciente en B. stormiae. A nivel de poblaciones, se obtuvo D de Tajima y Fs de Fu negativas y significativas para Zacapu (-1.22 y -3.44, P<0.05) y para San Juan Nuevo (-1.54 y -2.75, P<0.05), sugiriendo que hay desviaciones de la neutralidad esperada y expansión poblacional. Los resultados de la prueba de Mantel, incluyendo todas las poblaciones, indicaron una alta relación positiva (no significativa) entre la distancia geográfica y la distancia genética calculada (R²=0.790, P>0.05), obteniendo una mucho menor correlación si se excluye la población de Jalisco (R²=0.222, P>0.05).

Diversidad genética con ISSR

La matriz ISSR produjo 97 loci (bandas) oscilando entre 11 (UBC807) y 22 (UBC841) por marcador y se generaron 81 genotipos diferentes. Se encontró un alto porcentaje de loci polimórficos (80.41%). La diferenciación genética fue moderada, como indicó la F_ST (0.239, P<0.001), presentándose los valores más altos entre la población de Manantlán y las de Michoacán (Cuadro 5). La población de Arocutín fue la más diversa, ya que presentó el mayor porcentaje de loci polimórficos (60.82%), el mayor número de alelos observados (Na=1.6082) y alelos efectivos (Ne=1.3953), así como los valores más altos de diversidad genética de Nei (h=0.2271) y de Shannon (I=0.3359). Por el contrario, la población de Manantlán fue la que presentó menores niveles de diversidad (Cuadro 6).

El AMOVA con ISSR indicó que 76% de la variación fue explicada dentro de las poblaciones (F_ST=0.239, P<0.001), seguido por 17% entre las regiones de Michoacán y Jalisco (F_CT=0.169, P<0.001), mientras 7% de la variación restante fue entre poblaciones (F_SC=0.084, P<0.001) (Cuadro 4). El delta K del método de Evanno en STRUCTURE Harvester indicó que el número óptimo de grupos genéticos fue cuatro. La estructura genética coincidió con la distribución geográfica de B. stormiae en Jalisco y Michoacán, pero con ancestría mezclada en la mayoría de los individuos (Figs. 4A, B). La prueba de Mantel con ISSR, incluyendo todas las poblaciones, indicó una correlación positiva significativa (R²⁼0.921, P<0.05), contrastando con la nula correlación obtenida (R²=0.054, P>0.05) para esta misma prueba si se excluye la población de Jalisco.

Figura 4: Agrupación de las poblaciones de Balmea stormiae Martínez de Michoacán y Jalisco obtenidas con STRUCTURE. A. distribución de los grupos genéticos obtenidos; B. asignación de los grupos genéticos a partir de los datos ISSR de las poblaciones estudiadas para el delta K=4, donde cada color representa un diferente grupo genético. 

Estructura genética e historia demográfica

Los resultados obtenidos en el AMOVA con ETS e ISSR indican que los porcentajes de variación genética coinciden en ser bajos entre poblaciones, intermedios entre regiones y altos dentro de las poblaciones, indicando flujo genético entre poblaciones (al menos entre las poblaciones de Michoacán) y bajo entre regiones. Mientras tanto, la alta variación genética dentro de las poblaciones puede ser consecuencia de una expansión demográfica, detectada por la D de Tajima, Fs de Fu y la “mismatch distribution”, con datos del ETS.

La ausencia de estructura filogeográfica estimada con G_ST y N_ST con los datos del ETS es congruente con su red de haplotipos, que muestra una estructura filogeográfica muy simple reflejando principalmente la separación entre las poblaciones de Michoacán y Jalisco. Por otro lado, la combinación de diversidad haplotípica alta y diversidad nucleotídica baja de este marcador en conjunto con la red de haplotipos con forma de estrella, en la que hay un haplotipo notoriamente predominante, también sugiere un crecimiento poblacional reciente (^{Avise, 2000}).

Es de resaltar que los altos valores de loci polimórficos, diversidad genética de Nei y diversidad de Shannon obtenidos con ISSR para la especie en estudio son similares a los reportados en otras especies de rubiáceas, como Psychotria hastisepala Müll. Arg. y Morinda tomentosa B. Heyne ex Roth (^{Arya et al., 2014}; ^{Silva et al., 2014}). Para P. hastisepala por ejemplo, se sugirió que los niveles de variabilidad genética se atribuyen a su autoincompatibilidad y reproducción vegetativa (^{Silva et al., 2014}). Estos factores deberán ser descartados para B. stormiae en estudios venideros.

Otra hipótesis por comprobar en futuros estudios es la concordancia de los niveles de variación genética para las poblaciones de Michoacán obtenidos tanto con ISSR como con la π del ETS, con la orogénesis reportada para los flujos de lava de mayor a menor antigüedad: Arocutín, Zacapu, San Andrés Coru y Nuevo San Juan (^{Bocco et al., 1998}; ^{Martínez y Ramírez, 1998}; ^{Chevrel et al., 2016}; ^{García-Quintana et al., 2016}; ^{Osorio-Ocampo, 2018}; ^{Reyes-Guzmán, 2018}). Es decir, las poblaciones con mayor diversidad genética (Arocutín y Zacapu) se encuentran en los flujos de lava más antiguos, mientras que las poblaciones menos diversas se encuentran sobre los pedregales más recientes. Si esto se comprueba en estudios posteriores, será congruente con el argumento de que la variación genética de poblaciones ancestrales es mayor que la presente en las poblaciones más recientes, ya que estas últimas han tenido menos tiempo para acumular diferencias (^{Mastretta-Yanes et al., 2015}). Sin embargo, al parecer esto no aplica a la población de Manantlán que presenta la menor variación genética y se encuentra en la erupción más antigua del Mioceno-Pleistoceno (^{Martínez y Ramírez, 1998}). Esta observación podría confirmar que esta población está aislada de las demás poblaciones, tanto por distancia como por la orografía del lugar, lo que incide negativamente en el flujo genético y aumenta su endogamia.

Otra hipótesis aventurada, por confirmar, tiene que ver con los valores negativos y significativos de los estimadores de neutralidad que indican expansión poblacional de B. stormiae, lo que podría estar relacionado con los supuestos de la dinámica de la vegetación durante el Cuaternario Tardío, los cuales indican que los cambios en la estructura de la vegetación en el Holoceno fueron provocados por las actividades humanas durante los últimos 4000 años (^{Caballero-Rodríguez et al., 2018}). En ese sentido, se reporta que los primeros asentamientos humanos en Zacapu datan de aproximadamente 1500 a. C. (^{Reyes-Guzmán et al., 2018}). Se estima que, en su apogeo, la población prehispánica de esta localidad albergó alrededor de 20,000 habitantes, abarcando 1100 hectáreas sobre el pedregal (^{Michelet, 1998}; ^{Migeon, 1998}; ^{Smith, 2005}; ^{Forest, 2014}). Si suponemos que se utilizaron los recursos forestales disponibles en la región para la subsistencia de estos asentamientos humanos, esto pudo haber originado una recurrente perturbación al hábitat, incidiendo en la variación genética de esta especie. Sin embargo, el abandono posterior del pedregal durante la conquista pudo propiciar el crecimiento poblacional de B. stormiae; esto antes de nuevamente ser perturbada en tiempos actuales (^{Martínez-Salas et al., 2017}; ^{Reyes-Guzmán et al., 2018}).

Implicaciones para la conservación

Este es el primer estudio que contribuye al entendimiento de la genética poblacional de B. stormiae, especie considerada en peligro de extinción, distribuida en regiones aisladas. No obstante, es necesario llevar a cabo futuros análisis ampliando el tamaño del muestreo e incluyendo poblaciones de Colima, Guerrero, Oaxaca, Nayarit, Guatemala y El Salvador. Además, se sugiere explorar otros marcadores moleculares para fechar con precisión eventos de contracción y expansión demográfica para corroborar las hipótesis aquí expuestas y tener más elementos para establecer estrategias de conservación adecuadas para esta especie. Con base al trabajo de campo realizado en este estudio, considerando la perturbación observada en las diferentes localidades y el bajo número de individuos censados, se constata a esta especie en peligro de extinción (^{Martínez Salas et al., 2017}).