Navigation – Plan du site
2015
722

Identification des variables expliquant la distribution spatiale d’oiseaux de la forêt boréale et modélisation de tendances futures : une approche multivariée

Jonathan Gaudreau, Liliana Perez et Pierre Legendre

Résumés

Les changements climatiques prennent une importance grandissante dans l’étude de la dynamique des populations animales. Plusieurs experts s’entendent pour affirmer que les changements climatiques seront un des principaux moteurs de changement écologique dans les prochaines décennies. L’objectif de cette recherche est l’identification des principaux facteurs responsables de la distribution spatiale d’oiseaux de la forêt boréale au Québec, afin de proposer des modèles de distribution d’espèces suivant les changements climatiques prévus. Deux approches multivariées sont employées : l’analyse de redondance canonique (RDA) et le partitionnement de variation. Au total, 39 espèces d’oiseaux sont sélectionnées, en plus de variables bioclimatiques, anthropiques, écoforestières et d’élévation. Les variables bioclimatiques sont responsables de 53% de la variation dans la répartition spatiale des oiseaux étudiés, les variables ressources, c’est-à-dire l’altitude et le pourcentage de milieu humide, de 5%, et les variables anthropiques de moins de 1%. Les résultats des modèles démontrent que les deux espèces modélisées verront leur répartition spatiale fortement modifiée par les changements climatiques et se déplaceront vers les latitudes plus septentrionales ou en altitude, suivant l’intensité du réchauffement.

Haut de page

Texte intégral

Les auteurs voudraient remercier le Centre d’Étude de la Forêt (CEF) et l’Université de Montréal, qui ont offert leur support financier, les ornithologues amateurs et professionnels ayant mis au point la très complète base de données eBird, tous les organismes ayant fourni des données, et ce, complètement gratuitement. Les auteurs voudraient aussi remercier particulièrement tous les chercheurs contribuant à la démocratisation des outils scientifiques, libres et participatifs. Les auteurs sont aussi très reconnaissants pour les commentaires fournis par deux réviseurs anonymes dans une version précédente de ce manuscrit.

Introduction

1Les changements climatiques prennent une importance grandissante dans la recherche en environnement et pourraient être un des moteurs principaux de changement écologique dans les prochaines décennies (Parmesan 2006; Lawler et al. 2009; Virkkala et al. 2008). Dans son rapport de 2007, le Groupe d'experts intergouvernemental sur l'évolution du climat (GIEC) conclut que la plupart des systèmes naturels seront touchés à différentes intensités (IPCC, 2007). Les écosystèmes terrestres doivent s'adapter aux printemps hâtifs et la faune et la flore migrent, dans la plupart des cas, vers de plus hautes latitudes ou altitudes (Forister et al., 2010 ; Hurlbert & Liang, 2012 ; IPCC, 2007 ; Walther et al., 2002). Les conclusions du GIEC, basées sur plus de 29 000 tableaux de données provenant de 75 études, affirment qu'un nombre important d'écosystèmes seraient perturbés (IPCC, 2007). Les impacts des changements climatiques sur le milieu physique comprennent, entres autres, la fonte des calottes glaciaires, le dégel du pergélisol et la vulnérabilité et l'instabilité des versants montagneux en zones de pergélisol. Parmesan (2006) fait état des multiples conséquences des changements climatiques anthropiques sur la faune et sur la flore. Parmi celles-ci se retrouvent les modifications dans la répartition spatiale des espèces.

2Des études similaires, réalisées en Amérique du Nord, convergent et concluent que des espèces migrent vers de plus hautes latitudes ou en altitude, lorsque la topographie le permet (Forister et al. 2010 ; Hurlbert and Liang 2012 ; Walther et al. 2002). Cette hypothèse est confirmée par des recherches portant sur l'Antarctique, l'Arctique, ainsi que sur les milieux tempérés de l'hémisphère nord (Parmesan 2006 ; Virkkala et al. 2010). Ainsi, l'augmentation des températures, en particulier dans les régions boréales et polaires, induit-elle une pression considérable sur les écosystèmes forestiers?

3La province de Québec, située dans le Centre-Est du Canada, montre un gradient latitudinal important et comprend de nombreux domaines bioclimatiques, allant des forêts tempérées à la toundra arctique. Environ 50% de son territoire est couvert par la forêt boréale, qui comprend quatre domaines bioclimatiques distincts (MERN, 2013). La forêt boréale du Québec, s'étendant sur plus d'un million de kilomètres carrés et abritant de 150 à 300 espèces d'oiseaux migrateurs, sera fortement affectée dans les prochaines années (MERN, 2013).

4Ainsi, afin de prévoir les répercussions de changements climatiques sur la répartition spatiale de la faune aviaire boréale, il importe d’identifier les principaux facteurs de cette répartition. Cette recherche vise donc, dans un premier temps, l’identification des variables explicatives influençant la répartition spatiale de 39 espèces d’oiseaux de la forêt boréale via l’utilisation de méthodes de modélisation statistique appliquées aux données multivariées, c’est-à-dire l’analyse de redondance canonique (RDA) et le partitionnement de la variation. Suite à ces analyses, deux espèces voient leur répartition spatiale modélisée via des approches de modélisation bioclimatique. De plus, la présente recherche vise à approfondir les connaissances sur les réponses fauniques aux changements climatiques, dans le but de mieux adapter les politiques de conservation. L’hypothèse émise est que la répartition spatiale des espèces est fortement corrélée au climat, plus particulièrement aux précipitations et à la température moyenne annuelle. Ainsi, les prédictions climatiques devraient favoriser une migration des espèces boréales vers de plus hautes latitudes ou en altitude en présence de topographie variable.

Matériel et préparation des données

5Les données nécessaires aux analyses proviennent toutes de bases de données ouvertes d’organismes (Global Forest Watch, 2009 ; Gouvernement du Canada, 2000 ; Hijmans et al., 2005 ; CanFI, 2001 ; eBird, 2013). Les outils nécessaires à l’analyse sont, pour la plupart, le langage de programmation et de statistiques R 3.1.0 (R Core Team, 2014) et la suite ESRI ArcInfo 10.1 (ESRI, 2012) pour certaines transformations cartographiques.

Aire d’étude et sites échantillonnés

6La recherche porte sur la forêt boréale du Québec, Canada (figure 1). Les deux types de modèles nécessitent deux types d’intrants différents. La modélisation par analyse canonique de redondance (RDA) nécessite un tableau de données dont les observations sont répertoriées par site. Ainsi, puisque les données explicatives sont disponibles pour l’ensemble du Québec en format raster, un échantillonnage aléatoire est effectué, afin de réduire le nombre de points d’observation. Au total, cinquante sites sont aléatoirement sélectionnés, dans l’ensemble des domaines bioclimatiques du Québec (figure 1), en utilisant l’outil de création de points aléatoires distribués de façon uniforme d’ArcInfo (ESRI, 2012). La table d’attribut, à laquelle sont liées les valeurs des abondances d’oiseaux et des variables explicatives, est ensuite exportée et utilisée dans le logiciel statistique R (R Core Team, 2014) pour les analyses statistiques multivariées. Pour la seconde partie, c’est-à-dire la prédiction via des approches de modélisation bioclimatique, les données pour l’ensemble du territoire étudié seront directement utilisées comme intrants. En effet, la librairie biomod2 (Thuiller et al., 2013) exige, en intrant, des couches raster.

Figure 1 : Sites échantillonnés

Figure 1 : Sites échantillonnés

Variables réponses : abondances d'oiseaux

7Les données d’abondance d’oiseaux proviennent de la base de données ouverte eBird (eBird, 2013), qui contient, pour le Québec uniquement, plus de 1.3 millions d'entrées. Pour chaque observation, plusieurs informations sont disponibles, notamment l'heure de l'observation, l'espèce, le sexe, la localisation, entre autres. Seuls les champs concernant les coordonnées spatiales, soit la latitude et longitude, et le nombre d’observations, sont conservés pour l’analyse. Les analyses statistiques multivariées (RDA et partitionnement de variation) nécessitent une étape supplémentaire. En effet, afin de permettre un échantillonnage aléatoire pour l'ensemble des espèces en utilisant les mêmes sites, l’aire maximale de répartition est d’abord interpolée pour chaque espèce. Cette interpolation produit une carte raster d’abondance d’espèces basée sur des observations réelles, simulant un polygone de répartition convexe, mais fonctionnant dans un espace matriciel, c’est-à-dire composé de cellules. La méthode utilisée est l’interpolation des points d’observations, en abondances, via la méthode nearest neighbors d’ArcInfo (ESRI, 2014). Le résultat donne une couche raster représentant l’aire maximale de distribution de l’espèce, fait en pixélisant l’abondance, plutôt que la présence/absence. L’exercice est répété pour chacune des espèces choisies, totalisant 39 cartes matricielles. Cette interpolation des valeurs est requise pour l’échantillonnage aléatoire des sites nécessaire aux analyses statistiques multivariées, qui demandent des sites communs pour les 39 espèces. Évidemment, lors de la seconde étape, c’est-à-dire la modélisation bioclimatique, les observations réelles de la base de données sont entrées directement dans le modèle sous forme de coordonnées latitude/longitude, puisqu’une seule espèce à la fois est utilisée comme variable réponse.

Variables explicatives

8Le choix des variables explicatives, ou variables indépendantes, ou variables prédictives, est primordial à ces analyses. En effet, le but des analyses multivariées est l’identification de variables spatialisées influençant la répartition spatiale des oiseaux de la forêt boréale. De ces variables, seules les variables significatives sont utilisées dans les prévisions. Les variables explicatives sont projetées selon le système géodésique mondial de 1984, ou WGS84 et consistent en cartes matricielles, ou raster, de 2112 lignes par 2721 colonnes. Le cadastre s’étend de 51.1 degrés ouest à 79.775 degrés ouest en longitude et de 44.98 degrés nord à 62.58 degrés nord en latitude. La résolution est de 0.0083 degré décimal par cellule, ou 787m par 787m, soit de moins d’un kilomètre carré. Les données bioclimatiques proviennent de la base de données WorldClim et s’étendent de BIO1 à BIO19 (voir l'annexe 1 pour la signification des sigles et des noms de variables) (Hijmans et al, 2005). Ces données, de nature quantitative discrète et à faible résolution spatiale (<1km), furent dérivées de températures et précipitations mensuelles pour engendrer des variables plus pertinentes en écologie, notamment la température moyenne estivale, l’étendue annuelle des précipitations, la saisonnalité des précipitations, etc. Les données représentant les conditions climatiques actuelles consistent en la moyenne de la période de référence de 1950 à 2000. Les données futures représentent l'an 2070 (moyenne de la période 2061-2080) pour les quatre profils d’évolution de la concentration, ou scénarios RCP (boîte 1), établis par le Groupe d'experts intergouvernemental sur l'évolution du climat (GIEC), dans son cinquième rapport (IPCC, 2013).

Boîte 1 : Scénarios RCP

À l'occasion de son 5e rapport, le GIEC définit un nouvel ensemble de scénarios climatiques, portant le nom de Representative Concentration Pathways, ou RCPs. Ces scénarios sont nommés par référence à leur forçage radiatif total pour l'année 2100 par rapport à 1750, année de référence préindustrielle. Les scénarios, nommés RCP2.6, RCP4.5, RCP6.0 et RCP8.5 représentent respectivement des forçages de 2.6Wm-2, 4.5Wm-2, 6.0Wm-2 et de 8.5Wm-2. Les simulations (figure 2) furent effectuées avec des concentrations de CO2 atmosphérique, pour l'année 2100, de 421 ppm (RCP2.6), 538 ppm (RCP4.5), 670 ppm (RCP6.0) et 936 ppm (RCP8.5) (IPCC, 2013).

Figure 2 : Tendances dans l’évolution des concentrations de gaz à effet de serre selon les 4 profils d’évolution de la concentration. De gauche à droite, on retrouve : le dioxyde de carbone (CO2) le méthane (CH4), et les oxydes nitreux (N2O) (Vuuren et al. 2011)

Figure 2 : Tendances dans l’évolution des concentrations de gaz à effet de serre selon les 4 profils d’évolution de la concentration. De gauche à droite, on retrouve : le dioxyde de carbone (CO2) le méthane (CH4), et les oxydes nitreux (N2O) (Vuuren et al. 2011)

9Les données d’élévation proviennent de la Base Nationale de Données Topographiques du Canada (BNDT) (Gouvernement du Canada, 2000). Cette carte matricielle affiche l’élévation, en mètres, pour l’ensemble du Québec. Aux berges, rives, océans et estuaires est attribuée une valeur de zéro. Le niveau marin moyen est assigné selon le système de référence altimétrique géodésique du Canada de 1928 et la résolution verticale est d’un mètre (Gouvernement du Canada, 2000). Les données de perturbations anthropiques du Québec proviennent de Global Forest Watch Canada et sont composées de quatre différents types de perturbation : 1) linéaires (routes, pipelines, lignes de haute tension), 2) de forme polygonale (forêt exploitée, puits de pétrole/gaz, terres agricoles), 3) réservoirs et 4) mines actives (Global Forest Watch, 2009). Aucun réservoir, aucune mine et aucune perturbation de forme polygonale ne se trouvent sur les sites aléatoirement sélectionnés, si bien que seules les perturbations linéaires sont conservées pour l’analyse. Finalement, les données de foresterie proviennent de l’inventaire forestier canadien de 2001 et indiquent le pourcentage de territoire boisé par hectare, en plus de contenir les informations relatives à la proportion des différentes essences d’arbres (CanFI, 2001). Les données sont de nature quantitative discrète (CanFI, 2001). De cette carte multicouche, deux rasters sont extraits, soient : 1) une carte montrant le pourcentage de forêt par cellule et 2) une carte représentant le pourcentage de milieux humides par cellule. Finalement, ces variables peuvent être regroupées en trois catégories : les variables ressources, comprenant les données forestières et d’élévation, les variables bioclimatiques, BIO1 à BIO19, et les variables anthropiques, comprenant les perturbations de type linéaire.

Méthodes

10Afin d’être en mesure de prédire la répartition spatiale future d’espèces d’oiseaux de la forêt boréale, une première étape consiste à identifier les variables responsables de cette répartition spatiale communes à la plupart des espèces. La première partie consiste en des analyses multivariées destinées à identifier les prédicteurs communs des 39 espèces d’oiseaux sélectionnées. Lorsque ces variables sont sélectionnées, une seconde partie consiste en la prévision de la répartition spatiale future pour deux espèces sélectionnées parmi les 39 d’origine.

Analyses multivariées pour identifier les variables responsables de la répartition spatiale des oiseaux de la forêt boréale du Québec

Transformation des données

11Les résultats de l’interpolation des points d’observation des oiseaux donnent des abondances allant de 0 (absence) à 720 (plus grande observation répertoriée, il s'agit de macreuses à tête blanche observées au sanctuaire d'oiseaux de l'île aux Basques). Un test de normalité de Shapiro-Wilks, dont l’hypothèse nulle consiste à affirmer que les données suivent une distribution normale et où l'hypothèse alternative implique une distribution asymétrique ou anormale, permet de conclure que les observations sont de distribution très asymétrique (Royston, 1982). Une transformation logarithmique des données réponses est donc effectuée afin de réduire l’asymétrie des données et, ainsi, d’améliorer la puissance des tests statistiques.

12L'analyse de redondance canonique est utilisée d'abord sur l’ensemble des tableaux de données (Legendre et Legendre, 2012). Par la suite, les variables colinéaires sont successivement enlevées du modèle, afin d'atteindre un facteur d'inflation de la variance (VIF) inférieur à 10 pour chacune des variables retenues. Afin de conserver les variables les plus représentatives, une sélection progressive bidirectionnelle est ensuite menée. Ce type d’approche effectue une sélection des variables pertinentes en fonction du critère d’information d’Akaike, ou AIC et de la représentativité (p-value) (Oksanen et al., 2013). Une sélection progressive par ajout, en fonction du R-carré ajusté, est aussi utilisée, afin de comparer le résultat de la sélection (Oksanen et al., 2013). Ces trois analyses sont répétées jusqu’à ce que seules les variables significatives soient intégrées au modèle.

1) Analyse de redondance canonique

13Une analyse de redondance canonique est un type d’analyse asymétrique, c’est-à-dire où on retrouve un tableau réponse Y et un tableau explicatif X. Les données subissent d’abord une régression linéaire multiple, afin de calculer la matrice des valeurs ajustées Ŷ, avec l'équation suivante (Legendre et Legendre, 2012. eq.11.11) :

Ŷ = X[XʹX]-1XʹY

où X représente le tableau explicatif et Y, le tableau réponse. Par la suite, les valeurs ajustées font l’objet d’une analyse en composantes principales (ACP) et sont projetées dans un plan cartésien (ter Braak, 1986). La RDA permet de projeter trois types d’objets, soient les sites, les variables réponses et les variables explicatives, dans un seul graphique d’ordination, rendant la lecture beaucoup plus simple et permettant une meilleure compréhension des relations entre Y et X (Legendre et Legendre, 2012). Plusieurs tests sont effectués sur les résultats de la RDA afin de vérifier la signification des résultats. Les axes canoniques et les prédicteurs peuvent être soumis à des tests par permutations pour en calculer la statistique F et la probabilité associée.

2) Réduction de la colinéarité par suppression de variables

14Afin de bien identifier les variables expliquant la répartition spatiale d’espèces d’oiseaux, il est nécessaire de bénéficier de variables explicatives indépendantes les unes des autres. Ainsi, il est crucial de réduire la colinéarité. On parle de colinéarité lorsque deux ou plusieurs variables explicatives sont corrélées, partageant une relation linéaire. Une des mesures utilisées pour identifier la colinéarité est le facteur d’inflation de la variance (ou VIF), qui se calcule en menant, successivement, des régressions linéaires de chacune des variables explicatives sur l’ensemble des variables explicatives. Les variables bioclimatiques offertes par WorldClim étant basées sur des données de température et de précipitation, elles sont donc toutes construites sur les mêmes données, et donc très corrélées. Le facteur d'inflation de la variance est obtenu avec l'équation suivante (Legendre et Legendre, 2012. Eq. 10.7) :

où (j) est la variable explicative dont on désire connaître le VIF. Le seuil maximal est fixé à 10, puisqu’un seuil de 5 ne permet pas, dans le cas présent, de conserver assez de variables environnementales. Le retrait des variables à haut facteur d'inflation de la variance permet de réduire la colinéarité en supprimant les variables colinéaires, permettant ainsi de bien identifier l'influence individuelle des variables explicatives. Les variables abandonnées sont celles ayant le plus haut VIF et étant le moins pertinentes du point de vue écologique.

3) Partitionnement de variation pour séparer les variables en groupes

15Afin d’identifier les facteurs influençant davantage la répartition spatiale des oiseaux de la forêt boréale, les variables sont regroupées en trois groupes, soient birds.envClim, comprenant les variables bioclimatiques, birds.envRes, comprenant les variables ressource, soient l’élévation, le pourcentage de forêt et le pourcentage de milieu humide, et birds.envAnth, comprenant les perturbations anthropiques de type linéaire (Oksanen et al., 2013).

16Le partitionnement de variation (Borcard et al., 1992) est représenté par un diagramme de Venn, tel que présenté à la figure 3. Ce diagramme permet d’identifier la contribution individuelle de chaque groupe de variables. La partie [a] représente la part de la variation expliquée exclusivement par le groupe de variables représenté dans le cercle de gauche, [c] représente la part de la variation expliquée uniquement par le groupe de variables représenté dans la sphère de droite. Finalement, la partie [b] représente l’intersection entre les deux groupes, soit la portion de la variation expliquée conjointement par les deux groupes de variables.

Figure 3 : Diagramme de Venn pour partitionnement de la variation de Y en fonction de deux groupes de variables explicatives

Figure 3 : Diagramme de Venn pour partitionnement de la variation de Y en fonction de deux groupes de variables explicatives

Modélisation bioclimatique pour modéliser les distributions d’espèces

17Des approches de modélisation bioclimatique sont utilisées pour produire des modèles de distribution d’espèces (de l’anglais Species Distribution Models, SDMs), en utilisant la librairie R biomod2, une plateforme de modélisation, de calibration et d’évaluation de modèles de distribution d’espèces, ou de niche (boîte 2) permettant d’avoir recours à plus de dix algorithmes de classification et de régression (Thuiller et al., 2013).

Boîte 2 : Postulat d’équilibre et concept de niche

Il est important de définir les deux concepts de ces modèles statiques : le postulat d’équilibre et le concept de niche. Le postulat d’équilibre est défini comme suit : puisque les données environnementales et d’abondance d’espèces sont normalement échantillonnées sur une période de temps limitée ou sur des sites d’étude limités, les modèles ne reflètent qu’un cliché temporel de la relation attendue entre les conditions environnementales et la répartition des espèces. Ainsi, pour plus de commodité, on assume que les espèces modélisées se trouvent en pseudo-équilibre avec leur environnement (Guisan et al., 2005, traduit de l’anglais) Le second concept à définir est celui de niche. Les modèles de distribution d’espèce, plus particulièrement les cartes résultantes, reposent sur une combinaison entre la niche écologique de Grinnell, où une espèce se trouve partout où les conditions environnementales sont favorables à sa reproduction et la niche réalisée de Hutchinson, où une espèce est exclue de sa niche fondamentale par prédation ou compétition (Guisan et al., 2005). Cette combinaison est due au fait que les modèles de distribution d’espèce, ou modèles de niche, utilisent des observations réelles pour identifier la niche écologique d’une espèce. Or, les observations réelles sont le fruit de la compétition et de la prédation des espèces; la niche réalisée sert donc indirectement à la conception des modèles.

18Pour la modélisation des distributions des deux espèces sélectionnées, cinq algorithmes furent utilisés, soit : RF, MARS, MAXENT, FDA et GLM (tableau 1).

Tableau 1 : Algorithmes de classification et de régression utilisés pour la génération des modèles de distribution d’espèces

Tableau 1 : Algorithmes de classification et de régression utilisés pour la génération des modèles de distribution d’espèces

19Tout d'abord, pour chaque espèce, les données de présence sont combinées à 10 000 pseudo-absences générées de façon aléatoire sur l'ensemble du territoire. Puis, ces présences/absences sont superposées aux variables environnementales et divisées en un ensemble contrôle (80%), utilisé pour générer un modèle de niche, et en un ensemble témoin (20%) servant à valider le modèle de niche par le biais de métriques de discrimination. Pour éviter un trop important taux d’erreur, le seuil minimal d’aire sous la courbe ROC (de l’anglais area under the receiver operating characteristic curve (AUC)) est définit à 0.7 (rappel : une classification aléatoire aura un score AUC de 0.5 et une classification parfaite aura un score AUC de 1). Puis, les futures conditions environnementales sont ajoutées au modèle pour remplacer les conditions actuelles et créer des cartes de distribution futures. Pour définir la présence ou l'absence d'espèce, une approche de modélisation par ensemble (de l’anglais Ensemble Modeling) est utilisée (Thuiller et al., 2013). Cette méthode permet de calculer la moyenne des prédictions des différents algorithmes statistiques utilisés. Par la suite, un algorithme réalise une décision par comité (de l’anglais committee averaging), où les probabilités des différents modèles subissent une transformation en données binaires selon le seuil d’AUC > 0.7. Ainsi, les cellules ayant une valeur située entre 0 et 0.7 se voient attribuer la valeur de zéro et les cellules ayant une valeur supérieure ou égale à 0.7 se voient attribuer la valeur de un. Puis, chaque méthode vote en faveur ou en défaveur de la présence de l’espèce. Finalement, ce score est ramené entre 0 (absence unanime) et 1000 (présence unanime) et est attribué comme valeur finale. Ainsi, cette méthode permet de concilier les décisions binaires de présence/absence et l'incertitude reliée aux prédictions dans une même carte, témoignant de l’incertitude inhérente aux différents algorithmes de modélisation bioclimatique (Pearson and Thuiller 2006). Une valeur de 500 signifiera que 50% des modèles prédisent une présence et 50%, une absence (Georges et Thuiller, 2014). Le schéma 1 présente une version simplifiée du processus de modélisation bioclimatique utilisé pour une espèce.

Schéma  : Adapté de Berteaux et al., 2010. Représentation simplifiée de la méthode de modélisation bioclimatique utilisée.

Schéma  : Adapté de Berteaux et al., 2010. Représentation simplifiée de la méthode de modélisation bioclimatique utilisée.

Les chiffres entre parenthèses représentent les étapes dans l’exécution du code. Les parallélogrammes représentent des entrées/sorties, les rectangles, des processus et le cylindre bleu, une base de données. Le rectangle rouge représente le cœur du modèle, soit le modèle de niche.

Résultats

Sélection de variables

20Une première RDA menée sur l’ensemble des variables explicatives permet d’identifier plusieurs choses. Tout d'abord, le calcul des facteurs d'inflation de la variance permet d'identifier les variables les plus colinéaires et de les retirer manuellement. En incluant toutes les variables explicatives, certaines variables montrent un VIF extrêmement élevé, proche de 2000 (bio1, bio11, bio6). Le R-carré est de 0.854, mais le R-carré ajusté est de 0.735. Le test par permutation effectué sur les variables indique que les variables pct_for, bio13, bio2, bio3 et bio5 sont significatives (p<0.05). Le test par permutation effectué sur les axes canoniques indique que seuls les axes 1, 7, 11 et 15 sont significatifs (p<0.01) après 299 permutations.

21Afin de réduire le nombre de variables et augmenter la puissance des analyses, les variables colinéaires sont supprimées une à une, et ce de façon manuelle. Les variables conservées sont l'altitude, les perturbations linéaires, le pourcentage de forêt, le pourcentage de milieu humide et les variables bioclimatiques bio7, bio11, bio15 et bio16, représentant respectivement l'étendue annuelle des températures, les températures moyennes du trimestre le plus frais, la saisonnalité des précipitations et les précipitations du trimestre le plus arrosé. Le VIF moyen des huit variables sélectionnées est de 5.24, avec 11.7 comme valeur maximale (bio11), ce que nous jugeons acceptable. La RDA est ensuite effectuée une seconde fois, cette fois sur les variables non-colinéaires (figure 4). Le R-carré ajusté du modèle est de 0.58.

22Par la suite, une sélection bidirectionnelle de variables mène à la suppression des variables linear et pct_for. L'analyse finale est effectuée sur six variables, soient les quatre variables bioclimatiques (bio7, bio11, bio15 et bio16), le pourcentage de milieu humide et l'élévation (figure 5). Sur ce graphique, il devient clair que la variable bio11, soit la température moyenne du trimestre le plus frais, est la variable qui influence la distribution spatiale de la majorité des espèces étudiées. La seconde variable en importance est la variable bio16, ou les précipitations du trimestre le plus arrosé. À l'opposé, deux variables influencent négativement la présence et l'abondance des espèces : bio15 et bio7, soient la variation saisonnière des précipitations et l'étendue annuelle des températures.

Figure 4 : Résultats de la RDA sur les variables non-colinéaires (p<0.005, R2adj= 0.58)

Figure 4 : Résultats de la RDA sur les variables non-colinéaires (p&lt;0.005, R2adj= 0.58)

Figure 5 : Résultats de la RDA sur les variables sélectionnées (p<0.005, R2adj= 0.59)

Figure 5 : Résultats de la RDA sur les variables sélectionnées (p&lt;0.005, R2adj= 0.59)

23Un test par permutations effectué sur les variables explicatives permet de tester la signification des six variables sélectionnées (Tableau 2). Toutes les variables intégrées au modèle final sont hautement significatives (p=0.01), à l'exception de pct_wet, qui est faiblement significative (p=0.07).

Tableau 2 : Résultat du test de permutation sur les prédicteurs

Tableau 2 : Résultat du test de permutation sur les prédicteurs

24Finalement, le même test, mais réalisé sur les axes canoniques, indique que les trois premiers axes sont hautement significatifs (p<0.005) (Tableau 3) et que le quatrième axe est significatif (p=0.071).

Tableau 3 : Résultats du test de permutation sur les axes canoniques

Tableau 3 : Résultats du test de permutation sur les axes canoniques

Partitionnement de variation

25Pour le partitionnement de variation, les variables furent séparées en groupes. Les quatre variables bioclimatiques forment le groupe de variables bioclimatiques et les variables pct_wet et elevation formèrent un groupe de variables ressources. Puisque les perturbations linéaires sont rejetées du modèle, le groupe de variables anthropiques est abandonné. Les variables bioclimatiques sont responsables de 53% de la variation totale, soit directement (40%), soit dans l’information expliquée conjointement avec les variables ressources (13%) (Figure 6). Les variables ressources, quant à elles, ne sont responsables que de 18% de la variation totale, soit directement (5%), soit via l’information expliquée conjointement avec les variables bioclimatiques (13%). Finalement, les résidus représentent 41% de la variation.

Figure 6 : Partitionnement de la variation des données des 39 espèces d’oiseaux (log) par deux groupes de variables explicatives : les variables bioclimatiques et les variables ressource

Figure 6 : Partitionnement de la variation des données des 39 espèces d’oiseaux (log) par deux groupes de variables explicatives : les variables bioclimatiques et les variables ressource

Modélisation bioclimatique

26Les deux espèces sélectionnées pour les prédictions sont la mésange à tête noire (BCCH) et le gros-bec errant (EVGR), de par leur très haut R-carré ajusté (R2adj = 0.7943 et R2adj = 0.79 respectivement) obtenu lors de régressions linéaires multiples sur les variables explicatives sélectionnées lors de la RDA. Dans les deux cas, les espèces verront leur limite nordique s'étendre vers de plus hautes latitudes à une intensité différente selon le scénario climatique utilisé (figures 7 et 8).

27Les cartes sont obtenues par consensus entre les différentes techniques. Ainsi, une valeur de 1000 représente l'accord total sur la présence de l'espèce, une valeur de zéro le désaccord total et toutes les valeurs intermédiaires représentent le degré d'incertitude. Ces cartes peuvent être lues comme des cartes binaires de présence et absence, mais tenant compte des incertitudes dans les prédictions. Les cartes démontrent, pour la mésange à tête noire, un déplacement rapide vers le nord, surtout dans l'Ouest de la province, et plus lent vers l'Est, soit les vers régions de la péninsule gaspésienne et la Côte-Nord (AUC = 0.982). Pour les gros-becs errants, les résultats démontrent un déplacement de l'aire de répartition vers le nord (AUC = 0.959). La limite nordique augmente pour les scénarios 2.6, 4.5, 6.0 et augmente fortement pour le scénario 8.5. Pour les scénarios 6.0 et 8.5, les résultats démontrent un déplacement de la limite méridionale vers le nord, plus particulièrement dans le sud-est du Québec et dans la péninsule gaspésienne. Toutefois, l'incertitude pour ces deux secteurs laisse sous-entendre que certains modèles prédisent une présence et d'autres, majoritaires, prédisent une absence.

Figure 7 : Distribution actuelle et prévue de la mésange à tête noire pour 2070 selon les quatre scénarios RCP.

Figure 7 : Distribution actuelle et prévue de la mésange à tête noire pour 2070 selon les quatre scénarios RCP.

Figure 8 : Distribution actuelle et prévue du gros-bec errant pour 2070 selon les quatre scénarios RCP.

Figure 8 : Distribution actuelle et prévue du gros-bec errant pour 2070 selon les quatre scénarios RCP.

Discussion

28Les résultats de l’analyse de redondance canonique (RDA) démontrent que la répartition spatiale des 39 espèces d’oiseaux de la forêt boréale échantillonnés est liée en grande partie à des facteurs bioclimatiques, notamment les températures maximale et minimale, ainsi que les précipitations du trimestre le plus arrosé. Alors que l’analyse initiale inclut plusieurs variables, telles les perturbations anthropiques, la couverture forestière et les milieux humides, seules quatre variables bioclimatiques, l'élévation et la présence de milieux humides sont responsables de la distribution spatiale d’oiseaux. Les variables bioclimatiques expliquent à 53% la distribution spatiale des espèces aviaires, les variables ressources, 5%, et les variables anthropiques, moins de 1%. Ainsi, il est possible d’affirmer que, pour nos données, les perturbations anthropiques semblent n’avoir aucun lien avec la distribution spatiale des oiseaux de la forêt boréale. Toutefois, il importe de rester vigilant avec ces résultats, puisque les observations d’oiseaux sont habituellement prises en bordure des routes ou en milieu anthropisé. Ainsi, nous ne bénéficions pas ou très peu d’observations pour le nord du Québec, ni pour les endroits situés loin des villes et villages. Il serait plus avisé de trouver de meilleures variables représentant l’impact anthropique, notamment des données quant à l’utilisation du sol et aux occupations humaines. De plus, il faut rester vigilant à la sensibilité des analyses à la résolution spatiale des données. En effet, l’échantillonnage aléatoire donne des sites aux propriétés différentes selon la résolution des données. Ces résultats permettent toutefois de considérer les variables bioclimatiques comme un important moteur influençant la répartition spatiale d’oiseaux de la forêt boréale québécoise.

29Ainsi, pour la seconde étape de ce projet, soit la modélisation de distribution spatiale d’espèces d’oiseaux, seules les variables bioclimatiques et l'élévation sont utilisées. Toutefois, des régressions linéaires multiples calculées pour chacune des espèces de façon individuelle permettent d’identifier des variables différentes selon l’espèce. Certaines espèces, comme le viréo aux yeux rouges (REVI), a un R2adj de 0.78 pour quatre variables (elevation, bio11, bio16 et bio7), alors que d’autres, tel le roitelet à couronne rubis (RCKI) a un R2adj de 0.26, avec une seule variable significative (bio11). Les résultats des modèles appuient l’idée que les changements climatiques favoriseraient un changement dans la répartition spatiale des espèces d’oiseaux de la forêt boréale. Cette idée est en accord avec le paradoxe de biodiversité nordique, qui consiste à affirmer que bien que les changements climatiques seront une grande cause d’extinction d’espèces, les régions boréales, elles, subiront un accroissement dans la richesse des espèces et dans la biodiversité (Berteaux et al., 2010). Les résultats d’une étude similaire démontrent que pour 80% des scénarios envisagés, l’ensemble de l’Amérique du Nord devrait subir une perte nette d’espèces animales de 11% sous les scénarios climatiques B1 (conservateurs) ainsi qu’une perte d’espèces animales de 17% sous les scénarios A2 (moyennes à hautes émissions de gaz à effet de serre) (Lawler et al. 2009).

30Les modèles de distribution d’espèces tels que produits dans le cadre de cette recherche sont statiques et reposent sur un postulat de pseudo-équilibre, tel que défini par Guisan et Zimmermann (2000). En effet, toute modélisation statique repose sur la prémisse que les patrons de distribution des espèces modélisées sont en équilibre ou pseudo-équilibre avec l’environnement, puisque les modèles statiques ne permettent pas de gérer les conditions de déséquilibre ou d’équilibre dynamique que démontrent les écosystèmes réels. Les approches de modélisation de type dynamique et stochastique issues des théories des systèmes complexes, telle la modélisation basée sur l’individu ou l’automate cellulaire, permettent de pallier ces problèmes, puisqu’elles reposent sur une approche ascendante (bottom-up), modélisant à l’échelle de l’individu et permettant les rétroactions et la non-linéarité, alors que les approches mobilisées ici reposent sur une approche descendante (top-down), où les aires de distribution d’espèces sont identifiées par des statistiques et par des variables environnementales et non par les traits comportementaux et préférences de ces espèces.

31La méthode dressée ici, consistant en une combinaison d’analyses multivariées et de modélisation de distribution d’espèces pourrait aisément être adaptée à d’autres espèces et/ou à d’autres régions d’étude. La première étape, incluant l’analyse de redondance canonique, la suppression de colinéarité, la sélection bidirectionnelle puis le partitionnement de variation, permet de sélectionner les variables explicatives pertinentes et de les inclure dans la seconde étape, plutôt que de produire des modèles incluant plusieurs variables inutiles.

Conclusion

32Les analyses multivariées permettent d’identifier les variables explicatives responsables de la distribution spatiale de 39 espèces d’oiseaux de la forêt boréale du Québec. Les variables identifiées comme étant des variables clés comprennent l’étendue annuelle des températures, les températures moyennes du trimestre le plus frais, la saisonnalité des précipitations, les précipitations du trimestre le plus arrosé et l’élévation. La variable partageant le rapport le plus fort avec les deux espèces modélisées est bio11, ou les températures moyennes du trimestre le plus frais. Ceci s’explique aisément par le fait que les deux espèces modélisées sont des espèces boréales. Les résultats des analyses multivariées démontrent que les espèces d’oiseaux sélectionnées pour l’étude ont une distribution spatiale fortement corrélée au climat. Ainsi, comme le démontrent les résultats des modèles bioclimatiques, les changements climatiques risquent-ils d’induire des changements dans la répartition spatiale de ces espèces. Des modèles mécanistes et stochastiques prenant en considération la fragmentation forestière et les perturbations naturelles et anthropiques seraient en mesure de rendre justice à la complexité inhérente aux écosystèmes forestiers et, ainsi, de prendre en considération les interactions et rétroactions possibles avec d’autres espèces ou avec le milieu.

Haut de page

Bibliographie

Bertaux, D., de Blois, S., Angers, J.-F., et al., 2010, "The CC-Bio Project: Studying the Effects of Climate Change on Quebec Biodiversity", Diversity, vol. 2, 1181-1204.

Breiman, L., 2001, "Random forests", Machine Learning, vol. 45, 15-32.

Borcard, D., P. Legendre, P., Drapeau, P., 1992, "Partialling out the Spatial Component of Ecological Variation", Ecology, vol. 73, No.3, 1045-1055

Canada's Forest Inventory (CanFI), 2001, "Canada's Forest Inventory for the boreal eco-region - percent forested", Quebec. Jeu de données en ligne. Url : http://databasin.org/datasets/00b4cf27796d4a32be21df7716103e3f

eBird, 2013, "Basic Dataset", Version: ebd_CA-QC_198701_201412_relNov-2013, Cornell Lab of Ornithology, Ithaca, New York. Nov 2013, Jeu de données en ligne.

ESRI, 2012, "ArcGIS Desktop: Release 10.1", Redlands, CA, Environmental Systems Research Institute. Logiciel

Forister, M.L., McCall, A. C., Sanders, N.J. et al., 2010, "Compounded Effects of Climate Change and Habitat Alteration Shift Patterns of Butterfly Diversity", Proceedings of the National Academy of Sciences of the United States of America, vol. 107, No.5, 2088–92.

Friedman, J. H, 1991, "Multivariate Adaptive Regression Splines", The Annals of Statistics, vol. 19 No.1, 1-67

Georges, D., Thuiller, W., Engler, R., 2014, "BIOMOD_EnsembleModeling {biomod2}: Create and evaluate an ensemble set of models and predictions", Librairie R en Ligne. Url : http://www.inside-r.org/packages/cran/biomod2/docs/BIOMOD_EnsembleModeling

Global Forest Watch Canada, 2009, "Cumulative Anthropogenic Access, Eastern Canada". Jeu de données en ligne. Url : http://www.globalforestwatch.ca/data/20100913B

Gouvernement du Canada, Ressources Naturelles Canada, Centre d'information topographique, 2000, "Données numériques d'élévation du Canada", Sherbooke, Québec, Canada. Jeu de données en ligne. Url : http://www.geobase.ca/geobase/en/data/cded/description.html

Guisan, A., Zimmermann, N.E., 2000, "Predictive habitat distribution models in ecology", Ecological Modelling, vol. 135, 147-186.

Guisan, A., Thuiller, W., 2005, "Predicting species distribution: offering more than simple habitat models", Ecology letters, vol. 8, 993-1009

Hastie, T.,Tibshirani, R., Buja, A., 1994, "Flexible Discriminant Analysis by Optimal Scoring", Journal of the American Statistical Association, vol. 89, No.428, 1255-1270.

Hijmans, R.J., Cameron, S.E., Parra, J.L., Jones, P.G. and Jarvis, A., 2005, "Very high resolution interpolated climate surfaces for global land areas", International Journal of Climatology, vol. 25, 1965-1978.

Hurlbert, A.H., and Liang, Z, 2012, "Spatiotemporal Variation in Avian Migration Phenology: Citizen Science Reveals Effects of Climate Change", PloS One vol. 7 No.2, 1-11

IPCC, 2007, "Climate Change 2007 : An Assessment of the Intergovernmental Panel on Climate Change", Edited by R K Pachauri and A Reisinger, Change vol. 446, 12–17.

IPCC, 2013, "Summary for Policymakers", Climate Change 2013: The Physical Science Basis. Contribution of Working Group I to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change, Cambridge University Press, Cambridge, United Kingdom and New York, NY, USA.

Oksanen, Y. F. Blanchet, G., Kindt, et al., 2013, "vegan: Community Ecology Package", Librairie R en ligne. Url : http://CRAN.R-project.org/package=vegan.

Lawler, Joshua J, Sarah L Shafer et al., 2009, "Projected Climate-Induced Faunal Change in the Western Hemisphere", Ecology vol. 90 No.3, 588–97.

Legendre, P., Gallagher, E.D., 2001, "Ecologically meaningful transformations for ordination of species data", Oecologia vol. 129, 271-280

Legendre, P, Legendre, L., 2012, "Numerical Ecology, Third English Edition", Elsevier, Developments in Environmental Modelling, 24.

Lescop-Sinclair, K., and Payette, S., 1995, "Recent Advance of the Arctic Treeline Along the Eastern Coast of Hudson Bay", Journal of Ecology, vol. 83, No.6, 929–936.

McCullagh, P., Nelder, J.A., 1983, "Generalized linear models", Monographs on Statistics and Applied Probability, London, England. Chapman and Hall.

MERN, 2013, "Bref Portrait de La Forêt Boréale Au Québec". Url : http://www.MERN.gouv.qc.ca/presse/feux-grands.jsp.

Nelson, T., Coops, N., Wulder et al., 2014, "Predicting Climate Change Impacts to the Canadian Boreal Forest", Diversity, vol. 6, 133–157.

Parmesan, C., 2006, "Ecological and Evolutionary Responses to Recent Climate Change", Annual Review of Ecology, Evolution, and Systematics, vol. 37, No.1, 637–669.

Phillips, S.J., Dudik, M., Shapire, R.E, 2004, "A Maximum Entropy Approach to Species Distribution Modeling", Proceeding of the Twenty-First International Conference on Machine Learning, 655-662

R Core Team, 2014, "R: A language and environment for statistical computing", R Foundation for Statistical Computing, Vienna, Austria. Url : http://www.R-project.org/.

Royston, P, 1982, "An extension of Shapiro and Wilk’s W test for normality to large samples", Applied Statistics, vol. 31, 176-180.

Ter Braak, C. J. F., 1986, "Canonical Correspondance Analysis: A New Eigenvector Technique for Multivariate Direct Gradient Analysis", Ecology, vol. 67, 1167-1179

Thuiller, W., Georges, D., and Engler, R, 2013, "biomod2: Ensemble platform for species distribution modeling". Librairie R en ligne. Url : http://CRAN.Rproject.org/package=biomod2

Virkkala, R., Heikkinen, R.K., Leikola, N., and Luoto, M, 2008, "Projected Large-Scale Range Reductions of Northern-Boreal Land Bird Species Due to Climate Change", Biological Conservation, vol. 141, No.5, 1343–1353.

Virkkala, R., Marmion, M., Heikkinen, R.K. et al., 2010, "Predicting Range Shifts of Northern Bird Species: Influence of Modelling Technique and Topography", Acta Oecologica, vol. 36, No.3, 269–281.

Walther, G.-R., Post, E., Convey, P et al., 2002, "Ecological responses to recent climate change", Nature, vol. 416, 389-395.

Pearson, RG, and W Thuiller et al., 2006, "Model‐based Uncertainty in Species Range Prediction.", Journal of biogeography, vol. 33, No.10, 1704-1711.

Vuuren, Detlef P., Jae Edmonds, Mikiko Kainuma, et al., 2011, "The Representative Concentration Pathways: An Overview", Climatic Change, vol. 109, No.1 and 2, 5–31.

Haut de page

Annexe

Codage des variables bioclimatiques (Hijmans et al., 2005)

BIO1 = Température moyenne annuelle

BIO2 = Étendue diurne moyenne (Moyenne mensuelle de (tmax – tmin)

BIO3 = Isothermalité (BIO2/BIO7) (* 100)

BIO4 = Saisonnalité des températures (écart-type * 100)

BIO5 = Température maximale du mois le plus chaud

BIO6 = Température minimale du mois le plus frais

BIO7 = Étendue annuelle des températures (BIO5-BIO6)

BIO8 = Température moyenne du trimestre le plus arrosé

BIO9 = Température moyenne du trimestre le plus sec

BIO10 = Température moyenne du trimestre le plus chaud

BIO11 = Température moyenne du trimestre le plus frais

BIO12 = Précipitations annuelles

BIO13 = Précipitations du mois le plus arrosé

BIO14 = Précipitations du mois le plus sec

BIO15 = Saisonnalité des précipitations (coefficient de variation)

BIO16 = Précipitations du trimestre le plus arrosé

BIO17 = Précipitations du trimestre le plus sec

BIO18 = Précipitations du trimestre le plus chaud

BIO19 = Précipitations du trimestre le plus frais

Codage des autres variables

elevation = élévation en mètres

pct_for = pourcentage de forêt par cellule

pct_wet = pourcentage de milieu humide par cellule

linear = perturbations de type linéaires

Haut de page

Table des illustrations

Titre Figure 1 : Sites échantillonnés
URL http://cybergeo.revues.org/docannexe/image/26969/img-1.png
Fichier image/png, 171k
Titre Figure 2 : Tendances dans l’évolution des concentrations de gaz à effet de serre selon les 4 profils d’évolution de la concentration. De gauche à droite, on retrouve : le dioxyde de carbone (CO2) le méthane (CH4), et les oxydes nitreux (N2O) (Vuuren et al. 2011)
URL http://cybergeo.revues.org/docannexe/image/26969/img-2.png
Fichier image/png, 33k
URL http://cybergeo.revues.org/docannexe/image/26969/img-3.png
Fichier image/png, 1,1k
Titre Figure 3 : Diagramme de Venn pour partitionnement de la variation de Y en fonction de deux groupes de variables explicatives
URL http://cybergeo.revues.org/docannexe/image/26969/img-4.png
Fichier image/png, 4,7k
Titre Tableau 1 : Algorithmes de classification et de régression utilisés pour la génération des modèles de distribution d’espèces
URL http://cybergeo.revues.org/docannexe/image/26969/img-5.png
Fichier image/png, 21k
Titre Schéma  : Adapté de Berteaux et al., 2010. Représentation simplifiée de la méthode de modélisation bioclimatique utilisée.
Légende Les chiffres entre parenthèses représentent les étapes dans l’exécution du code. Les parallélogrammes représentent des entrées/sorties, les rectangles, des processus et le cylindre bleu, une base de données. Le rectangle rouge représente le cœur du modèle, soit le modèle de niche.
URL http://cybergeo.revues.org/docannexe/image/26969/img-6.png
Fichier image/png, 321k
Titre Figure 4 : Résultats de la RDA sur les variables non-colinéaires (p<0.005, R2adj= 0.58)
URL http://cybergeo.revues.org/docannexe/image/26969/img-7.png
Fichier image/png, 29k
Titre Figure 5 : Résultats de la RDA sur les variables sélectionnées (p<0.005, R2adj= 0.59)
URL http://cybergeo.revues.org/docannexe/image/26969/img-8.png
Fichier image/png, 27k
Titre Tableau 2 : Résultat du test de permutation sur les prédicteurs
URL http://cybergeo.revues.org/docannexe/image/26969/img-9.jpg
Fichier image/jpeg, 52k
Titre Tableau 3 : Résultats du test de permutation sur les axes canoniques
URL http://cybergeo.revues.org/docannexe/image/26969/img-10.jpg
Fichier image/jpeg, 48k
Titre Figure 6 : Partitionnement de la variation des données des 39 espèces d’oiseaux (log) par deux groupes de variables explicatives : les variables bioclimatiques et les variables ressource
URL http://cybergeo.revues.org/docannexe/image/26969/img-11.png
Fichier image/png, 12k
Titre Figure 7 : Distribution actuelle et prévue de la mésange à tête noire pour 2070 selon les quatre scénarios RCP.
URL http://cybergeo.revues.org/docannexe/image/26969/img-12.png
Fichier image/png, 585k
Titre Figure 8 : Distribution actuelle et prévue du gros-bec errant pour 2070 selon les quatre scénarios RCP.
URL http://cybergeo.revues.org/docannexe/image/26969/img-13.png
Fichier image/png, 658k
Haut de page

Pour citer cet article

Référence électronique

Jonathan Gaudreau, Liliana Perez et Pierre Legendre, « Identification des variables expliquant la distribution spatiale d’oiseaux de la forêt boréale et modélisation de tendances futures : une approche multivariée », Cybergeo : European Journal of Geography [En ligne], Systèmes, Modélisation, Géostatistiques, document 722, mis en ligne le 11 mai 2015, consulté le 23 août 2017. URL : http://cybergeo.revues.org/26969 ; DOI : 10.4000/cybergeo.26969

Haut de page

Auteurs

Jonathan Gaudreau

Département de Géographie, Université de Montréal, Canada
Jonathan.gaudreau@umontreal.ca

Liliana Perez

Département de Géographie, Université de Montréal, Canada
l.perez@umontreal.ca

Pierre Legendre

Département de Sciences Biologiques, Université de Montréal, Canada
Pierre.Legendre@umontreal.ca

Haut de page

Droits d’auteur

© CNRS-UMR Géographie-cités 8504

Haut de page