Les Sentinelles du Climat

L’objectif des analyses statistiques préalables suivantes est de vérifier que les conditions d’observation de terrain ont été similaires entre les différents suivis afin d’étudier spécifiquement les effets du changement climatique sur les abondances des espèces.

L’analyse de Redondance (RDA) a été utilisée pour déterminer le pourcentage d’explication des variables de réponse des abondances par espèces de lépidoptères et d’ascalaphes par type de milieu par rapport aux variables explicatives suivantes : codeObs (nom de l’observateur), site_Ref_Courte (dénomination du site), dateDebut (date de début de l’observation), horaireDeb (horaire de début de l’observation), latitudeY1, longitudeX1 (localisation géographique du début des transects), temperatur (température de l’air en °C), humidite (humidité relative de l’air en %), ventBeaufo (échelle de Beaufort de 0 à 12), couvNuage (couverture nuageuse en %). Cette analyse a permis de déterminer les variables les plus explicatives par type de milieu, les variables climatiques et les espèces qui semblent être influencées par ces variables. L’association entre les variables climatiques et les abondances de ces espèces a été ensuite évaluée par le test de Spearman (p=0.05).

Pour chaque site, les dispersions des données climatiques via le coefficient de dispersion ont été comparées entre les différents types de milieux par mois d’observation de mai à août (Test de Kruskall et Wallis p=0.05) et entre les différents mois de suivis par année (Test de Friedman p=0.05). Les conditions climatiques sont présentées en moyenne plus ou moins l’écart type. Sur l’ensemble des sites, les conditions climatiques ont été comparées sur les années 2017-2018 pendant la période de mai à août pour chaque type de milieu (Test de Wilcoxon p=0.05).

47 sites de milieux humides (appelés « MH »), secs (appelés « MS ») et de montagne (appelés « MM ») ont été sélectionnés et validés en 2018. Les milieux humides correspondent aux landes humides, les milieux secs aux pelouses sèches et les milieux montagnards aux pelouses thermophiles montagnardes.

Quatre à huit passages pour les 47 sites sont effectués par an (de 2016 à 2019) et sont répartis entre mai et août (avec la possibilité d’avoir quelques passages fin avril à début septembre), soit 1 à 2 passages par mois. Les suivis sont espacés d’au moins 15 jours. Ils sont réalisés d’une année sur l’autre aux mêmes semaines . Les lépidoptères des landes humides dites « MH » (pour milieu humide) ont été suivis en 2019 sur 14 sites, contenant 36 transects parcourus de mai à août. Les transects sont parcourus lors de 73 journées dans l’année. Cela conduit à une longueur de 36 km. La largeur d’observation est de 5 m soit 18 ha de surface observée sur 4 mois. Les lépidoptères de pelouses sèches dites « MS » (pour milieu sec) ont été suivis en 2019 sur 23 sites contenant 59 transects parcourus de mai à fin août. Les transects sont parcourus lors de 114 dates dans l’année soit au total 27,1 km parcourus ou 13 ha de surface dénombrée dans les pelouses sèches. Les lépidoptères de pelouses de montagne dites « MM » (pour milieu de montagne) ont été suivis sur 10 sites contenant 12 transects de 150 m soit 6 km parcourus ou 3 ha de surface échantillonnée.

À partir de ces données, le nombre d’individus observés N_i est réparti dans des espèces à l’échelle du site, du milieu, de la région. Ces valeurs sont résumables par des proportions relatives p_i = N_i/N, N nombre total de maille à l’échelle de l’analyse. À partir de cette base de données, l’utilisation des « nombres de Hill » permet de généraliser les indices de diversité classiquement utilisés en écologie (diversité de Shannon-Wiener, indice de Simpson, Berger-Parker, etc.) et de les inclure. Le tracé du profil de diversité (proportions relatives des espèces rares, communes, abondantes) est une conséquence de ce choix d’analyse. L’utilisation des nombres de Hill s’impose peu à peu en méthode dans la littérature depuis quelques années . Le profil de diversité permet d’évaluer l’effort d’échantillonnage par le taux de couverture aux différentes échelles. En 2017 – 2019, les taux de couverture pour caractériser les milieux à partir des diversités des sites secs, humides et de montagne de la région sont de l’ordre de 99 % donc le nombre de sites est suffisant pour obtenir des résultats significatifs par type de milieu. De plus, les conditions d’observation sont réalisées dans des conditions météorologiques similaires entre les suivis à environ 90% (Test de permutation : Pseudo F=0.89, 0.84, 1.1, p<0.002) .

Le profil de diversité est décrit par une courbe décroissante contenant un paramètre « q » variant de 0 à l’infini, dit « d’indifférence aux espèces rares » . Parmi les espèces, celles qui sont rares sont potentiellement les plus fragiles dans le milieu. Ces profils sont des représentations graphiques exprimés en nombre d’espèces S(q) partant du nombre total d’espèces observés S(0) avec q=0. Puis plus la valeur de « q » est grande et moins les espèces rares y sont prises en compte. Le paramètre q supérieur à 2 (par exemple q = 4) est proche de l’équivalent traduit en nombre d’espèces de la part en individus de l’espèce la plus abondante S(q=infini). Avec q = 2 (équivalent à l’indice de Simpson), il s’agit en nombre équivalent de la part des espèces communes et abondantes S(2). Pour q = 1 (soit l’indice de Shannon), le nombre d’espèces renferme les espèces rares, communes et abondantes S(1). Enfin, pour q = 0, cela revient au nombre total d’espèces observées très rares, rares, abondantes .

Entre la richesse intrinsèque en nombre d’espèces observées S(0) = D₀ et la part en nombre d’espèces équivalentes induite par l’espèce la plus abondante mesurée par D(infini)= D_inf (indice de Berger-Parker) la forme du profil allant des espèces très « rares » à l’espèce la plus nombreuse est assimilable à un indice d’état de la fragilité du milieu. La forme générale retenue est donnée par une fonction F générale qui dépend de q :

3 variables principales du profil sont définies : nombre d’espèces D₀ ; part de l’espèce la plus abondante D_inf ;b répartition des abondances relatives des espèces. L’intérêt du lissage de cette courbe est de la résumer alors dans l’intervalle principal 0<q<4 par 3 paramètres de courbure au lieu de 3 indices ponctuels (q=1 Shannon, q=2 Simpson, q=4). Les profils évoluant vers des formes très concaves près de q=0, donc avec beaucoup d’espèces très rares, ont une fragilité intrinsèque : une espèce peu nombreuse a plus de chance de disparaitre statistiquement simplement par son nombre d’individus dans les températures extrêmes . Plus le paramètre b est important et plus le profil est concave et plus les espèces rares y sont relativement peu nombreuses donc vulnérables par exemple pour les températures >35°C ou inférieures à -4°C. Le milieu sec se situe à b = 0.76, le milieu de montage b = 0.52 et le milieu humide b = 0.11. Ensuite, les paramètres m₀ et m_inf corrigent et caractérisent encore plus finement la courbure pour détecter où se trouve le point spécifique de fragilité à surveiller. Plus m₀ est grand et plus les espèces très rares parmi les espèces « rares » sont concernées. Ces valeurs m₀ et m_inf sont proches pour les 3 milieux. Les milieux sont ainsi caractérisés par des paramètres principaux et d’autres secondaires qui s’emboitent.

Les tracés de profils de diversité et les figures d’interpolation et extrapolation avec une incertitude à 95 % de confiance estimée à partir de 50 bootstrappings sont obtenus dans ce texte avec le logiciel iNEXT Online : « software for interpolation and extrapolation of species diversity » .

Sur chaque site, des stations météorologiques ont été installées. Les analyses des données sont illustrées dans ce texte sur les années 2017 à 2019 entre le 1^er janvier jusqu’au 31 décembre pour les sites de milieux humides, de milieux secs et de montagne.

L’année civile correspond à un cycle thermique de source d’énergie de croissance des lépidoptères. Le cumul de températures dans le temps permet de renouveler les individus et de passer du stade latence de début d’année aux différents stades de développement de leur cycle biologique (œuf, chenille, chrysalide, imago) avec reproduction des individus au cours de l’année jusqu’à retrouver une nouvelle plage de latence en décembre. La sortie de léthargie est calculée lors du premier trimestre et l’entrée en léthargie correspond à la fin décembre correspondant à 99% des degrés-jours de l’année. La relation entre humidité relative mesurée et le calcul du point de rosée est utilisée.

Les données horaires d’une sonde de terrain représentent 24 données de température, d’humidité relative de l’air et donc du point de rosée sur 365 jours soit 24x3x365=26280 données annuelles.

La démarche de l’analyse de ces données consiste à les traiter statistiquement. Afin de comparer les signatures climatiques entre les différents types de milieux naturels, les données manquantes sont complétées par la méthode d’imputation multiple de Markov Chain Monte Carlo (MCMC) aléatoirement par année et par site : pour l’année 2017 milieu MS 3 MCMC, MM 4 MCMC ; pour l’année 2018 MH 9 MCMC, MM 13 MCMC ; et l’année 2019 MH 9 MCMC, MM 13 MCMC. Le test de normalité Test de Shapiro-Wilk montre que les données météorologiques températures, humidité relative et point de rosée ne suivent pas une loi normale (p<0.0001). Les tests statistiques appliqués par la suite sont des tests non paramétriques. L’ensemble des tests statistiques sont réalisés sur Excel, avec le logiciel R (version 3.6.0) et XLSTAT (version 19.02).

Nous nous concentrons ensuite sur un plan biologique sur une analyse à partir des degrés-jours de croissance. C’est une mesure couramment utilisée de l’accumulation thermique, c’est le lien thermo-physico-chimique entre le changement climatique et la phénologie des espèces liée aux réactions enzymatiques .

Les degrés-jours caractérisent l’énergie nécessaire pour l’évolution des espèces par exemple pour un papillon passant d’un stade à un autre (œuf, chenille, chrysalide, imago). D’après la littérature scientifique, ce développement se met en place principalement dans la plage de 10 à 35°C . Les différentes classes de température agissant sur la faune et la flore se résument alors à :

• En dessous des températures moyennes des jours inférieures à 5°, le froid met les réactions biologiques en latence,
• 5 à 10°C est la plage d’initiation des réactions biologiques,
• 10°C à 30°C correspond au développement de la faune et la flore avec une vitesse exponentielle peu prononcée assimilable à une droite justifiant l’approximation par la cumul de degrés-jours,
• 30°C à 35°C entraine un déclin du développement et une probabilité de létalité,
• >35°C correspond à la plage létale. La létalité d’un organisme est amorcée vers 30°C et est quasi-complète à 40°C.

En résumé, le taux de développement (nombre d’individus généré) par unité de temps des insectes soumis à la température obéit à des règles générales quasiment indépendantes des types d’espèces. L’ensemble d’un taxon réagit quasiment comme un système. Les vitesses de développement s’activent à partir de 5°C, croissent en fonction de la température entre 10°C et 30°C et décroissent ensuite . Les températures en dessous de 10°C ont en marge un rôle particulier nécessaire de « repos ».

Les températures minimales et maximales des journées sont lissées sur l’année. Les cumuls des degrés-jours de l’année sont calculés à partir des températures horaires. Elles sont aussi déductibles des lissages de ces températures minimales, maximales des journées sur l’année auxquelles il faut ajouter un profil moyen des températures de ces journées. Le lissage non linéaire de la courbe de cumul des degrés jours sur l’année donne un coefficient de corrélation R² de 0.9. La forme utilisée pour le cumul des degrés-jours est celle de Weibull :
J jour de l’année, t₀ jour de début d’activité, t_m jour d’apport maximal d’énergie.

La courbe du cumul des degrés-jours est ainsi résumée par 4 paramètres principaux (t₀, m, t_m, DJC_max):

• t₀ est le déclenchement du développement des individus; en pratique t₀=0,
• la composante m est liée à la vitesse de montée en température de la période chaude (forme) ;
• t_m le jour de la vitesse de montée maximale (échelle) ;
• DJC_max est la valeur asymptotique de fin d’année.

Les données météorologiques des sites d’étude sont celles de la maille de 1km/1km (maillage utilisé dans la modélisation précédente) : entre 1950 et 2020 ce sont les données dites Météo France, puis en prolongement dans le temps celles des scénarios climatiques RCP26, RCP45 et RCP85 entre 2021 jusqu’à 2100. Les variables quotidiennes d’analyses retenues sont notées : température minimale Tmin (Tn), température maximale Tmax(Tx), température du point de rosée (Td), température moyenne de journée Tmoy(Tm), et Degrés jours DJ, dont « DJC » les degrés-jours de croissance entre 10 et 35°C, « DJ30 » les degrés-jours au-delà de 30°C létalité statistique des individus et « DJ0 » les degrés-jours en deçà de 0°C avec une léthalité principale en dessous de -4°C .

L’analyse « mécanistique » consiste ainsi à agréger les données météorologiques en suivant une démarche sous-jacente de compréhension des phénomènes. Il s’agit de se rapprocher de variables intensives, en lien premier avec la croissance potentielle des individus (DJC) donc à leur abondance observée sur le terrain dans le temps et l’espace, puis avec le nombre d’espèces identifiées, leurs proportions relatives en individus, soit une variable de la richesse du milieu et leur fragilité dans le milieu. Ces proportions mesurent aussi le niveau de l’effort de terrain. Les degrés-jours >30°C et < à 0° interviennent en complément comme une pression de mortalité stochastique modifiant le résultat de la croissance. La démarche d’analyse propose systématiquement le développement de lissage des données pour alimenter l’analyse mécanistique. Les variables des températures, (ou degrés-jours) sont ainsi d’abord lissées année par année. Puis les 4(ou 2) paramètres de la fonction de lissage choisie, résumant l’année, sont à leur tour lissés dans le temps entre 1950 et 2100 par un polynôme de degré 1 ou 2 identifiant d’autres paramètres d’évolution moyenne entre 2021 et 2100. L’année thermique de développement des espèces est donc obtenue à partir des données horaires des températures. L’évolution selon les scénarios climatiques utilisent les données journalières principales. Pour chaque ensemble de points des températures maximales et des températures minimales de ces journées, une fonction sinusoïdale modifiée (à la puissance x) est retenue . Les degrés-jours utilisent une forme de cumul de Weibull à 2 ou 3 paramètres.

Les 8 paramètres climatiques sur l’année sont mis en corrélation (corrélation de Spearman p=0.05) avec les 6 paramètres qui résument le profil de diversité sur l’année pour chaque site ou chaque milieu. Ces derniers représentent sa richesse avec le nombre d’espèces, sa fragilité par le profil d’abondance (nombre d’individus par espèce). Les différences entre les valeurs des paramètres sont également comparées par année par le test de Wilcoxon sur données appariées (p=0.05). La mise en place du réseau de stations météorologiques à partir de l’année 2016 permet en exercice pour cette année 2019 d’analyser 7 sites de suivis en landes humides et en pelouses calcicoles sur les années 2017 et 2018. A la fin du programme en 2021, 10 sites pourront être analysés sur 4 années, 26 sites sur 3 années et 8 sites sur 2 années.

L’ensemble des tests statistiques de ce chapitre sont réalisés avec le logiciel R (version 3.6.0) et XLSTAT (version 19.02).

Les données de présence des espèces en Nouvelle-Aquitaine proviennent de l’Observatoire Aquitain de la Faune Sauvage (OAFS) . Elles comportent un total de 1413 espèces de lépidoptères, sous-espèces ou genres répertoriés.

Une large partie cependant ne comporte pas assez d’observations (nombre de mailles de présence < 20) pour permettre de réaliser un modèle corrélatif avec n=1188 espèces. Un premier tri de la base de données a été effectué afin d’éliminer les observations dupliquées ainsi que les observations imprécises, qui s’arrêtent au genre par exemple. Les observations de l’Apollon (Parnassius apollo) ont été fusionnées avec les observations de la sous-espèce Parnassius apollo pyrenaicus car nos données sont localisées dans les Pyrénées. Parmi les espèces restantes, 114 ont été sélectionnées car elles font partie des cortèges d’espèces de terrain observées dans le cadre du programme sentinelles du climat. Afin de limiter les effets du biais d’échantillonnage et de l’autocorrélation spatiale, les données de présence ont été agrégées par maille de 1km², les informations sont donc des présences par mailles.

Dans cette étude, ces espèces sont associées aux traits biologiques liés à l’habitat, la taille, la couleur, la valence, la capacité de vol, l’alimentation des chenilles, les périodes de vol, le nombre de générations, le risque climatique et le classement de la Liste rouge Aquitain.

La sélection des variables environnementales est une étape cruciale, car la pertinence des modèles en dépend. De nombreuses méthodes existent pour extraire celles qui participeront au modèle. Un premier choix de variables pouvant influencer les répartitions des espèces a été effectué en combinant les données climatiques, les données sur l’occupation du sol (végétation) et les données topographiques. Les variables les plus corrélées entre elles (corrélation de Spearman, p=0.05) ont été supprimées au profit de celles qui avaient le plus de sens écologique du point de vue de nos connaissances actuelles sur ces espèces et qui pouvaient être pertinentes pour l’ensemble des espèces étudiées.

Les variables climatiques proviennent des simulations Aladin 52 CNRM 2014 téléchargées via le site de Météo France DRIAS à partir desquelles des indices ont été calculés. Les données sont à une résolution de 8 km² redécoupée en mailles de 1km² et les valeurs sont des moyennes mensuelles. La période du présent est définie par les années allant de 1991 à 2020 et les horizons futurs sont définis comme suit : Horizon 1 (H1) = 2021-2050, Horizon 2 (H2) : 2051-2070 et Horizon 3 (H3) : 2071-2100. Ces horizons ont été déterminés d’après les recommandations de Météo-France qui préconise une durée de 30 ans afin de lisser les « bruits » inclus dans les valeurs des simulations climatiques . Les simulations dans la période dite du présent contiennent pour les années 2006-2020 des données du scénario RCP 8.5 et les données précédentes de 1991 à 2005 sont des données historiques. Pour chaque « Horizon », un scénario climatique « RCP » (Representative Concentration Pathway) : Profils représentatifs d’évolution possible de forçage radiatif lié entre autres aux concentrations en équivalent CO₂ est associé. Pour les simulations Aladin, ils sont au nombre de trois RCP2.6, RCP4.5, RCP8.5 et représentent chacun une trajectoire différente d’émissions et de concentrations de gaz à effet de serre, d’ozone et d’aérosols. Ces projections prennent aussi en compte des évolutions socio-économiques (adaptation, réduction, stabilisation ou dépassement).

Une première sélection de variables est basée sur les connaissances de biologie et d’écologie des espèces, puis afin d’éviter d’utiliser des variables corrélées entres elles un test de corrélation de Spearman (non-paramétrique : les variables climatiques ne suivent pas la loi normale) est effectué. Ainsi cinq variables sont conservées : Moyenne de l’Humidité Relative en été (MoyHR_ete), le nombre de jour de neige dans l’année (NjN_annee), le nombre de jour de pluie au printemps et en été (NjP_print_ete), la somme des degrés-jours entre 10°C et 35°C du 1er janvier à la fin du printemps (TotDJ1035_print) et la somme des degrés jours supérieurs à 30°C de l’année (TotDJ30_annee). Les degrés-jours caractérisent l’énergie nécessaire pour la croissance des individus, par exemple pour un papillon le passage d’un stade à un autre (œuf, chenille, chrysalide, imago). Ces degrés-jours sont calculés à partir des formules mathématiques décrites par F. Mallard approximant linéairement l’évolution en température du taux de croissance.

En plus des variables climatiques, celles décrivant l’occupation du sol ont été utilisées (exprimées en pourcentage de couverture sur un pixel de 1km²) : le pourcentage de forêt, le pourcentage de culture, le pourcentage de prairies, le pourcentage de landes ligneuse et de pelouses, le pourcentage d’eau. Les variables d’occupation du sol sont calculées à partir du « CES Occupation des sols (OSO) » et produites à partir d’images satellites (2019) d’une résolution de 25m. Les 23 nomenclatures sont fusionnées en 11 classes. La variable topographique sélectionnée est l’altitude (données IGN).

Toutes les analyses sont réalisées à partir du logiciel R (version 3.6.3, R Core Team, 2020). Les modèles ont été effectués avec le package BIOMOD2 (version 3.4.6)  qui permet de créer des modèles d’ensemble à partir de différents algorithmes et simulations. Pour cette étude, les algorithmes GLM, RandomForest et Maxent  ont été sélectionnés et pour chaque algorithme trois répétitions ont été effectuées afin d’évaluer la performance (80 % des données pour l’entraînement du modèle et 20% pour l’évaluation). Les points de fond ont été générés aléatoirement sur la zone d’étude Nouvelle-Aquitaine (n=10 000). Le seuil pour la création des modèles d’ensemble a été choisi comme le quantile 0,5 des valeurs TSS des 24 modèles individuels créés. Si tous les modèles ont une valeur TSS inférieure à 0,4 alors aucun modèle d’ensemble n’est créé car les modèles individuels ne sont pas assez performants. Les modèles d’ensemble représentent la moyenne des prédictions des modèles individuels pondérée par la performance de chacun d’eux. Les prédictions pour la région Nouvelle-Aquitaine ont été réalisées pour le présent et pour chaque scénario et « Horizon ». Elles sont exprimées en probabilité relative comprise entre 0 (plus mauvais habitat) et 1 (meilleur habitat). Des cartes de présence/absence sont aussi créées avec le seuil TSS (fonction rangesize du package BIOMOD2). Ainsi, le nombre de mailles occupées par les espèces est calculé pour le présent et pour chaque scénario- « Horizon ». Ces cartes sont utilisées pour connaître le nombre de mailles passant de « présent » à « absent » (perte), d’« absent » à « présent » (gain).

Le taux d’espèces en diminution et stabilisation/progression modélisé dans le programme sur la région Nouvelle-Aquitaine (SDCR) est comparé aux résultats de l’atlas du risque climatique européen des lépidoptères sur la zone européenne (EU), sur la région Nouvelle-Aquitaine (EUR) . Cette comparaison est réalisée pour les scénarios les plus pessimistes : pour le SDCR période 2071-2100 pour le RCP 8.5 (réchauffement compris entre 2.6°C–4.8°C) ; et pour le EU et EUR période 2051-2080 pour GRAS (GRowth Applied Strategy) se rapprochant du scénario A1FI (augmentation moyenne attendue de la température de 4,1 °C) . Ces résultats sont également comparés aux résultats à dire d’expert (PYG) réalisé à partir de l’estimation de la répartition des espèces dans la région pour le scénario le plus pessimiste en 2100 d’après les connaissances de terrain de l’expert P. Y. Gourvil Coordinateur du Plan Régional d’Actions en faveur des lépidoptères patrimoniaux en Nouvelle-Aquitaine.