1. IntroducciÓnTop
Los Puertos de Tortosa-Beceite presentan una flora muy rica en endemismos, así como también albergan numerosas plantas que encuentran en este macizo el límite de su área de distribución. El Decreto 172/2008 de 26 de Agosto, de creación del Catálogo de Flora Amenazada de Cataluña, incluye 18 taxones que se dan en el ámbito del macizo de los Puertos y sierras vecinas.
La localización de poblaciones de especies raras, amenazadas o mal conocidas es un aspecto importante en la conservación de flora, tanto para evaluar su grado de amenaza como para desarrollar estrategias efectivas de conservación in situ y ex situ (Menon et al., 2010). Debido a que los datos de especies raras son generalmente escasos, es importante contar con formas eficientes para obtener datos adicionales. Los métodos de muestreo tradicionales resultan poco prácticos para la localización de especies raras, pues existe una gran porción del territorio con una probabilidad muy baja de ocurrencia. En este sentido, las predicciones espaciales realizadas a partir de modelos de nicho ecológico se pueden utilizar para el diseño de muestreos de campo. Dirigir las prospecciones hacia localidades con valores altos de idoneidad del hábitat es un método aconsejable para incrementar la eficiencia y reducir los esfuerzos y costes del muestreo (Guisan et al., 2006).
Habitualmente se dispone de datos precisos de localización de las especies amenazadas, aunque suelen ser escasos por la propia rareza de estas. La escasez de datos supone una limitación importante a la hora de modelizar la distribución de una especie. El poder predictivo de los modelos está claramente influenciado por el número de datos usados en su construcción (Carroll & Pearson, 1998; Cumming, 2000; Pearce & Ferrier, 2000; Kadmon et al., 2003; Hernández et al., 2006). Los niveles de incertidumbre asociados a los parámetros estimados decrecen generalmente con el incremento del tamaño de la muestra (Wisz et al., 2008). Varios autores han comprobado, no obstante, que si la pobreza de datos se debe a la rareza intrínseca de la especie y no a un muestreo sesgado, pocos registros pueden ser suficientes (Pearce & Ferrier, 2000; Kadmon et al., 2003; Hernández et al., 2006). Particularmente es así cuando se trata de modelizar el nicho de especies de hábitats muy específicos con escasas tolerancias ambientales y área de distribución restringida (Hernández et al., 2006; Williams et al., 2009).
Algunos de los que han estudiado los efectos del tamaño de la muestra (e.g. Pearson et al., 2007; Wisz et al., 2008; Williams et al., 2009), concluyen que el algoritmo de máxima entropía (Maxent) es uno de los que tiene una capacidad predictiva más alta con muestras pequeñas. Otros trabajos recientes (Gastón & García-Viñas, 2011) han demostrado, sin embargo, que el éxito de Maxent, sobre todo con muestras pequeñas, se debe a que usa regularización para evitar el sobreajuste. Estos autores afirman que pueden obtenerse resultados muy similares a los de Maxent con modelos lineales generalizados (GLM) que apliquen técnicas de regularización.
Por otro lado, para demostrar que un modelo predice de manera fiable es necesario proporcionar un elevado porcentaje de localidades de test, ya que en caso contrario el área predicha como adecuada es estadísticamente indistinguible de una predicción al azar (Anderson et al., 2002). El método más simple y común es partir los datos disponibles en muestras de entrenamiento y de test, utilizando particiones al azar o estratificadas espacialmente. Este método no resulta apropiado cuando los datos disponibles son limitados, ya que tanto los datos de entrenamiento del modelo como los de test resultan insuficientes (Anderson et al., 2002; Anderson & Martínez-Meyer, 2004). Cuando el tamaño de la muestra es muy pequeño el modelo es muy sensible a las localidades individuales, y un solo dato puede aportar mucha información al modelo (Pearson et al., 2007). Afortunadamente, existen métodos como la validación cruzada dejando uno fuera (Leave one out), que no implican la partición de los datos disponibles (e.g. Pearson et al., 2007).
En este trabajo se explora la utilidad de los modelos de nicho ecológico para la localización en campo de nuevas poblaciones de seis especies raras y amenazadas a escala global o regional, dentro del Parque Natural Els Ports. Los objetivos son (1) comprobar cómo afecta la limitación de datos y el tipo de hábitat ocupado por las especies en las predicciones realizadas por los modelos, (2) evaluar la utilidad de los modelos a la hora de diseñar un muestreo de campo de un conjunto de especies, y (3) incrementar el número de localidades conocidas dentro del Parque Natural de las especies de estudio.
2. Material y mÉtodosTop
2.1.Área de estudio
El ámbito de este estudio corresponde a la parte catalana del macizo de los Puertos de Tortosa-Beceite (NE Península Ibérica), si bien se usaron datos de presencia de una especie que ocurre en la vecina sierra de Cardó y los modelos se proyectaron también en esta zona (Figura 1). Las coordenadas más extremas son al noroeste 0º 8’ 58’’ y al sureste 0º 39′ 12″, el área total es de 1.138 km2, aunque el muestreo se llevó a cabo solo dentro del Parque Natural de Els Ports (PNP) de 354 km2.
|
Figure 1: Study área. |
Los Puertos de Tortosa-Beceite se encuentran situados entre las estribaciones más septentrionales del Sistema Ibérico y las sierras más meridionales del Sistema Costero catalán. El macizo de Cardó, separado de los Puertos por el río Ebro, forma también parte de la Cordillera Prelitoral meridional. Ambas sierras están formadas por materiales calcáreos mesozoicos que determinan un relieve muy abrupto. En la parte septentrional destaca el frente de conglomerados volcado hacia la depresión del Ebro. Los Puertos presentan una importante área de altura superior a los 1.000 m. La cima más alta es el Mont Caro, con 1.447 m. El macizo de Cardó culmina en la Creu de Santos (942 m).
El clima es bastante heterogéneo, ya que factores como la altitud y la orientación tienen una especial relevancia. La vertiente sureste de los Puertos y el macizo de Cardó, de influencia mediterránea, presentan un régimen de temperaturas más suave que la vertiente interior, que soporta unas condiciones más continentales. La temperatura media al pie de la sierra es de 14ºC, mientras que alrededor de los 1.000 m de altitud disminuye hasta los 11-12ºC. La distribución de las lluvias varía sustancialmente, mientras que la vertiente de levante recibe precipitaciones que oscilan entre los 700-900 mm/año, las de la vertiente continental son del orden de 450-600 mm/año (Royo et al., 2008). El viento es otro factor fundamental en la caracterización del clima del territorio, especialmente el viento del NO que domina con un 37% de recurrencia y alcanza velocidades de 50 a 110 km/h (Royo et al., 2008).
2.2. Especies de estudio
Se escogieron 6 especies consideradas amenazadas a escala nacional o regional y protegidas legalmente en Cataluña (Figura 2). Dos de ellas (Aquilegia paui Font Quer y Antirrhinum pertegasii Rothm.) son endemismos exclusivos del área de estudio. Otras dos (Erodium celtibericum Pau y Salix tarraconensis Pau) son endemismos ibéricos de área de distribución reducida que tienen en el área de los Puertos-Cardó una parte importante de sus poblaciones (Tabla 1). Finalmente, Atropa baetica Willk. y Prunus prostrata Labill. tienen distribuciones que sobrepasan los límites de la Península Ibérica pero que se consideran amenazadas en el ámbito estatal (solamente A. baetica) o regional (ambas).
|
Figure 2: Habitat of the target species: a) Aquilegia paui; b) Antirrhinum pertegasii; c) Erodium celtibericum; d) Prunus prostrata; e) Salix tarraconensis; f) Atropa baetica. |
Tabla 1: Forma vital, distribución general y porcentaje incluido dentro del PNP (para Salix se incluye también Cardó) y hábitat de las 6 especies de estudio.
Table 1: Life-form, general distribution and percentage included within the PNP (Cardó is also included for Salix) and habitat of the 6 target species.
| Especie | Forma vital* | Distribución | PNP % | Hábitat | Altitud (m) |
|---|---|---|---|---|---|
| Salix tarraconensis | NP | Montañas de Tarragona y N de Castellón | 33 | Pedregales fijados al pie de riscos o fisuras de rocas sombreadas | 450-1430 |
| Aquilegia paui | H | Puertos de Tortosa-Beceite (Tarragona) | 100 | Pedregales al pie de riscos orientados al N | 900-1400 |
| Antirrhinum pertegasii | Ch | Puertos de Tortosa-Beceite (Tarragona, Castellón y Teruel) | 70 | Fisuras de acantilados calcáreos, generalmente desplomadas | 300-1200 |
| Erodium celtibericum | Ch | Sistema Ibérico oriental (Tarragona, Teruel, Castellón, Cuenca y Guadalajara) | 15 | Crestas rocosas en sitios expuestos al viento | 850-1180 |
| Prunus prostrata | Ch | Disperso por el S de Europa, NW de África y SW de Asia | <1 | Grietas de roca en sitios expuestos al viento | 850-1330 |
| Atropa baetica | H | CE y SE de la Península y N Marruecos | <5 | Claros de bosque, matorrales abiertos, pedregales, sobre suelos ligeramente perturbados | 1200-1310 |
* Forma vital: Nanofanerófito (NP), Hemicriptófito (H), Caméfito (Ch).
Los datos de presencia para generar los modelos se extrajeron de la Cartografía de la flora Amenazada del PNP, en la que se localizan las poblaciones conocidas mediante polígonos o puntos con un error generalmente inferior a 100 metros. Para S. tarraconensis también se incluyeron algunas localidades de la sierra de Cardó, fuera del Parque. Todos los datos de presencia se convirtieron a puntos situados en el centro de una malla de 250 m de lado (resolución de trabajo, véase apartado 2.3). Solo se utilizaron los datos que tenían una precisión menor que la resolución de las variables. El número de celdas con datos de presencia que se obtuvo para cada especie (Tabla 2) varió entre 8 y 33.
Tabla 2: Catalogación en Cataluña (Decreto 172/2008) y en la Lista Roja Nacional (Bañares & al. 2008), datos de presencia originales (celdas de 250 × 250 m) y tamaño poblacional estimado en el PNP (se incluye Cardó para Salix).
Table 2: Classification in Catalonia (Decree 172/2008) and National Red List (Bañares & al. 2008), original presence data (250 × 250 grid cells) and estimated population size in the PNP (Cardó is included for Salix).
| Especie | Leg. CAT | Lista Roja España |
N. celdas 250 m |
UTM 1×1 |
Individuos estimados |
|---|---|---|---|---|---|
| Salix tarraconensis | VU | VU | 20 (11 PNP) | 13 (6 PNP) | <2500 |
| Aquilegia paui | EN | CR | 8 | 5 | <2500 |
| Antirrhinum pertegasii | VU | EN | 33 | 28 | >5000 |
| Erodium celtibericum | VU | – | 15 | 5 | >5000 |
| Prunus prostrata | VU | – | 16 | 8 | <2500 |
| Atropa baetica | VU | EN | 8 | 4 | 25 |
Fuentes: Royo et al., 2009; Sáez et al., 2009, 2010; Buira et al., 2009; Mayoral et al., 2013.
2.3. Resolución de trabajo y variables predictoras
La elección del tamaño de celda (resolución) de las variables predictoras es un factor importante que afecta a las predicciones (Guisan et al., 2007). Depende en gran medida del tamaño de celda de las capas ambientales disponibles (Graham et al., 2004), así como de la precisión de los datos de presencia de las especies a modelizar. El tamaño de celda condiciona también las prospecciones, la unidad de muestreo no debe resultar demasiado grande para poder ser explorada de forma eficiente en el campo (Guisan et al., 2006). Tanto la resolución de las capas climáticas como la precisión de los datos de presencia impedían en nuestro caso el uso de un tamaño de celda inferior a 200 m. Por ello, se optó por utilizar una resolución espacial de 250 m, que es proporcional a las divisiones de la cuadrícula UTM. Se considera un tamaño de celda suficientemente pequeño para representar condiciones de microhábitat, así como para ser prospectado íntegramente o casi en el campo.
Las variables predictoras seleccionadas deberían ser idealmente aquellas que se consideren las causantes directas de la distribución de la especie, aunque lo habitual es utilizar solo aquellas disponibles (Benito de Pando & Peñas de Giles, 2007). Por otra parte, las variables deberían afectar a la distribución de las especies a una escala adecuada a la extensión geográfica y el tamaño de celda empleada. Las variables climáticas son apropiadas a escala global o media; las variables topográficas afectan a la distribución de las especies a escala media o reducida, y las variables de cubierta de suelo influencian la distribución de las especies a escala reducida (Phillips et al., 2006). Teniendo en cuenta que el ámbito de este estudio es relativamente poco extenso y la resolución espacial pequeña, se consideraron las variables topográficas y de porcentaje de cubiertas del suelo como las más importantes a la hora de generar los modelos. No obstante, debido a que existen diferencias climáticas más o menos marcadas entre una vertiente y otra del macizo, se consideraron también algunas variables climáticas (Tabla 3).
Tabla 3: Variables predictoras.
Table 3: Predictive variables.
| Grupo | Variables | Fuente |
|---|---|---|
| Topográficas | Elevación (media, coeficiente de variación) Pendiente (máxima, media) Orientación (% plano, N, NE, E, SE, S, SW, W, NW) |
Modelo Digital de Elevaciones 30×30 m |
| Climáticas | Temperatura (min. T mes más frio, máx. T mes más cálido) Precipitación (media, coeficiente de variación anual) Radiación solar |
Atlas Climático Digital de Cataluña (ACDC) (Ninyerola et al., 2000) |
| Cubiertas de suelo (% superficie) |
Cobertura de arbolado Cobertura de matorral Cobertura de prados y herbazales Cobertura de roquedos |
3ª edición del Mapa de Cubiertas del Suelo de Cataluña (MCSC-3) (CREAF, 2005-2007) |
2.4. Análisis
2.4.1. Análisis de Componentes Principales
La colinealidad entre variables es un factor a tener en cuenta, ya que añade información redundante al algoritmo de modelización y puede generar resultados de difícil interpretación por sobreajuste del modelo (Benito de Pando & Peñas de Giles, 2007). La colinealidad entre algunas variables puede llegar a ser muy elevada en nuestro caso, teniendo en cuenta que las variables climáticas derivan de modelos de interpolación espacial creados a partir de los datos de las estaciones meteorológicas y las variables topográficas (Ninyerola et al., 2000).
Para reducir el número de variables (22) y de esta manera paliar los efectos de la colinealidad, se realizó un Análisis de Componentes Principales (PCA) (e.g. Pérez & Font, 2012). Se seleccionaron las 7 primeras componentes como variables predictoras para la construcción del modelo, que por sí solas explicaban más del 84% de la varianza total. Estas componentes principales o factores representan una combinación lineal de las variables originales, y tienen la ventaja de que son independientes entre sí. El uso de componentes, por el contrario, no permite conocer cuáles son las variables que más contribuyen a los modelos, pero teniendo en cuenta que nuestro objetivo era usar los modelos para dirigir el muestreo de campo, se consideró una cuestión poco relevante.
2.4.2. Obtención de los modelos binarios
El algoritmo que se utilizó para modelizar la distribución de las especies es el de máxima entropía (Maxent), el cual, a partir de la entrada de un conjunto de coordenadas de presencia y una serie de capas de información climática y geográfica, computa un modelo continuo de idoneidad del hábitat para la especie a lo largo del área de estudio. Más concretamente, el algoritmo calcula la distribución más uniforme de los puntos de ocurrencia, asumiendo que el valor esperado para cada variable predictora según esta distribución coincide con su media empírica, obtenida a partir de la media los puntos de presencia (Phillips, 2006).
Para ayudar a la validación e interpretación de los modelos, resulta conveniente generar mapas binarios fijando un umbral de corte por encima del cual los valores del modelo se predicen como “presencia” y por debajo como “ausencia” (Pearson et al., 2004). Multitud de propuestas diferentes se han utilizado para establecer umbrales de corte (Liu et al., 2005), pero cuando los modelos se generan sólo con datos de presencia, una opción es utilizar el umbral del valor mínimo en las presencias de entrenamiento (UMP) (e.g. Pearson et al., 2007; Guisan et al., 2006). Este umbral predice como presencia toda el área que queda por encima del valor más bajo obtenido en los datos de presencia utilizados para entrenar al modelo, y como ausencia el resto. Otro umbral más restrictivo es el que rechaza el 10% de las presencias de entrenamiento con los valores más bajos en el modelo (UP10). Se consideró a priori este umbral más oportuno para nuestro objetivo de reducir el esfuerzo en el muestreo de campo, ya que predice un área menor con valores de idoneidad más altos.
Para la obtención de los modelos se utilizó el software Maxent versión 3.2.19 (disponible en http://www.cs.princeton.edu/~schapire/maxent; Phillips et al., 2006). Se introdujeron los datos de presencia de las especies de estudio y las variables obtenidas del PCA como variables ambientales. Con el objetivo de evaluar la capacidad predictiva de los modelos (véase el apartado 2.4.3), se configuró la opción Crossvalidate, introduciendo tantas replicas como n datos de presencia tenía cada especie. El modelo final se obtiene a partir de la media de los n modelos generados. El formato de salida elegido fue el logístico, que representa un valor continuo entre 0 y 1 de idoneidad del hábitat para cada especie. Todos los mapas se transformaron en binarios aplicándoles el umbral de corte UP10 para el posterior diseño del muestreo.
2.4.3. Validación cruzada de los modelos
Dado que en nuestro caso el tamaño de las muestras era generalmente muy reducido (dos muestras de 8, tres inferiores a 21 y una de 33) se optó por utilizar el método de validación cruzada dejando uno fuera o Leave one out (Pearson et al., 2007), el cual no obliga a partir los datos disponibles. Este método consiste en retirar una localidad de muestra del conjunto de datos y construir el modelo utilizando las n-1 localidades restantes. Así, para una especie con n localidades observadas, se construyen n modelos separados para su evaluación (opción Crossvalidate de Maxent, véase el apartado 2.4.2). La capacidad predictiva del modelo se evalúa basándose en la habilidad de cada modelo para predecir la única localidad excluida del conjunto de datos de entrenamiento.
A partir de los resultados obtenidos para cada especie, se calculó la ratio de éxito como el cociente entre el número de veces que el modelo predijo la localidad excluida y el número total de datos de entrenamiento (n). Este cálculo se realizó tanto para los modelos binarios clasificados con el UP10 como con el UMP. Para evaluar si esta ratio era o no significativa se consideró el máximo de los P-valores calculados por el programa en cada test. Este P-valor es el resultado del test Chi cuadrado, en el que la hipótesis nula es que la tasa de omisión (1 o 0 en este caso) obtenida en el modelo con el área proporcional dada como presente (ambos valores dependen del umbral de corte fijado), no es diferente al de una predicción al azar con la misma área presente.
2.4.4. Diseño del muestreo de campo
Algunas de las especies de estudio comparten el hábitat; Salix tarraconensis, Aquilegia paui y Antirrhinum pertegasii ocupan siempre ambientes rupícolas sombreados y muy a menudo se encuentran próximas en el espacio. Por su parte, Prunus prostrata y Erodium celtibericum viven en crestas venteadas y también suelen ocurrir juntas en el espacio. Teniendo en cuenta lo anterior, y con el objetivo de reducir al máximo el esfuerzo, el muestreo de campo se dirigió hacia las zonas donde los modelos predijeron un mayor número de especies de estudio.
Para localizar las zonas con un mayor número de especies predichas se generó un mapa (Figura 3) a partir de la suma de los mapas binarios (con valores de 0 o 1) de las seis especies de estudio. Partiendo de la información de esta capa, se escogieron 8 áreas (A1-A8) para establecer itinerarios de muestreo. Los criterios para la elección de las zonas a prospectar fueron: agregación de celdas con especies predichas, presencia de celdas con 3 o más especies predichas, repartición entre la parte norte y sur, distribución más o menos equitativa del número de celdas predichas para cada especie (aproximadamente entre 15 y 30 celdas por especie) y accesibilidad. A pesar de todo, es importante tener en cuenta que no se puede despreciar cierto grado de subjetividad en la elección de las zonas y las celdas a prospectar debida a la experiencia del autor. El muestreo se restringió al ámbito del PNP, y se realizó durante los meses de Junio a Agosto de 2010, con una dedicación de unas 60 horas aproximadamente. Debido a la escabrosidad del terreno, en algunos casos se utilizaron prismáticos para localizar las plantas y técnicas de rápel para acceder a determinadas áreas.
|
Figure 3: Sum of the binary maps classified with UP10 and sampling areas: A1) Tossal de Montclí; A2) Rases del Maraco; A3) Tossa de la Reina; A4) Barcina; A5) Serrassoles; A6) Faixes Tancades; A7) Castell de l’Airosa; A8) Catinell. |
2.4.5. Evaluación de los modelos con los datos de campo
Para evaluar la capacidad predictiva de las clasificaciones binarias se calcularon los valores de sensibilidad y especificidad con los datos obtenidos en el muestreo de campo, y se representaron en el espacio ROC (acrónimo de Receiver Operating Characteristic). La sensibilidad o ratio de verdaderos positivos (RVP), mide la proporción de presencias nuevas encontradas en el campo que habían sido predichas por el modelo. Mientras que la especificidad o ratio de verdaderos negativos (RVN), mide la proporción de ausencias que habían sido predichas por el modelo como tal. Finalmente, para evaluar la eficacia en la localización de poblaciones nuevas en campo se calculó la precisión o valor predictivo positivo (PPV), que equivale a la proporción de celdas predichas por el modelo como presencia en las que se encontró la especie. Existen varias razones, que se discuten más adelante, por las que no debe tomarse este valor como una ratio de ocurrencia real respecto a lo predicho por el modelo.
RVP=VP/VP+FN
RVN=VN/VN+FP
PPV= VP/(VP+FP)
VP= Verdaderos positivos (celdas donde el modelo predijo presencia y se encontró la especie)
FP= Falsos positivos (celdas donde el modelo predijo presencia y no se encontró la especie)
VN= Verdaderos negativos (celdas donde el modelo predijo ausencia y no se encontró la especie)
FN= Falsos negativos (celdas donde el modelo predijo ausencia y se encontró la especie)
3. ResultadosTop
3.1. Umbrales de corte y validación cruzada de los modelos
Tal como se observa en la Tabla 4, los valores de UMP y UP10 coinciden en Aquilegia paui y en Atropa baetica, ya que al tener en ambos casos un tamaño de muestra inferior a 10, no se puede rechazar el 10% de los datos de entrenamiento con valores más bajos. Por lo general, se observa que los modelos que parten de un tamaño de muestra mayor, como es el caso de Antirrhinum pertegasii, predicen un mayor número de celdas. El rango de condiciones favorables para las especies suele ampliarse al incrementar el número de muestras, lo que generalmente conlleva un aumento del área predicha como idónea. Excepcionalmente, A. baetica presenta un área predicha relativamente grande a pesar de partir de un tamaño de muestra muy pequeño.
Tabla 4: Umbrales de corte, número de celdas predichas y ratios de éxito obtenidas en la validación cruzada.
Table 4: Thresholds, number of predicted cells and cross-validation success rates.
| Especies | n | UMP | N. celdas | Réxito | P-valor | UP10 | N. celdas | Réxito | P-valor |
|---|---|---|---|---|---|---|---|---|---|
| S. tarraconensis | 20 | 0,28 | 1.131 | 0,8 | 5,1e-02 | 0,56 | 313 | 0,7 | 1,5e-02* |
| A. paui | 8 | 0,31 | 92 | 0,75 | 5,2e-03* | 0,31 | 92 | 0,75 | 5,2e-03* |
| A. pertegasii | 33 | 0,03 | 2.833 | 0,97 | 1,4e-01 | 0,32 | 597 | 0,82 | 2,7e-02* |
| E. celtibericum | 15 | 0,22 | 379 | 0,87 | 1,9e-02* | 0,35 | 188 | 0,79 | 1,9e-02* |
| P. prostrata | 16 | 0,32 | 257 | 0,69 | 1,1e-02* | 0,42 | 191 | 0,69 | 7,8e-03* |
| A. baetica | 8 | 0,41 | 242 | 0,75 | 1,2e-02* | 0,41 | 242 | 0,75 | 1,2e-02* |
n: tamaño de la muestra; UMP: Umbral del valor Mínimo en las Presencias de entrenamiento; UP10: Umbral que rechaza el 10% de las Presencias de entrenamiento con los valores más bajos.
Las ratios de éxito obtenidas en la validación cruzada de los modelos (Tabla 4), tanto si se toma el UMP o el UP10, se hallan por encima de 0,69 en todas las especies. Esto quiere decir que al menos en un 69% de los casos la localidad excluida en los sucesivos n modelos cayó dentro del área predicha. Sin embargo, los modelos de Salix tarraconensis y A. pertegasii cuando se aplica el UMP predicen áreas muy grandes y se obtienen valores de P no significativos (> 0,05). Por lo tanto, en estos casos no se puede afirmar que los puntos de test caigan en las áreas predichas más a menudo que lo esperado por azar. Por el contrario, los valores de P obtenidos en los tests utilizando el umbral UP10 se encuentran muy por debajo de 0,05, y las ratios de éxito siguen siendo altas para todas las especies. Estos resultados apoyan, por tanto, la decisión de aplicar el UP10 en lugar del UMP a los modelos que se utilizan para dirigir el muestreo.
3.2. Resultados del muestreo de campo
La Figura 4 muestra los mapas binarios obtenidos al aplicar el UP10 a los modelos continuos. Solo se representa la parte correspondiente al PNP, aunque en el modelo de Salix tarraconensis también se incluyeron datos del macizo de Cardó. La Figura 5 muestra el incremento de datos de presencia en cuadrados de 250 m y 1 km de lado después del muestreo de campo. Se ha optado por representar también los resultados del muestreo en UTM de 1 × 1 km porque esta escala da una idea más real del incremento en términos de área de distribución conocida para la especie.
|
Figure 4: Binary maps classified with the UP10: a) Salix tarraconensis; b) Aquilegia paui; c) Antirrhinum pertegasii; d) Erodium celtibericum; e) Prunus prostrata; f) Atropa baetica. Original data (training data) and new data after the sampling in the PNP are shown. |
|
Figure 5: Total presence data before and after the survey by species, shown in 250 × 250 m and 1 × 1 km grid cells. |
Los valores más altos de especificidad, sensibilidad y precisión (Tabla 5) se han obtenido para el modelo de Aquilegia paui, por lo que es el que más se acercaría a una “clasificación perfecta” (Figura 6). El incremento de localidades conocidas (60% en UTM de 1 km) es también destacable, sobre todo si se tiene en cuenta que se trata de un endemismo con un área de distribución muy reducida. No obstante, todas las localidades nuevas consistieron en poblaciones con pocos individuos y cercanas a otras ya conocidas, a excepción de una localizada en A8 (Catinell) que se encuentra en una nueva UTM de 10 × 10 km (31TBF61) (Figura 4).
Tabla 5: Valores de sensibilidad, especificidad y precisión calculados a partir de los datos de presencia y ausencia obtenidos en el muestreo de campo.
Table 5: Sensitivity, specificity and precision values calculated on the basis of presence/absence data from field sampling.
| Especie | VP | VN | FP | FN | RVP | RVN | VPP |
|---|---|---|---|---|---|---|---|
| S. tarraconensis | 8 | 16 | 19 | 1 | 0,89 | 0,46 | 0,29 |
| A. paui | 5 | 30 | 8 | 0 | 1 | 0,79 | 0,38 |
| A. pertegasii | 6 | 25 | 10 | 0 | 1 | 0,71 | 0,37 |
| E. celtibericum | 5 | 27 | 15 | 1 | 0,83 | 0,64 | 0,25 |
| P. prostrata | 1 | 28 | 16 | 1 | 0,5 | 0,64 | 0,06 |
| A. baetica | 0 | 28 | 17 | 0 | 0 | 0,62 | 0 |
RVP: Ratio de Verdaderos Positivos - Sensibilidad; RVN: Ratio de Verdaderos Negativos - Especificidad; VPP: Valor Positivo de Predicción - Precisión.
|
Figure 6: Position of each binary classification in the ROC space (sensitivity vs (1-specificity)). The (0,1) point is a perfect classification, the diagonal line is the line of no-discrimination (random guess). |
El mayor número de presencias nuevas se ha obtenido para Salix tarraconensis, con un total de 9 localidades en cuadrados de 250 m (Figura 5), 8 de las cuales habían sido predichas por el modelo. Buena parte de las localidades nuevas se detectaron relativamente cercanas a otras ya conocidas, exceptuando una población localizada en A7 (Castell de l’Airosa), que se encuentra a más de 4,5 km de las hasta ahora conocidas. El valor de especificidad en este caso es relativamente bajo debido a que hay una gran cantidad de celdas donde el modelo predijo presencia y no se encontró la especie, lo que aleja su clasificación binaria del punto de “clasificación perfecta” (Figura 6).
Las prospecciones de Erodium celtibericum han dado como resultado valores razonablemente altos en la clasificación binaria. El número de localidades conocidas de Erodium después del muestreo se ha incrementado considerablemente (hasta un 60% en cuadrados de 1 km); con una población hallada a más de 4,5 km de las ya conocidas (A1, Tossal del Montclí), y las otras relativamente cercanas. Los valores asociados a la clasificación de Antirrhinum pertegasii son incluso mayores, pero al tratarse de un endemismo bastante frecuente en el macizo, su incremento en presencias resulta mucho más bajo, apenas de un 7% en UTM de 1 km.
En cuanto a Prunus prostrata, solo se localizaron dos presencias nuevas en el muestreo, de las cuales solo una había sido predicha por el modelo. De ahí que su modelo de clasificación binaria se sitúe cerca de la línea de clasificación aleatoria (Figura 6), es decir que según estos datos, no sería mucho mejor que una predicción hecha al azar. Finalmente, no se encontró ninguna localidad nueva de Atropa baetica en las prospecciones de campo. Los posibles motivos del fracaso en el muestreo de estas dos especies se discuten más adelante.
En síntesis se puede decir que el incremento de localidades de presencia ha sido notable. El número de UTM de 1 km conocidas dentro del PNP ha aumentado en al menos un 60% para A. paui, S. tarraconensis y E. celtibericum. Sin embargo, las poblaciones nuevamente reportadas estaban compuestas generalmente por un número bajo de efectivos, motivo por el cual hayan podido pasar inadvertidas hasta el momento. El mayor número de localidades nuevas se ha encontrado en A6 (Faixes Tancades), que es una de las zonas de mayor interés botánico del PNP, y además, ha sido relativamente poco explorada. El número de falsos negativos ha sido muy bajo en todas las especies (como máximo 1), por lo que se obtienen valores altos de sensibilidad en todos los casos excepto para A. baetica y P. prostrata. Asimismo, los valores de precisión para las cuatro primeras especies (Tabla 5), aunque no sobrepasan de 0,4 en ningún caso (el mayor es el de A. paui de 0,38), se pueden considerar altos al tratarse de especies raras (véase el apartado 4.3).
4. DiscusiÓnTop
4.1. Algoritmo y umbral de corte
Los resultados obtenidos en la validación cruzada (Leave one out) demuestran que los modelos de las seis especies tienen una capacidad predictiva razonablemente alta. Por otra parte, el muestreo de campo, dirigido hacia las áreas con un mayor número de presencias predichas, ha dado como resultado la localización de 28 nuevas presencias de cinco especies, el 89% de las cuales habían sido predichas por los modelos. En vista de estos resultados se considera adecuada la elección de Maxent dado que, en general, se han obtenido modelos aceptables con muestras pequeñas. Cabe mencionar, no obstante, que los autores que han comparado Maxent con modelos lineales generalizados (GLM) con regularización para muestras pequeñas, han obtenido resultados muy similares (e.g. Gastón & García-Viñas, 2011). Por tanto, muy probablemente el uso de GLM regularizados en nuestro caso hubiera dado resultados igualmente satisfactorios.
En cuanto al umbral de corte, nunca debe elegirse arbitrariamente sino que depende de la finalidad de los modelos (Wilson et al., 2005). Si la prioridad es reducir el error de comisión (cuando el modelo predice presencia y en realidad es ausencia), resulta conveniente utilizar un umbral de corte más alto a expensas de obtener un error de omisión mayor (Hernandez et al., 2006). Teniendo en cuenta que nuestro objetivo en el muestreo era encontrar presencias nuevas con el menor esfuerzo posible, se consideró más indicado aplicar un umbral más alto, en este caso el UP10, con el fin de obtener un área predicha menor con valores de idoneidad más altos. En cualquier caso, los datos obtenidos en el muestreo de campo demuestran que, aun aplicando un umbral restrictivo como el UP10, los errores de comisión son muy superiores a los de omisión (hay muchos falsos positivos y apenas falsos negativos).
Los resultados obtenidos en la validación cruzada consolidan también la elección del UP10 frente al UMP. Las ratios de éxito aplicando el UP10 han sido razonablemente altos (≥0,69) y significativos para todas las especies. Además, los modelos con un área predicha muy grande son poco realistas cuando se trata de especies raras.
4.2. Área de distribución potencial versus área de distribución real
Tanto los resultados del muestreo de este estudio como algunas revisiones corológicas (e.g. Aparicio & Balada, 2005; Sáez et al., 2009), parecen demostrar que Salix tarraconensis, Aquilegia paui y Antirrhinum pertegasii faltan en la mitad septentrional del macizo. Uno de los factores que podría limitar su presencia es la marcada continentalidad que afecta a la parte norte de los Puertos. La mayor distancia al mar y la presencia del macizo de Cardó, que actúa como una barrera impidiendo la entrada de nieblas húmedas procedentes del mar, determinan unas condiciones de mayor sequedad y temperaturas más extremas en esta zona. Sin embargo, los modelos de estas especies, y especialmente el de Salix, proporcionan valores altos de idoneidad del hábitat en algunas localidades rocosas de la mitad norte, posiblemente debido a una contribución débil de la componente climática al modelo. A ello se deben en gran parte la elevada cantidad de falsos positivos obtenidos en las prospecciones de campo de Salix, que se han dado mayormente en el sector norte del PNP. En el caso de Antirrhinum, los valores idóneos se concentran principalmente en la vertiente marítima, que es donde la especie es más frecuente en realidad. A pesar de que existe correspondencia con el área real, y de que los resultados en la prospección han sido satisfactorios, se considera que el área predicha en el modelo binario (597 celdas) es desproporcionadamente grande al tratarse de una especie relativamente poco abundante.
El modelo de Aquilegia paui proporciona un área idónea muy reducida, que se ciñe básicamente a algunos puntos culminantes de la vertiente marítima (Figura 4). Esto concuerda con la especificidad del hábitat y la rareza de este endemismo estricto del macizo. Se había considerado una planta extremadamente localizada y escasa (30 individuos en Sáez & Guàrdia, 2003; 339 en Martinell et al., 2007), hasta que en las prospecciones realizadas por el personal técnico del Parque se localizó en una nueva UTM de 10 km y se advirtió que era bastante más abundante de lo que se pensaba (Buira et al., 2009; Sáez et al., 2010). Esto puede atribuirse a la inaccesibilidad de la mayor parte de las canales rocosas donde habita. A pesar de todo, se han obtenido resultados buenos en el muestreo de campo de este trabajo, con un incremento considerable de localidades conocidas y una eficacia razonablemente alta (precisión de 0,38). Esto sugiere que su modelo tiene una capacidad predictiva alta, y por otra parte, que en las pocas áreas donde se reúnen las condiciones favorables la especie suele ocurrir.
La distribución conocida de Erodium celtibericum en los Puertos es prácticamente inversa a la de las especies anteriores, aparece solo en algunas crestas ventosas de las montañas de la mitad septentrional del macizo. En este caso, la distribución potencial arrojada por el modelo se corresponde bastante con la real, las celdas de mayor idoneidad del hábitat se concentran en zonas culminantes de la zona norte (Figura 4). Por otra parte, el área predicha en la clasificación binaria es relativamente reducida, lo que concuerda con la rareza de la especie en el macizo. La distribución real de Prunus prostrata coincide en gran parte con la de E. celtibericum y, en consecuencia, sus modelos son muy similares. Sin embargo, los resultados del muestreo para esta especie sugieren que podría ser mucho más escasa o bien que su modelo es poco fiable. La mayor parte de zonas visitadas ofrecían aparentemente un hábitat favorable para la especie, por lo que el fracaso podría atribuirse, en parte, a la escasa detectabilidad de las plantas (sin flores y en algunos casos sin hojas) durante el periodo de prospección.
Atropa baetica posee con diferencia el área de ocupación más reducida del grupo de plantas estudiado. Sólo se conocen 5 poblacionales (20-30 individuos en total) distribuidas por la zona culminante de la Barcina-Caro (A4), en un área inferior a 100 ha (Mayoral et al., 2013). A pesar de que su modelo partió de un conjunto de datos de presencia muy escasos (8) y agregados, proyectó un área relativamente grande de condiciones potencialmente idóneas en las zonas más elevadas de la mitad meridional del PNP (Figura 4). Este resultado no parece del todo desacertado si se tiene en cuenta que el hábitat de la especie en el macizo es poco específico. Los pocos individuos conocidos viven en claros de pinar, matorrales abiertos e incluso pedregales, ocupando una zona quemada hace décadas, en diferentes exposiciones y pendientes. Según Martínez et al., (2004) es una planta extremadamente rara en la Península, que parece ir asociada a comunidades inestables, integradas por una mezcla de especies de bosque y de sus etapas de sustitución. Todo esto sugiere que la presencia de A. baetica puede estar más ligada a determinadas interacciones bióticas que a unas condiciones ambientales muy específicas, y que por tanto, los modelos de nicho ecológico aplicados con el objetivo de este trabajo tienen una utilidad escasa.
4.3. Utilidad de los modelos
Los resultados obtenidos en el muestreo de campo han permitido incrementar considerablemente el número de localidades conocidas de algunas especies. En el caso de Aquilegia paui, Salix tarraconensis y Erodium celtibericum, el incremento ha sido de al menos un 60% en UTM 1 × 1 km. De A. paui, catalogada como en peligro, se conocen actualmente unas 12 poblaciones con alrededor de 2.500 individuos. Cabe decir, que las poblaciones nuevamente reportadas aquí constaban, por lo general, de un número bajo de efectivos, por lo que el incremento en abundancia conocida no sería tan significativo como el de cuadrículas de presencia. Tampoco se han ampliado las áreas de distribución conocida a gran escala (solo A. paui se ha citado en una UTM de 10 × 10 km nueva), pero se ha afinado en el conocimiento de las áreas de ocupación real. La delimitación precisa de las poblaciones de plantas amenazadas permite gestionar los usos que se dan en el Parque Natural y puedan afectarlas, realizar un seguimiento temporal de su extensión, elaborar modelos de selección de hábitat más precisos, evaluar el estado de conservación mediante algunos criterios UICN e incluso estimar tamaños poblacionales en algunos casos (Goñi et al., 2006).
La eficacia en la localización de las tres especies arriba mencionadas y de Antirrhinum pertegasii ha sido considerablemente alta tratándose de plantas en general muy escasas, puesto que los valores de precisión han variado entre 0,38 y 0,25. No obstante, existen varias razones por las que no debe tomarse este valor como una ratio de ocurrencia real respecto a lo predicho por el modelo. La experiencia del autor en la elección previa de las celdas a muestrear y la disposición agregada y con múltiples especies predichas de estas (en comparación con las celdas aisladas), podría haber incrementado la ratio de éxito en la prospección de estas cuatro especies. Por el contrario, la baja detectabilidad de Prunus prostrata durante el periodo de muestreo podría haber determinado el fracaso en las prospección de esta especie (precisión próxima a 0). La probabilidad de detección es un factor a tener en cuenta, ya que a menudo varía con las mismas covariables que determinan la probabilidad de ocurrencia (Tyre et al., 2003; Yackulic et al., 2013).
La capacidad predictiva de los modelos, evaluada tanto a través de la validación cruzada como de los resultados del muestreo de campo, ha sido más o menos satisfactoria según la especie. El modelo de Aquilegia es el que tendría una capacidad predictiva más alta, a pesar de partir de un tamaño de muestra de tan solo 8 presencias. Estos resultados apoyan a otros estudios (e.g. Raxworthy et al., 2003; Engler et al., 2004; Guisan et al., 2006; Hernández et al., 2006), que ponen de manifiesto que las predicciones de los modelos pueden resultar útiles en la priorización de los esfuerzos de muestreo de especies amenazadas con pocos datos de presencia. El éxito en el uso de los modelos se debe en gran parte a la elevada especificidad del hábitat de las especies utilizadas en nuestro caso, todas excepto Atropa, viven en ambientes rupícolas con condiciones ambientales muy particulares. Generalmente, los modelos de especies con área de distribución reducida y especialistas del hábitat son más certeros que aquellos de especies de amplia distribución y de rango de tolerancia ecológica amplio (Kadmon et al., 2003; Thuiller et al., 2004; Elith et al., 2006; Hernández et al., 2006; Williams et al., 2009).
Los resultados de este estudio muestran también que otro método ventajoso para optimizar el muestreo de campo es llevarlo a cabo para un conjunto de especies a la vez. Dentro del grupo de las plantas endémicas y amenazadas, por ejemplo, son abundantes las especies casmofíticas y glareícolas (Domínguez et al., 2003; Goñi et al., 2006). En este sentido, los modelos pueden utilizarse para dirigir las prospecciones de grupos de plantas amenazadas con hábitats y áreas de distribución similares, con el objetivo de incrementar la probabilidad de hallazgo.