Omnidirectional and omnifunctional connectivity analyses with a diverse species pool

Connectivity among habitat patches in both natural and disturbed landscapes needs to be accounted for in conservation planning for biodiversity maintenance. Yet methods to assess connectivity are often limited, because simulating the dispersal of many species is computationally prohibitive, and current simulations make simplifying assumptions about movement that are potentially erroneous. Here we show how these limits can be circumvented and propose a novel framework for the assessment of omnifunctional and omnidirectional connectivity in a 28000 km2 area in the Laurentian region of Québec, Canada. Our approach relies on (i) the use of Omniscape, an improved version of Circuitscape which allows omnidirectional simulations that better emulate animal movement and (ii) the synthesis of large volume of species-level dispersal simulations through a posteriori clustering of the current intensity. Our analysis reveals that the movement of 93 species evaluated can be clustered into three functional dispersal guilds, corresponding to mostly aquatic species, terrestrial species able to use aquatic environments, and strictly terrestrial species. These functional guilds do not share connectivity hotspots, suggesting that corridor planning would need to account for the multiplicity of dispersal strategies. Although this approach requires a large volume of computing resources, it provides richer information on which landscape features are critical to maintain or need to be regenerated for broader biodiversity maintenance goals.

intensity. Our analysis reveals that the movement of 93 species evaluated can be clustered into 33 three functional dispersal guilds, corresponding to mostly aquatic species, terrestrial species able 34 to use aquatic environments, and strictly terrestrial species. These functional guilds do not share 35 Introduction 41 7 weighted means of different landscape features to assess link probabilities, and two different 105 summation types for nodes in the network. These studies have demonstrated that it is possible to 106 create connectivity analyses for a large number of species, which may be ecologically more 107 relevant. Indeed the results could be meaningful for decision-makers, helping them to target new 108 protection areas, restore others, and inform more sustainable strategies in city development 109 projects (Ersoy et al., 2019). Despite these successes, they remain limited by three key aspects. 110 First, even a pool of 20 species cannot represent the entire functional diversity of realistic 111 communities. Second, selected species pools continue to under value the contribution of semi-112 aquatic species and their habitats to landscape connectivity. Furthermore, freshwater biodiversity 113 is considered the most imperiled globally (Hendriks, 2016), suggesting that the need to include 114 aquatic and semi-aquatic species is urgent. Finally, because species clustering is done before the 115 simulations, this assumes that species connectivity depends on habitat requirements alone and 116 neglects differences that could emerge due to landscape configuration. 117 In this study we present a general framework for multi-species connectivity analysis, built 118 around the needs to (i) integrate terrestrial and aquatic ecosystems (Herbert et al. for the way species are predicted to move across the landscape by looking for similarities in their 125 8 dispersal simulations, rather than in the cost matrix of these species. We argue that this new way 126 of interpreting and post-processing connectivity simulations will result in a more integrative 127 view of landscape connectivity to favour biodiversity maintenance in the long-term. 128 129 9

131
We assessed our new connectivity approach using a 27 994 km 2 area located on the north shore 132 of the St. Lawrence River, near the city of Montréal, Québec, Canada (Fig. 1). Seven major 133 rivers are included in this region, and the area sits on two strikingly different geological 134 provinces, the Saint Lawrence lowlands, dominated by agriculture and urban areas, and the 135 Canadian Shield, covered largely by pristine temperate forests with a very high density of lakes 136 and rivers (19.8% ). This high aquatic surface supports the need to account for more semi-137 aquatic species in our simulations. This region was chosen because it shows a strong gradient of 138 anthropogenic disturbance, with heavy urbanization and agricultural development to the south, 139 major highways on the north-south axis, many roads, large forested areas particularly on the 140 Shield where there is an active forestry industry, some major National Parks (Mont Tremblant,  141 Oka, la Mauricie), and several small cities that service a large cottage industry. Moreover, 142 biological surveys conducted by conservation agencies and different government bodies resulted 143 in the creation of a curated list of species, which we integrated in this analysis. 144

145
In order to select the species for this study, two lists of recorded species from the main National 146 Parks were used (SEPAQa, 2019; SEPAQb, 2019). We kept 45 mammals, all amphibians (16) 147 and 9 reptile species, plus 15 avian species with a conservation status. We added five more avian 148 species, two amphibians and one reptile, which have been recorded in GBIF and Ebird database 149 since 2000, for a final list containing 93 species (see Appendix 2). 150 Collecting species habitat preference data 151 From the list of focal species, we got information on habitat preference and the way they move 152 across the landscape. Since these data are very heterogeneous, we gathered them from several The complete cost matrix for all species is given in APPENDIX 2. 157 Creating the resistance maps

158
Resistance maps are a mix between energetic landscape and habitat quality modelling. They 159 represent the energetic and ecological cost to cross a certain plot as a function of its 160 configuration (topography, human structures) and composition (land cover) ( herbaceous wetland, etc.). To transform it into a resistance layer, we integrated the resistance 170 score of each species, to each type of land cover in the layer's data (Appendix 2). This is the 171 basic landscape features which influence species dispersal. In order to assess landscape 172 complexity to its fullest extent, we added the effect of slope, buildings presence, waterbodies, In order to create a resistance map by summing the layers for each species, we converted vector 180 resistance layers to raster layers. We first chose to use square cells at a resolution of 2 meters 2 , 181 which represents the smallest disturbance width considered in this study, hiking trails. Then, for 182 each individual species, we summed the resistance layers corresponding to their canvas, which 183 produced the 93 species-specific resistance maps. The resistance score of each cell represents the 184 sum of all landscape features, according to the species dispersal behaviour canvas. Analyzing landscape connectivity using Omniscape 2008). We used three different radiuses (30, 60 or 120 blocks, each block being a square parcel 195 of 180m 2 ), and averaged the quantiles of the cumulative current map for all. This resulted in 196 connectivity scores ranging from close to 0 (specifically the inverse of the number of simulated 197 pixels) to 1. The choice to examine quantiles as opposed to raw cumulative current values was 198 13 driven by the fact that although different species will disperse at different rates (and therefore 199 will have different total cumulative currents), quantiles provide information as to their relative 200 use of the landscape. This would therefore facilitate the search for functional dispersal guilds. 201 Note that for some applications, analyzing the normalized current (i.e. the cumulative current 202 divided by the potential flow, as in e.g. McRae (2008), which is to say measuring how surprising 203 the observed current is compared to a null expectation) is preferable, notably for corridor 204 mapping or the identification of dispersal bottlenecks. In this analysis, we focus on how 205 organisms could disperse in the landscape, and therefore use the cumulative maps. 206 All simulations were performed on the beluga supercomputer, operated by Calcul Québec; the 207 project required an estimated total of 3 core-years to complete, being roughly equivalent to 3 208 years of time on a single-core machine. 209 A posteriori functional dispersal guilds 210 After all the simulations were complete, we aggregated the maps using fuzzy C-means (C. 211 Bezdek, 1981; Dunn, 1973). This approach allows to have species contribute to more than a 212 single cluster (for example, species that are terrestrial but have an affinity for water will have 213 landscape use that borrows from strictly aquatic and strictly terrestrial species), as captured by 214 the fuzziness parameter. Following community-established best practices, we set this parameter 215 to 2. Specifically, we applied fuzzy C-means by transforming all the connectivity scores into a s 216 by p matrix, with s the number of species and p the number of pixels with non-null connectivity 217 values, then by calculating the s by s Euclidean distance matrix between all pairs of species. We 218 extracted the centers of every cluster, and the species with the highest weight was the most 219 representative species for the dispersal guild. It may therefore be used as an indicator species in 220 the future. In addition, the weight of each species for each cluster can be used to measure species 221 specificity; we used PDI (Poisot et al., 2012) as a specificity estimator, which has the desirable 222 property of returning values that are lower than 0.5 for generalist species, and values above 0.5 223 for specialist species. terrestrial habitats are more difficult to reconcile (fig 4b, 4c). This includes, trivially, the St 272 Lawrence river itself, but also areas surrounding large water bodies such as the Lac Taureau to 273 the North. It should be noted that although both the coefficient of variation and Pielou's evenness 274 provide a similar information (namely, the homogeneity in the use of every pixel by the three 275 FDGs), the later is a more appropriate measure. Although widely used, the coefficient of 276 variation performs better on ratio scales, whereas the quantiles of dispersal, being bounded, are 277 expressed on an interval scale. 278 Best-connected pixels

279
In fig. 5a, we represent the 17% of pixels with the highest average connectivity, which can be 280 thought of as forming the basis for a dispersal network across the region, independent of guilds. 281 18 This information should be contrasted to fig. 4b, in which we show the spatial distribution of the 282 best 17% of pixels for each of the different FDGs and superimposed them. Because the guilds 283 are optimized to maximize the differences between them, their most connected areas rarely 284 overlap. As a result, protecting the top 17% of pixels for all guilds combined would cover up to 285 three times more surface than the best 17% pixels on average. Furthermore, note that very few 286 aquatic habitats would be protected using the average classification as a target.

291
We have simulated the dispersal of 93 species across a 220 by 230 km area in the Laurentian 292 region of Québec. We then classified the multi-model average of connectivity maps using Fuzzy 293 C-means, to provide a regional average (Fig. 4a). Based on this result, we identified areas of 294 higher uncertainty (Fig. 4b, 4c) as well as candidate indicator species for each of the three 295 Functional Dispersal Guilds (FDGs). We found that connectivity is lower in the St Lawrence 296 Lowlands, which is more impacted by anthropogenic activities, and generally better in the 297 forested areas of the Canadian Shield. Although this result was expected, our approach was the 298 first to assess connectivity across such a large heterogeneous landscape, by accounting for the 299 full taxonomic diversity of both terrestrial and semi-aquatic species. Furthermore, these species 300 belong to different guilds representing three emergent dispersal landscapes that coexist within 301 the study region. Conservation actions that may neglect the protection of a specific guild due to 302 poor representativity of indicator species could have overall negative consequences on the 303 biodiversity maintenance of the entire region. Indeed, semi-aquatic, which form a guild on their 304 own, are often neglected. 305 We have used fuzzy C-means to cluster the species into FDGs, as it represents a more objective 306 approach in which species are able to contribute relatively among dispersal guilds. Indeed, we 307 identified about 16% of species classified as generalists that contribute in similar proportions to 308 all FDGs. Under a non-fuzzy scheme, this information would be lost as species would be forced 309 into one guild only. This being said, as we show in fig 6, the results of the fuzzy C-means 310 20 algorithm are directly comparable with a hierarchical clustering using Ward's agglomeration on 311 the pre-squared distance matrix. Although simple hierarchical clustering would have led to the 312 same dispersal guilds, the rich information of the movement of generalists species among guilds 313 could only be achieved using the fuzzy C-means approach. 314 Clustering returns groups of species based on predicted dispersal, and therefore represents a 315 process to select indicator species within each guild. By contrast, the a priori selection of 316 indicator species risks either underrepresenting or over representing different dispersal modes. 317 This is a particularly concerning issue, since as we illustrate in fig 4, the identification of the 318 most connected pixels varies greatly depending on whether the FDGs, or the average 319 connectivity, are used. Knowing the importance of the choice of representative species, the more 320 ecologically relevant a posteriori clustering should be favored, especially given that this amount 321 of simulations is computationally feasible. Indeed, among the 93 species used in this study, 322 several that would be considered "terrestrial" only, were found to benefit from freshwater 323 habitats for travel and therefore clustered differently. For example, this was notably the case for 324 the clustering of moose (Alces alces, Linnaeus 1758) with the common garter snake 325 (Thamnophis sirtalis, Linnaeus 1758) and many avian species, among others. This 326 counterintuitive association demonstrates that with regards to landscape connectivity, our 327 preconceived ideas of species groupings do not apply. 328 Despite the large diversity of species and their habitat preferences, fuzzy C-means identified an 329 optimal grouping of three functional dispersal guilds only. Surprisingly this is a lower number of 330 classes than most multi-species simulations have been using, and so it calls into question whether 331 21 using more species would have led to pseudo-replication, or to over-representing the needs of 332 some groups of organisms. Previous methods to suggest surrogate species for connectivity 333 modeling also resulted in a larger number of species (e.g. between 5 to 7 for Meurant et al. 334 2018). We argue that the number of species guilds should not be determined a priori, since these 335 classifications would rely only on between-species differences in habitat preferences, and would 336 not be able to account for how these differences interact with the spatial configuration and 337 habitat diversity in the landscape. Therefore we suggest our approach be used with the largest 338 number of species available, particularly in biophysically complex and species rich landscapes. 339 Once initial clustering is made, the entire species list may no longer be required for subsequent 340 simulations carried out for the same region. In a context such as a predicting the effect of land-341 use or climate change on connectivity, which would involve many thousands of simulations, one 342 strategy would be to identify the most representative species for the optimal number of FDGs, 343 and then use these indicator species in the simulations. While this results in the loss of some 344 information, this offers a reasonable compromise between the effort to identify relevant species 345 in initial computationally intensive analysis, and the effort to perform projections of future 346 connectivity maps. Additional schemes to pick a subset of species to work on should be 347 evaluated on a case by case basis and can involve picking a mix of specialist and generalist 348 species. 349 Here we provide, to our knowledge, the most robust and fully integrated terrestrial and aquatic 350 multi-dispersal connectivity assessment at a regional scale. We used a broad suite of diverse taxa 351 in a heterogeneous landscape structured by a number of different anthropogenic disturbances 352 22 across two major geological provinces. The integration of the multiple dispersal pathways 353 provides a robust connectivity conservation plan and our approach offers a road map for future 354 connectivity assessments. Given the ongoing biodiversity crisis, identifying and implementing 355 corridors to maintain habitat connectivity for conservation for as many species as possible is 356 increasingly urgent. Yet, current conservation strategies often neglect the diversity of ways in 357 which species actually use the landscape, since they focus largely on terrestrial species. For 358 example, as a worst case scenario, habitat patches crucial to the dispersal of a specific functional 359 group would not be adequately protected; this connectivity loss for an entire subgroup of species 360 could have a cascading effect with the potential to destabilize the entire ecosystem. 361   features got scored between 0 and 100 for their resistance (1). Then, the layers were converted in 381 a raster format to be summed, which created one resistance map by species (2). Following this 382 step, the maps were analysed using Omniscape, which created connectivity maps (3). We 383 assigned to each pixel its quantile in the observed map of the cumulative current for the species 384 (4) and conducted a clustering fuzzy C-means analysis on the unfolded 93 connectivity map to 385 get the a posteriori dispersal guilds (5). After the clustering, we extracted the weighted centroids 386 of each cluster (6) and measured the additional indices (evenness, variation, location of the most 387 connected pixels) (7). 388