Les Sentinelles du Climat

Les observations sont reportées sur fond cartographique. La carte de chaleur ou carte de densité est un outil sous QGIS 2.18.14 qui simule la densité d’individus observés dans un périmètre donné. Plus le nombre d’observations est important, plus la couleur est intense. Cet outil permet ainsi de visualiser les zones de concentration au sein de chaque « population ». Les cartes de chaleur sont établies sur l’ensemble de chaque site, selon la somme des observations annuelles de chaque espèce. Les altitudes et coordonnées moyennes par espèce sont également calculées.

L’estimation d’abondance par la méthode N-mélange requière des suivis répliqués dans le temps et dans l’espace (Royle, 2004). Pour construire un modèle de mélange N, il est nécessaire de définir K, un paramètre de réglage sans sens biologique et de choisir la loi de distribution de densité la plus appropriée . Lorsque K augmente, l’abondance de l’estimateur augmente aussi jusqu’à sa stabilisation. Nous avons donc choisi la valeur de K au-delà de laquelle l’abondance se stabilisait. Parallèlement, nous avons défini la meilleure loi de distribution de densité pour les données de comptage étudiées. Dans le package unmarked utilisé , trois lois de densité ont été mises en œuvre, la loi de Poisson (P), la loi de Poisson avec inflation de zéros (ZIP) et la loi binomiale négative (NB). La distribution ZIP gère la quantité élevée de zéro dans les données de comptage (Tu, 2006) et la distribution NB gère la sur-dispersion des données . En utilisant le modèle constant, p(.)lambda(.), nous avons défini la distribution la plus appropriée entre P, ZIP et NB à partir du critère corrigé d’Akaike (AICc). Ensuite, nous avons calibré la valeur K et réalisé des tests d’ajustement avec un test paramétrique bootstrap avec 1000 simulations. Ensuite, le modèle variable p(t)lambda(.) et les autres prenant en compte des covariables standardisées ont été construits. Tous les modèles ont été comparés selon les critères de l’AICc. L’ajustement des modèles avec un delta AICc inférieur à 10 a été testé. Le meilleur des meilleurs modèles a été sélectionné pour estimer l’abondance sur chaque placette.

Cette année, nous avons réalisé cette analyse pour les deux sites suivis sur placettes, Mendive et Esterençuby. Les covariables de site sont le pourcentage d’occupation du sol par les roches et l’altitude ; les covariables d’échantillonnage sont la température, la nébulosité et la visite, comme en 2019. Ces covariables n’ont été relevées qu’en 2019 et 2020, les modèles réalisés sur les données de 2017 et 2018 se restreignaient aux modèles simples p(.)lambda(.) et p(visite)lambda(.).

Lors de cette modélisation, les indications du livre de Guisan et al. Guisan , en particulier celles du chapitre 19.1 sont respectées scrupuleusement. La modélisation corrélative se réalise selon plusieurs étapes, présentées brièvement ci-dessous, détaillées dans la méthodologie décrite dans la thèse de Florèn HUGON et dans le chapitre du rapport précédent .

La période de référence est 1991-2020 et les projections futures sont étudiées sur trois horizons, proche 2021-2050 (horizon 1), moyen 2041-2070 (horizon 2) et lointain 2071-2100 (horizon 3) pour trois scénarios RCP, 2.6, 4.5 et 8.5. Pour décrire ces périodes futures, nous avons utilisé le terme « couple scénario-horizon », nous distinguons ainsi les couples RCP26h1, RCP26h2, RCP26h3, RCP45h1, RCP45h2, RCP45h3, RCP85h1, RCP85h2 et RCP85h3.

Les points d’occurrence résultent du suivi de l’abondance des deux lézards par le programme et de la base de données GBIF. Pour éviter la pseudo-réplication des données, un seul point d’observation est sélectionné dans chaque maille .

Les variables climatiques issues des simulations ALADIN52 (CNRM 2014) sont utilisées pour le calcul d’indices climatiques cohérents avec la biologie et le cycle de vie des espèces . Ainsi, à partir de neuf variables – températures minimale et maximale journalière à 2 m en °C, précipitations liquides en mm, chutes de neige en mm, rayonnements infra-rouge et visible incident à la surface enW.m-2, humidité spécifique en g.kg-1, vitesse du vent à 10 mètres en m.s-1, vitesse maximum du vent à 10 m en m.s-1 – 61 indices climatiques dont 59 mensuels et 2 annuels ont été calculés.

Nous avons sélectionné des indices climatiques qui nous semblaient pertinents par rapport au cycle de vie des espèces. « Le cycle de vie des lézards que nous étudions est scindé en deux périodes, une période d’activité et une période d’hivernation définies par les conditions environnementales . Chez le Lézard de Bonnal, la période d’activité débute dès la fonte des neiges et prend fin dès les premières chutes de neiges . Cette période varie selon l’altitude du site d’étude, nous avons défini la période d’activité du 01 mai au 31 octobre et la période d’hivernation du 01 octobre au 31 mai. Pour le Lézard catalan, la période d’activité est bornée par la fin de l’hiver et le milieu de l’automne , nous avons définis la période d’activité du 01 mars au 31 octobre et la période d’hivernation du 01 novembre au 31 mars. Nous n’avons pas distingué la période de reproduction et la période d’activité car la période post-reproduction permet le renouvellement des réserves énergétiques pour la période de reproduction de l’année suivante . De plus, chez le Lézard de Bonnal, il a été montré que l’activité n’est pas différente selon la période de reproduction et post-reproduction . Ainsi, pour les deux espèces, 29 indices climatiques ont été sélectionnés.

Les variables environnementales ont été choisies selon l’écologie de l’espèce. Le Lézard de Bonnal est présent sur des habitats soumis aux vents avec de forts contrastes thermiques et hydriques. La pluviométrie est élevée, la neige est présente sur 6 à 8 mois de l’année . La température de l’air et le rayonnement solaire expliqueraient en partie son activité . Le Lézard catalan est soumis à des conditions climatiques similaires . Ainsi, nous avons sélectionné pour ces espèces, les variables « altitude », « exposition », « pente », « ombrage » issues d’un MNT 100 mètres puis calculées pour la résolution 1km, la variable « pourcentage de fragments grossiers dans le sol » , les variables « zone humide » , « cours d’eau » issue de la BD CARTHAGE, « plans d’eau » issue de la BD TOPO et enfin les variables occupation du sol, « prairies », « forêts » et « pelouses », calculées par le CESBIO en pourcentage.

L’ensemble des 34 variables a été obtenu à la résolution 1km², ce qui représente un compromis entre les données climatiques et les données d’habitat et projeté dans le système géodésique mondial WGS84. Les tests de corrélation ont été réalisés en deux à deux en ne conservant que les variables ayant un coefficient inférieur à 0,7.

Pour modéliser la répartition, nous avons utilisé le package BIOMOD2  qui est une version plus aboutie de BIOMOD . Ce package, facile d’utilisation, permet de sélectionner les variables selon leur importance, i.e. leur pouvoir explicatif, et les méthodes selon les courbes de réponses des variables, puis de projeter les répartitions au présent et au futur selon différents scénarios .

Les méthodes de régression (GAM, GLM, MARS), de classification (ANN, FDA) et de ré-échantillonnage (CTA, GBM, RF) ont été utilisées. Toutes ces méthodes requièrent un jeu de pseudo-absence. À l’issue d’une étude bibliographique , nous avons décidé de générer 10 fois 1000 pseudo-absences afin de réduire l’effet du jeu de pseudo-absences .

Pour évaluer la qualité des modèles, il est nécessaire d’utiliser soit un autre jeu de données indépendant, soit la validation croisée, soit un jeu de données indépendant issu de notre jeu de données en scindant celui-ci en deux parties . Ne disposant pas d’un autre jeu de données, il a été décidé de scinder le jeu de données en deux, avec 70% pour l’apprentissage et 30% pour l’évaluation selon 10 runs. Les modèles ont été évalués au moyen du critère TSS – True Skill Statistic , un modèle étant considéré comme bon si ce critère est supérieur à 0,7 . Enfin, pour déterminer l’importance des variables, 10 permutations ont été réalisées.

Une première modélisation a été réalisée à partir des 8 méthodes et des variables sélectionnées selon le critère de corrélation maximal de 0,7. Les sorties obtenues ont permis d’effectuer une sous-sélection des variables et des méthodes selon les coefficients d’importance des variables et les courbes de réponses. Les méthodes dont les courbes de réponses étaient lisses ou en escalier ont été supprimées. Les 10 premières variables par ordre d’importance pour les méthodes retenues ont été sélectionnées.

Les différentes méthodes de modélisation sont régulièrement comparées entre elles . Les articles scientifiques comparatifs indiquent l’importance de tester plusieurs fois le modèle et ses erreurs de projection . Il est également souvent proposé de combiner plusieurs modèles entre eux pour obtenir un modèle d’ensemble et diminuer ainsi les incertitudes de chacun . Les modèles d’ensemble peuvent être obtenus en calculant la moyenne des probabilités des modèles, la médiane, la moyenne pondérée et par « commitee averaging » (transformation des probabilités de présence en système binaire 0/1 puis moyenne) . La méthode coefficient de variation des probabilités permet d’estimer l’incertitude, la variation entre tous les modèles pour la modélisation d’ensemble, une valeur élevée indique que l’incertitude est importante ; cette méthode n’est pas bornée dans l’intervalle [0,1] comme le sont les autres . Ici, les 5 méthodes ont donc été réalisées avec les modèles ayant un critère TSS au-dessus de 0,7.

À partir du modèle d’ensemble, les projections au temps présent et aux horizons futurs ont été obtenues. Enfin, les pourcentages de changement d’aire de répartition entre la projection au temps présent et les projections futures ont été calculés, c’est à dire les pourcentages de perte et de gain d’habitat ainsi que les pourcentages d’habitat qui restaient inoccupés ou occupés.

Les premières étapes de la modélisation mécanistique consistent à définir la méthode de calcul du temps d’activité / d’inactivité au cours de la saison de reproduction et définir le lien entre temps d’activité / d’inactivité et probabilité de persistance. Ensuite, le temps d’activité est lié à une variable environnementale, communément le maximum de la température de l’air . Cette équation permet de projeter le temps d’activité selon différents scénarios climatiques. À partir du temps d’activité, la probabilité de persistance peut ensuite être calculée pour obtenir une carte de probabilité de présence.

En premier lieu, la représentation des séries de températures opérantes et le calcul de statistiques descriptives permet l’observation des variabilités spatiales et temporelles. Les statistiques calculées sont l’étendue, les quartiles, les quantiles 0,025 et 0,975, l’écart type et la moyenne.

Pour le calcul du temps d’activité, nous utilisons la méthode « low-constraint » où le temps d’activité est défini lorsque la température corporelle est incluse dans la fenêtre thermique d’activité bornée par VTmin et VTmax .

Pour chaque journée de la période de reproduction, les heures d’activité journalières (HaDaily) sont obtenues par somme pour obtenir ensuite le nombre total d’heures d’activité (HaTot) (voir équation ci-dessous) . La méthode « organismal » permettant de considérer l’activité comme une grandeur continue est en cours de développement.

L’indice de persistance (PI) a été défini à partir d’un seuil de persistance (Tpersist) qui représente le temps d’activité total minimum nécessaire au recrutement sur un site au cours d’une période de reproduction. Tpersist est calculé en multipliant la durée d’activité quotidienne de référence (HaDailyref), c’est-à-dire le nombre moyen d’heures d’activité quotidienne, par le nombre de jours de la période de reproduction. PI est le rapport entre HaTot et Tpersist, il fournit une estimation ponctuelle, comparable au sein d’un site, entre les sites et entre les années, ce qui permet l’étude de la variabilité. S’il est compris dans l’intervalle [0,1], il représente une probabilité d’extinction (extP) alors que s’il est supérieur à 1, il représente une marge de sécurité de temps d’activité (HaSM), par rapport à l’augmentation prochaine de la température moyenne et de sa variabilité  :

Le seuil de persistance peut être défini sur la base d’observations naturalistes énoncées dans la littérature. Pour le Lézard de Bonnal, d’après les observations de terrain effectuées sur les populations proches de nos sites d’étude, HaDailyref serait environ égal à 4,5 heures. Arribas indique un pic d’activité de 10h30 à 12h et une baisse importante de 13h à 18h. L’étude de Pérez-Mellado  indique également un pic matinal puis un second pic en fin d’après-midi. Enfin, Pottier indique que les individus sont actifs entre 8h et 12h30 seulement avec un pic également de 9h à 11h. Ainsi, nous avons défini Tpersist = 225 heures (4,5 heures multipliées par 50, le nombre de jours de la période de reproduction).

Pour le Lézard catalan, Martin-Vallejo et al. montrent que le nombre de captures est plus important le matin de 8h à 12h puis aux alentours de 18h. Carneiro  souligne également que le pic d’activité du matin est plus long que celui de l’après-midi. D’après Matthieu BERRONEAU, le naturaliste qui suit les Lézards catalan et de Bonnal au cours du programme, le temps d’activité journalier du Lézard catalan devrait être plus long que celui du Lézard de Bonnal. Enfin, le développement d’un modèle de rythme d’activité journalier par Ferri-Yanez montrait un temps journalier de 8h34min où l’individu atteignait VTmin et de 5h45 où l’individu atteignait la température optimale Topt comprise entre VTmin et VTmax. Ces temps incluent tous deux les moments où les individus sont en dessous de VTmin et où l’activité n’est pas possible. Compte tenu de l’ensemble de ces éléments, nous pouvons faire l’hypothèse que HaDailyref serait égal à 5 heures et déduire Tpersist = 460 heures (5 heures multipliées par 92, le nombre de jours de la période de reproduction).

Ce travail a été réalisé cette année pour le Lézard de Bonnal  et sera prochainement réalisé pour le Lézard catalan. Les séries chronologiques HaDaily ont été comparées par paires pour étudier les différents types de variabilité. La variabilité spatiale a été étudiée au niveau intra-site (4 tests), inter-exposition (6 tests), inter-altitude (6 tests) et inter-site (8 tests). La variabilité temporelle a été étudiée uniquement au niveau inter-annuel (18 tests). Pour comparer les séries, nous avons utilisé les tests bilatéraux de Wilcoxon de rangs signés pour échantillons appariés avec la fonction wsrTest du paquet asht . Les p-values estimées étaient associées à l’estimateur de Hodges-Lehmann et à son intervalle de confiance à 95 %, il correspond à la pseudo-médiane de la série de différences . Toutes les p-values ont été corrigées par la méthode FDR pour identifier des comparaisons réellement significatives . Enfin, nous avons tracé les valeurs p transformées en base décimale logarithmique et les avons comparées au seuil de 1,30 correspondant à -log10(0,05). Les indices de persistance ont été comparés en intra-site, inter-pente, inter-altitude, inter-site et inter-année. Nous avons également calculé la moyenne et l’écart-type interannuel de l’indice de persistance pour quantifier le risque d’erreur dans les estimations de la persistance si les études sont menées sur une seule année.

L’intégration de la notion continue de l’activité est fastidieuse car elle requière la définition d’une probabilité d’activité en fonction de la température corporelle. Ce travail a été débuté en 2019 en utilisant de manière novatrice les résultats de la modélisation de l’abondance . Il est attendu une courbe en cloche avec un maximum d’activité sur la gamme de température préférée. L’activité serait possible dans la fenêtre d’activité thermique définie par VTmin et VTmax , VTmin étant la température qui signale à l’organisme qu’il faut rechercher des températures plus chaudes ou aller en refuge et VTmax celle qui permet d’éviter le comportement de surchauffe . Les estimations d’abondance répondent ainsi à deux questions, y a-t-il une variation des effectifs au cours du temps et quelle est la relation entre le taux d’activité et la température ?

Pour calculer les proxys d’activité, nous devons d’abord estimer la probabilité de détection pour chaque couple visite – placette (pij) et l’abondance pour chaque placette (Ni). Pour calculer le proxy sur la placette i à la visite j (pAij), nous avons défini l’équation suivante : pAij = (nij / pij) / Ni. Cette équation prend en compte le nombre d’individus disponibles pour la détection et la probabilité de détection. Par exemple, un comptage de dix individus avec une abondance estimée à cinquante et une probabilité de détection estimée à 0,25 donne un proxy d’activité de 0,8. Si nous observons le même nombre d’individus avec la même probabilité de détection et une abondance de cent individus, le proxy d’activité ne sera que de 0,4. De même, si nous comptons toujours le même nombre d’individus avec une probabilité de détection de 0,5 sur un site où l’on trouve cent individus, l’indicateur de l’activité ne sera que de 0,2. Ceci illustre la pertinence de quantifier l’abondance réelle sur chaque placette d’étude et la probabilité de détection pour chaque visite sur chaque placette. Pour obtenir les proxys de l’activité en termes de probabilités dans l’intervalle [0,1], notés pA(ij scaled, les pAij ont été divisés par la valeur la plus élevée.

Les proxys d’activité ont ensuite été associés à la température corporelle Tbij plutôt qu’à la température de l’air mesurée sur chaque placette à chaque visite. Pour cela, nous avons établi une relation linéaire Te=f(Tair) à partir des données de températures opérantes et de l’air mesurées sur Mendive à 730 mètres (modèle biomimétique) et 727 mètres (station météo), sur la période du 1er mai au 31 août 2018. À l’aide de cette relation, nous avons calculé les températures corporelles dépendantes des températures de l’air mesurée au moment de l’échantillonnage. Puis, nous avons couplé les proxy d’activité pA(ij scaled) avec la température corporelle Tbij pour chaque placette et visite. Nous avons modélisé cette courbe avec une fonction polynomiale du troisième degré. Nous n’avons utilisé que la première et la dernière valeur zéro car la pertinence des autres zéros pouvait être discutée (le zéro ne serait pas lié à l’activité, l’espèce serait réellement absente de la placette).

Cette méthodologie, dite « organismal » a été réalisée à partir des données de 2019 du site de Mendive, elle sera réitérée très prochainement à partir des données 2020. Cependant, elle est encore en développement et ne peut donc pas être utilisée pour la modélisation mécanistique dans notre étude. Une méthodologie plus simple à mettre en œuvre a été retenue (présentée dans la partie précédente), c’est la méthode dite « low constraint » où le temps d’activité est défini lorsque la température corporelle est incluse dans la fenêtre thermique d’activité bornée par VTmin et VTmax.

Pour le Lézard de Bonnal, aucun protocole au sein du programme ne permet de déterminer ces grandeurs. D’autres études ayant été menées sur des populations spatialement proches des nôtres, nous avons privilégié l’utilisation de leurs résultats en faisant l’hypothèse « raisonnable » que nos populations d’étude présentaient les mêmes caractéristiques thermiques que la population la plus proche spatialement étudiée dans la bibliographie . Nous avons défini VTmin et VTmax en utilisant l’étude d’Arribas menée sur le Lac Bleu du massif de Bigorre. Cette étude était le meilleur compromis entre la distance géographique avec nos sites d’étude et la disponibilité des données mesurées sur le terrain. Ainsi, VTmin = 20,8°C et VTmax = 35,2°C.

Pour le Lézard catalan, les températures préférées de 7 individus capturés sur le site du chemin de la Mature (altitude 850 mètres) ont été mesurées toutes les minutes pendant une heure en 2016 par Barry Sinervo et Frank D’Amico. Le calcul des quantiles 5ème et 95ème permet d’obtenir VTmin = 26,1°C et VTmax = 38,4°C. Ces résultats sont en accord avec d’autres études, Arnold (1987) avait déterminé que les températures corporelles d’individus actifs variaient de 25,4 à 38,4°C pour des populations vivant au nord de Madrid à 1400 mètres d’altitude et Carretero et al. avaient indiqué que les températures préférées variaient de 26 à 37.3°C pour des populations proches de Barcelone à 150 mètres d’altitude.

Dans le contexte sanitaire Covid 19 particulier de l’année 2020, les modèles biomimétiques n’ont pas pu être déployés sur les sites du Lézard catalan. Nous avons profité de leur disponibilité pour étudier les variations des températures opérantes (Te) dans les pierriers utilisés par le Lézard de Bonnal, là où se réfugient les individus lorsque les températures de l’air sont trop faibles ou trop élevées pour leur permettre de maintenir leur température interne dans la fenêtre de tolérance thermique, bornée par CTmin = 6,10°C et CTmax = 42,20°C. La fenêtre thermique d’activité, comme vue précédemment est définie par VTmin = 26,1°C et VTmax = 38,4°C. Selon une étude très préliminaire, la température d’activité minimale (Tb MinActivity sensu Sinervo) pourrait s’abaisser à 14.0°C (Sinervo unpub.) mais cette hypothèse n’est considérée pour le moment.

Les modèles ont été déployés du 21/07/2020 à 18h00 au 15/08/2020 à 06h00 sur un pierrier au nord du col d’Aubisque (Latitude : 42,985186° ; Longitude : -0,335362° ; Altitude : 1695 m) à la surface et à 3 profondeurs dans le pierrier, -10 cm, -20 cm et -30 cm. Deux réplicats ont été réalisés, soit un total de 8 modèles déployés.

Connaissant la gamme de températures au sein de laquelle un animal est actif naturellement, gamme encadrée par les températures volontaires d’activité mini et maxi (VTmin = 20,8°C et VTmax = 35,2°C), il est possible d’estimer le temps d’activité journalier théorique. Par définition ce temps d’activité correspond à la période de temps, pendant une journée du cycle d’activité (nuit exclue donc, puisque l’espèce est diurne), durant laquelle la température opérante enregistrée est comprise dans la fenêtre [VTmin, VTmax]. Il peut être exprimé en nombre d’heures d’activité, notées Ha ou en pourcentage. Le temps de restriction est alors le temps restant, lorsque Te>VTmax ou Te<VTmin. Ces valeurs aux différentes profondeurs du pierrier (versus celle de surface) ont été comparées pour la journée la plus chaude de la période, la journée la plus froide et une journée “moyenne”.

Nous avons d’abord caractérisé les variations de températures opérantes (Te) en calculant leurs statistiques descriptives. Puis nous avons étudié les éventuelles journées remarquables, celles où les températures atteignaient VTmax et celles où les températures descendaient en dessous de VTmin. Enfin, nous avons calculé les pourcentages de temps d’activité et de restriction encadré par VTmin et VTmax comme expliqué plus haut.