How does dispersal shape the genetic patterns of animal populations in European cities? A simulation approach

Abstract


Introduction
Most of humankind currently lives in cities and urban areas are predicted to cover three times the area they covered in 2000 by 2030 (Seto et al., 2012), and four times by 2050 (Angel et al., 2011).Urbanization is an important component of the anthropogenic pressures triggering the erosion and spatial redistribution of biodiversity (Diaz et al., 2019).Indeed, it favors the spread of invasive species, concentrates pollution, and urban lifestyles are important drivers of natural resource over-exploitation and greenhouse gas emissions (McDonald et al., 2020).Last but not least, it is a direct cause of habitat destruction and fragmentation (Beninde et al., 2015).As a result, many empirical studies have evidenced negative effects of urbanization on biodiversity (Aronson et al., 2014;Piano et al., 2020).However, this relationship is complex and variable across taxa, biological organization levels (ecosystems, species, genes), and across cities (Fidino et al., 2020).Urban ecology research is therefore needed and particularly crucial if humanity is willing to limit its impact on biodiversity and to maintain the ecosystem services it provides (Verrelli et al., 2022).
Urban ecology studies have shown that species are not all affected in the same way by urbanization (Aronson et al., 2014;Blair, 1996;Fanelli et al., 2022;Fidino et al., 2020).While some species are mostly, if not only, present in cities because they are specialist of anthropized environments (urban adapters), others cannot survive in these areas (urban avoiders).Some more ubiquitous species are present in both urban and non-urban areas.These urban tolerant species are reliable indicators of urban influences on population dynamics along rural-to-urban gradients.They may also be the most affected by the impact of urban planning on environmental conditions across such gradients.Species-specific responses to urbanization not only affect species diversity patterns, but also explain the variability of genetic patterns observed at the intra-specific level along rural-to-urban gradients.Although rapid genetic adaptations to urbanization have been observed in several species (Santangelo et al., 2022), urbanization also shapes neutral genetic patterns in a significant manner (Fusco et al., 2021;Miles et al., 2019).For instance, urban populations of human pests can exhibit higher levels of genetic diversity than non-urban ones (Miles et al., 2018).In contrast, Khimoun et al. (2020) and Schoville et al. (2013) did not detect any significant difference in genetic diversity nor any patterns of isolation by distance pattern when studying ant and butterfly populations, respectively, in urban and non-urban settings.Similarly, the relationship between amphibian genetic diversity or differentiation and the degree of urbanization of several North American cities was not significant in the study by Schmidt and Garroway (2021).On the contrary, Delaney et al. (2010) showed that urbanization decreased genetic diversity and increased genetic differentiation in three lizard species and one bird species, mainly due to higher road density in urban areas.Likewise, Stillfried et al. (2017) showed that urban populations of wild boars exhibited lower genetic diversity levels than suburban ones in Berlin.
The variability of the neutral genetic patterns observed in urban settings stems from the demographic dynamics determining the intensity of both genetic drift and gene flow (Frankham et al., 2004;Miles et al., 2019;Munshi-South and Richardson, 2020).On the one hand, genetic drift can lead to allele loss, especially in small-sized populations.On the other hand, gene flow following successful dispersal events between populations leads to genetic exchanges.This can compensate for the diversity loss due to drift and decrease the resulting genetic differentiation between populations.Consequently, both the size of urban populations and the permeability of urban environments to individual dispersal drive genetic diversity and genetic differentiation, because they determine genetic drift and gene flow.As such, the genetic patterns described above span the whole range of patterns that are theoretically expected from variations in the respective intensity of drift and gene flow across urban and non-urban areas (Frankham et al., 2004;Hutchison and Templeton, 1999).Under strong genetic drift and low gene flow, a low genetic diversity is expected.Besides, isolation by distance patterns should not be significant, due to the prevailing role of local stochasticity on pairwise genetic differentiation irrespective of pairwise distances among populations.On the contrary, under moderate genetic drift and gene flow, a significant isolation by distance pattern and a higher genetic diversity level are expected.Finally, at intermediate levels, isolation by distance patterns are expected, at least temporarily, until a pairwise distance threshold beyond which drift becomes the prevailing stochastic driver of genetic differentiation.However, because these predicted genetic patterns depend on the respective intensities of both drift and gene flow, disentangling the influence of these two processes is a complex task.
In urban tolerant species, genetic drift could explain most of the variability observed in genetic patterns given that some of these species only form small populations in urban areas (e.g., Lourenco et al. (2017)), thereby exacerbating drift effects, whereas others can maintain large populations (e.g., Miles et al. (2018)).Yet, gene flow, which is mainly driven by the ability of urban tolerant species to disperse across urban fabric, could also be the prevailing driver of genetic responses to urbanization.Accordingly, assessing how dispersal limitations and resulting gene flow reductions shape genetic patterns is crucial for several reasons.First, many biodiversity conservation programs rely on habitat connectivity conservation and restoration for maintaining genetic diversity in urban areas.They assume that dispersal limitation is the main cause of biodiversity loss.Besides, they often focus on Urban Green Spaces (thereafter referred to as UGS) as biodiversity sources, although their spatial location within cities could convert them into sink patches (Lepczyk et al., 2017;Pulliam, 1988;Verrelli et al., 2022), thus compromising the efficiency of habitat and ecological corridor restoration in urban cores.Determining to which point species movements in these areas are sufficient for preventing diversity loss is therefore needed for estimating the potential benefit of such measures.For that purpose, simulations have been commonly recommended in landscape genetics (Balkenhol et al., 2016;Munshi-South and Richardson, 2020), particularly for studying how neutral genetic patterns emerge from the interplay of drift and gene flow processes, independently of adaptive processes.They also proved to be efficient for reproducing empirically observed genetic patterns in urban areas (Rochat et al., 2017).
Considering this, our research objective in this study was to answer the following question: how does dispersal limitation explain the variability of genetic patterns in urban tolerant species?To this end, we simulated neutral genetic patterns resulting from multi-generational gene flow between populations of urban tolerant species located in both UGS and forest areas within and around 325 European cities.Using scenarios introducing variations in the ability of three virtual animal species to disperse across urban fabric, we assessed how this ability affects genetic patterns, independently from any other process.Across scenarios, genetic drift intensity was constant for a given set of populations.The simulated variations in dispersal, for each city, could thus affect genetic patterns by modifying the respective intensities of drift and gene flow.We also compared the connectivity of UGS and forest areas, and contrasted the genetic diversity Paul Savary et al. and differentiation levels observed in these habitats to shed light on how the connectivity and spatial configuration of habitats drive genetic responses to dispersal scenarios in cities.

Material and methods
To answer our research questions, we needed to assess the genetic responses to urbanization of urban tolerant species having different abilities to cross artificial areas, while being equivalent in terms of population density and dispersal abilities to cross favorable areas.We also needed to assess these responses at the level of entire urban areas and to replicate the analyses to ensure that our results were not merely due to the specificity of a single study area.To meet these conditions, we implemented a simulation approach and applied it to 325 European urban areas, hereafter referred to as cities.
For each city, we first modeled the connectivity of habitat patch networks for forest-dwelling species occupying both UGS and forest patches according to three dispersal scenarios, using a graph-based approach.Then, we simulated drift and gene flow processes in populations occupying both types of habitats, and analyzed the resulting genetic pattern at both the withinpopulation (genetic diversity) and between-population (genetic differentiation) levels.We provide details in the following sub-sections.

Habitat connectivity analyses
Land cover data.We used urban land cover data from the Urban Atlas 2018 database of the Copernicus European agency.These land cover data are available for 788 functional urban areas, sensu Organisation for Economic Co-operation and Development (OECD, Moreno-Monroy et al. (2021)), counting at least 50,000 inhabitants in 38 member or partner states of the European Union.These vector land cover data include 27 land cover types at a relatively fine spatial resolution (Minimum Mapping Unit: 0.25 ha inside urban core areas, and 1 ha in surrounding rural areas).They were reclassified into 10 land cover types (Supporting information 1 -Table S1) and rasterized at a spatial resolution of 5 m.The two "Forests" and "UGS" land cover types were considered for delineating habitat patches in the analyzes."Forests" tended to be peripheral (i.e., peri-urban) whereas "UGS" were mostly located within city cores, although in each city, some patches did not conform to this pattern.UGS included public green areas predominantly dedicated to recreational use (e.g., gardens, parks), as well as suburban natural areas that are used and managed as urban parks.We obtained the spatial coordinates of the center of every city from the Open Street Map database using Nominatim and the jsonlite R package (Ooms, 2014).In most cases, it coincided with the city hall, which is commonly used as the center of cities (Lemoy and Caruso, 2020)(Figure 1).
We delineated the cities under study by standardizing their proportion of artificial areas in order to assess the influence of dispersal on genetic patterns in areas having the same degree of urbanization, although differing in the configuration of their urban fabric and natural areas.Delineating study areas of the same spatial extent would have led the differences among cities to mainly reflect the effect of varying densities of artificial areas (i.e., proportion of sealed area in the study area).In contrast, we here assessed species genetic responses in areas which differ in the location and shape of their urban fabric.To that end, we calculated the proportion of artificial areas in disks of increasing radius centered on the city center.These radius ranged from 5 km to 40 km, with steps of 500 m.The target proportion of artificial areas was fixed at 20% ± 1% after preliminary analyses.This threshold allowed us to maximize the number of cities for which a radius encompassing a fixed proportion of artificial areas could be found between 5 and 40 km and to minimize the variance of this radius.Cities were included in our sample when at least 95 % of the delineated disk was covered by the Urban Atlas land cover data (Figures 1A and 1B).When required, we completed the remaining peripheral sectors with Corine Land Cover data, matching their typology with our land cover classification, as indicated in Table S2 (Supporting information 1).Our method selected 325 cities and mainly excluded coastal cities and very small or very large cities for which the proportion of artificial areas was never below or above 20 %, respectively, within the range of considered radii.
Dispersal cost scenarios.We wanted to assess whether genetic responses were affected by the capacity of species to cross the least favorable areas when dispersing from one habitat patch to another.We therefore made three dispersal cost scenarios, consisting of cost values associated with land cover types and representing the cost of species movements across pixels of each land cover type (Table 1).These costs were similar to the costs used by Sahraoui et al. (2017) and Tannier et al. (2016) for modeling habitat connectivity in urban areas for forest species such as rodents (e.g., Mustela putorius) or birds (e.g., Picus canus).They were minimal (1) in habitat areas and higher in land cover types supposed to highly affect forest species movements, such as grasslands (10), agricultural areas (50) or wetlands (100).These costs were constant across scenarios.In contrast, the costs associated with water, roads and artificial land cover types increased from scenario 1 to 3. They were equal to 90, 900 and 9,000 for water and roads, and to 100, 1,000 and 10,000 for artificial areas in scenarios 1, 2 and 3, respectively.This reflects the fact that built-up areas are more difficult to cross than unfavorable yet open areas such as roads or water bodies.A cost scenario based on similar assumptions has been empirically validated by Balbi et al. (2019) in an urban context.The cost variations in the least favorable land cover types for a forest species allowed us to distinguish three virtual species with different dispersal behaviors in urban areas.Nonetheless, the total cost they could endure during dispersal was kept equal for the three virtual species.This means that they had similar absolute capacities to disperse whatever the cost scenario.Yet, the spatial paths they followed when dispersing could vary substantially from one scenario to another depending on the configuration of land cover types.Using these cost scenarios, we thus aimed to isolate the effect of variations in spatial patterns of dispersal on genetic patterns across rural-to-urban gradients.
Landscape graphs.We modeled habitat connectivity with landscape graphs using the Graphab 2.8 software program (Foltete et al., 2021).Landscape graphs represent habitat patch networks as sets of habitat patches (nodes) connected by potential dispersal paths (links) (Urban and Keitt, 2001).We built them using the Urban Atlas land cover data and the three dispersal cost scenarios (one graph per scenario).Each forest or UGS habitat patch above 0.25 ha was a node of the graphs.Although a single type of habitat patch (node) is most often considered in graph-based connectivity analyzes, the special feature of our analyzes was the distinction between forests and UGS, considered as two distinct types of nodes in subsequent analyzes (Savary et al., 2024a).
We computed least-cost paths between every pair of habitat patches.This simple method estimates dispersal paths by finding the ones connecting patches while minimizing the total cost accumulated when crossing pixels with specified cost values.Despite the known limitations of this method (Zeller et al., 2012), it has proved to be relevant in empirical studies measuring Paul Savary et al.  species movements in urban areas (e.g., Balbi et al. (2019)).For each scenario, we obtained a set of potential dispersal paths and the accumulated cost along them, i.e., the cost-distance (Figure 1C).We then created three minimum planar graphs, sensu Fall et al. (2007), in every urban area.In these graphs, links correspond to least-cost paths, connect neighbor patches, and are weighted by the corresponding cost-distances.
We used these graphs to assess the contribution of each habitat patch to the connectivity of the habitat network and identify whether the distribution of forest and UGS patches could drive potential source-sink dynamics.To that purpose, we computed the 'Flux' (F ) connectivity metric for each patch to estimate the amount of habitat that is reachable from the focal patch.It considers the area of the other patches and the cost-distances associated with the shortest path to these other patches on the graph, as follows: with i the focal patch index, j the index of every other patch among the n habitat patches, d ij the cost-distance between patches i and j, and a j the area of each patch j.
α was set such that the probability of covering a cost-distance equivalent to an average path of 5 km (d 5km ) is equal to 0.05, i.e., p(d 5km ) = e −αd 5km = 0.05.This distance can be considered as the maximum dispersal distance of forest species with medium dispersal capacities (Sahraoui et al., 2017).To obtain d 5km , we assessed the relationship between (1) the length in metric units of the least-cost paths not crossing the least favorable areas and (2) the corresponding costdistances, using log-log linear regressions (Tournant et al., 2013).For that purpose, we only considered the set of paths (spatial trajectories and associated cost-distances) that never crossed any pixel of artificial area, road or water, as they are supposed to be representative of the most common species movements.We then used the estimated relationship between the length of these paths and their cost-distances to convert 5 km into cost units and computed the average converted value for each urban area and each scenario.d 5km averaged 20,000 cost-distance units.By using the same d 5km and α values for all three dispersal scenarios in this analysis and for parameterizing the genetic simulations (see below), we ensured that the virtual species we considered had similar absolute capacities to disperse whatever the cost scenario.
Paul Savary et al. 7 Peer Community Journal, Vol. 4 (2024), article e40 https://doi.org/10.24072/pcjournal.407 When computing the F metric, we distinguished cases where patches i and j were respectively either (i) both forest patches (F Forest↔Forest ), (ii) forest and UGS patches (F Forest↔UGS ), (iii) UGS and forest patches (F UGS↔Forest ) or (iv) both UGS patches (F UGS↔UGS ).For each urban area and scenario, forest patches associated with F Forest↔UGS values in the upper quartile of the distribution were considered as "Forest Interface" patches.Similarly, UGS patches associated with F UGS↔Forest values in the upper quartile of the distribution were considered as "UGS Interface" patches.Indeed, these habitat patches are the most important ones for the connectivity between several types of habitats.We expect them to be mostly, though not strictly, located at the periphery of city centers.They could therefore play a significant role in potential sourcesink dynamics.Their connectivity could also drive genetic diversity transfers from less to more anthropized habitat areas, and conversely, e.g., when UGS Interface patches are connected to other UGS located closer to the city center.
Finally, we computed the Equivalent Connectivity (EC ) metric for estimating the connectivity of the whole habitat patch network.This metric represents the area of the unique patch that would provide species with the same amount of reachable habitat as the whole habitat patch network, given its degree of subdivision and the resistance of the matrix (Saura et al., 2011).We used the following formula (see also Supporting information 2 -Figure S1): Given that we distinguished two types of habitat patches, we could estimate the contribution to EC of (i) the forest patches and the connections between them (EC Forest.Forest , both i and j are forest patches), (ii) the UGS patches and the connections between them (EC UGS.UGS , both i and j are UGS patches) and (iii) the connections between forest and UGS patches, weighted by their respective areas (EC Forest.UGS , i is a forest patch and j a UGS patch, or conversely).To make these values comparable among cities, we standardized them by the total area of each city.We also assessed the relative value of the three EC components here considered, standardizing each of them by their sum.

Genetic simulations
Population size and location.The large number of patches in each city prevented us from simulating a population in all habitat patches, due to limitations in computational capacities.Furthermore, we can reasonably assume that habitat patches are not all occupied in actual metapopulations.Nevertheless, we wanted the number of populations to reflect the subdivision of habitat areas.Therefore, for each type of habitat h, the number of populations Npop h,c in city c ranged from 10 to 400 and was proportional to the logarithm of the number of patches of that type in the city c (p h,c ), such that: with min ∀i p h,i and max ∀i p h,i the minimum and maximum number of patches of type h across all cities.
Then, for each city and habitat type, we randomly sampled a number of patches equal to Npop h,c among the habitat patches of type h in the city c, before populating them with individuals.As we wanted the total population in a given city to reflect the total amount of habitat of that city, the total number of individuals Nind c in each city ranged from 500 to 10,000 and was proportional to the logarithm of the total habitat area in each city (using the same formula as above).Log-transformations normalized the distributions of population and individual numbers.Each population (i.e., sampled habitat patch) was then populated with at least 10 individuals, in a way that reflected the area of its patch.We wanted to assign larger populations to large habitat patches while ensuring that the total number of individuals across populations was equal to Nind c .Thus, the numbers of individuals in each population were randomly drawn following a multinomial distribution.Each sampled patch was associated with a probability of being assigned supplementary individuals (beyond the minimum effective of 10 individuals) that was proportional to its area (see Supporting information -Figure S2 for an histogram of patch areas).For each city, these probabilities summed to 1. Therefore, cities with many habitat patches had many populations, cities with large habitat areas had many individuals overall, and large habitat patches on average contained many individuals.On average, there were 23 individuals per population (median: 19), and the largest population included 229 individuals in a 1000 ha patch.Consequently, population sizes varied among patches and among cities, thereby affecting drift intensities.Yet, drift intensity was constant across dispersal scenarios in every patch, making it possible to directly attribute genetic response variations across scenarios to variations in dispersal patterns.
Dispersal and reproduction parameters.Dispersal between populations depended on cost-distances computed for the three cost scenarios.10 % of the individuals of each population could disperse from their origin population to another at each generation, over a total of 250 generations.The dispersal probability from population i to population j decreased with cost-distance, such that: p d ij = e −αd ij .We used the same α value as specified in the previous section.
The population size was constant over time and the sex-ratio initially equal to 1.After birth, individuals potentially dispersed (see above) and could then only mate and reproduce with individuals from the same population.Each female had 3 offspring, with a sex-ratio averaging 1, and supernumerary individuals, either juveniles or adults, died to keep the local population constant.Initial genotypes were randomly generated by drawing 1 in 20 alleles, twice for each of the 20 simulated loci and for every individual.Mutations could happen at a 0.0005 rate.We carried out the simulations using PopGenReport R package (Adamack and Gruber, 2014).

Genetic data analyses
We wanted to assess the influence of the dispersal scenarios on the genetic responses simulated in each type of habitat.Thus, at the end of the simulations, we assessed both (i) intrapopulation genetic diversity (local) and (ii) inter-population genetic differentiation (pairwise), and their variations according to the scenarios and types of habitat (i.e., Forest, UGS, and their respective "interface" patches).Besides, to gain insights into gene flow influence on spatial genetic patterns in cities, we assessed the scaling properties of isolation by landscape resistance patterns, and whether populations formed genetic clusters coinciding with the spatial structure of habitat patches.Genetic analyses were performed using the graph4lg package (Savary et al., 2021b) in R (R Core Team, 2022).
Isolation by Landscape Resistance (IBLR).Isolation by distance (IBD) patterns have been known for a long time for providing insights into the relative influence of drift and gene flow on genetic differentiation (Hutchison and Templeton, 1999;Slatkin, 1993).We wanted to assess to which degree differences in spatial patterns of dispersal due to cost scenarios could affect the spatial genetic structure.We therefore analyzed isolation by landscape resistance patterns (IBLR) and their scaling properties.For that purpose, we iteratively computed Mantel correlation coefficients between D PS and cost-distances, while filtering population pairs according to increasing cost-distance thresholds until all population pairs became included.We identified the threshold at which this coefficient was maximal and called it Distance of Maximum Correlation (DMC), following Van Strien et al. (2015).To obtain comparable values across cities, we standardized the range of the DMC between 0 and 1 by dividing it by the maximum cost-distance between population pairs.A value of 1 is supposed to indicate cases where IBLR leads to a continuous and linear relationship between D PS and cost-distances at the scale of the whole study area (i.e., the equivalent of the case I IBD pattern sensu Hutchison and Templeton (1999), although we did not use the exact same framework as for classical Isolation By Distance analyses (Rousset, 1997)).On the contrary, values between 0 and 1 could indicate the presence of a plateau in the relationship (case IV IBD pattern sensu Hutchison and Templeton (1999)) (Van Strien et al., 2015).Besides, the correlation coefficient value between D PS and cost-distances at the DMC showed us to what extent genetic differentiation increased due to increases in cost-distance.
Module partitions.In each city and for each scenario (n = 325 × 3), we modeled population genetic structure using a genetic graph.Each node represented a population, and the links were weighted by D PS values.To identify genetic clustering patterns potentially emerging from dispersal limitations between sets of habitat patches, we identified modules in these graphs using the fast greedy modularity algorithm by Clauset et al. (2004).This algorithm identifies the partition of populations into modules maximizing a modularity index.This index takes genetic differentiation values into account and takes large values when populations from the same module are also genetically similar.
We wanted to determine (i) whether the spatial distribution of genetic modules reflected the dispersal constraints imposed by the different cost scenarios, and (ii) whether populations from different types of habitat tended to be assigned to different genetic clusters as well.We therefore compared the genetic modules with (i) modules computed in similar population graphs with links weighted by cost-distance values instead of D PS values, and (ii) the classification of habitat patches into forests or UGS.In the former case, we set the number of modules in each spatial cost-distance graph to be equal to the optimal number of modules identified in each corresponding genetic graph.Then, we compared these partitions using the Adjusted Rand Index (ARI, Hubert and Arabie (1985)), following the method described by Savary et al. (2022).This index takes its maximal value (1) when two nodes from the same module in one graph also belong to the same module in the other graph.It is equal to 0 when two partitions are comparable to

Genetic response modeling
We modeled the different genetic responses as a function of the dispersal cost scenarios (and habitat type or type of population pairs for the allelic richness and genetic differentiation, respectively) using mixed-effects models with random intercepts at the city level.This allowed us to take into account the lack of statistical independence between simulations made for three cost scenarios in the same city.In that way, we also accounted for the fact that all cities do not have the same habitat area, number of populations and individuals, which created overall differences in genetic structure, irrespective of cost scenarios.
Because all our genetic responses did not have the same range of values and distributions, we used mixed-effects models assuming various distributions and link functions (when generalized), as explained in the Results section and Supporting information 3. Models were fitted with a Restricted Maximum Likelihood approach using the lme4 (v1.1-30) (Bates et al., 2015) and glmmTMB (Brooks et al., 2017) R packages.The adequation between the distribution of the residuals and the models' assumptions were checked using a simulation-based approach implemented in the DHARMa R package (Hartig, 2022).We only interpreted the models whose residuals matched these assumptions.Because we modeled simulated values for which we controlled the number of replicates, we did not interpret the p-values (which took the lowest value reported by the R software program in most cases) (White et al., 2014).Mean values per fixed effect level and their confidence intervals (95 %) were obtained using emmeans (Russel, V. L., 2022).
When both the cost scenarios and the habitat type (4 levels: Forest, Forest Interface, UGS, UGS Interface) or type of pairs (3 levels: Forest-Forest, Forest-UGS, UGS-UGS) were considered as fixed effects, we included an interaction between these two variables in the models.Indeed, variations in dispersal were not supposed to affect the genetic responses similarly according to the type of habitat considered.

Assessing results' sensitivity to city size and habitat amount variations among cities
We wanted to check whether the results were consistent across cities regardless of prominent differences in terms of city size and habitat cover.To this end, we analyzed the results by considering separate groups of cities, created from the quartiles of the city sizes (study area radii) and of the total amount of habitat (UGS and forest total area).Additionally, we assessed whether allelic richness contrasts between habitat types were similar when the connectivity of UGS (EC UGS.UGS ) was larger than that of forests (EC Forest.Forest ).

Habitat connectivity
The 325 cities delineated in this study were covered at 20 % ± 1 % by artificial areas and had a radius ranging from 5 to 40 km, with a mean of 10.8 km and a median of 8.5 km.Their population averaged 300,000 inhabitants in 2018 (median: 140,000, maximum: 6,000,000).In these areas, the overall amount of reachable habitat (EC ) was mainly due to forest patches and much less to UGS, with EC Forest.Forest largely higher than EC UGS.UGS , and than EC Forest.UGS , although to a lesser extent (Figure 2A).This contrast was stronger when modeling connectivity according to dispersal  cost scenarios 2 and 3 (Supporting information 2 -Figures S3 and S4).However, in some urban areas where forest areas are very limited, the contribution of UGS to the amount of reachable habitat was the highest, as reflected by the proportional share of EC UGS.UGS (Figure 2B).

Genetic diversity
The simulated genetic diversity varied substantially both among cities and among cost scenarios and habitat types (Figure 3).We modeled the mean allelic richness across populations for each city and type of habitat occupied by the populations (Forest, UGS, Forest Interface, UGS Interface), as a function of habitat type and cost scenario using a linear mixed-effects model (LMM) with random intercept at the city level.The random effects explained 78 % of the overall variance (ICC: 0.783) due to large differences among cities.The residual distribution was satisfactory, as well as the model fit (conditional R 2 = 0.91, marginal R 2 = 0.56).Although population sizes, drift intensity, dispersal rates and dispersal distances were constant in a given city from one scenario to another, the overall allelic richness decreased substantially from scenario 1 to 3 in all types of populations when taking among-cities differences into account (Table 2), as expected from Figure 3.This resulted from stronger constraints on gene flow exerted by artificial areas, roads, and water bodies, which modified the respective intensities of drift and gene flow.The effect of the interaction between cost scenarios and habitat types was much lower than their main effects (χ 2 values from Wald test: main effect of cost scenario: χ 2 =16037.6,main effect of habitat type: χ 2 =5960.9,interaction: χ 2 =890.1).In the cost scenario 1, allelic richness was in all habitats relatively higher than in other scenarios, although it took slightly lower values in UGS (Table 2).For a given type of habitat (Forest, UGS), the highest allelic richness was observed in Interface habitats.It strongly decreased from scenario 1 to 3 but the decrease depended on the type of patches in which the populations were located.In Forest and Forest Interface patches, allelic richness was halved from scenario 1 to 3 (Forest: 10.68 to 5.93, Forest Interface: 11.68 to 6.60).It was almost divided by three in UGS Interface (10.95 to 4.19), and almost by four in UGS (9.37 to 2.39) (see CI in Table 2).
Table 2 -Results of the mixed-effects model of the simulated genetic diversity.Predicted values and confidence intervals of the mean allelic richness across populations at the city level as a function of dispersal cost scenario, habitat type and their interaction."Forest Interface" patches correspond to the forest patches most connected to UGS according to the F Forest↔UGS metric, whereas "UGS Interface" patches correspond to the UGS patches most connected to forests according to the F UGS↔Forest metric.Figure 3 -Distribution of the mean allelic richness of populations located in "Forest", "Forest Interface", "UGS" and "UGS Interface" patches in the 325 cities for the three dispersal cost scenarios."Forest Interface" patches correspond to the forest patches most connected to UGS according to the F Forest→UGS metric, whereas "UGS Interface" patches correspond to the UGS patches most connected to forests according to the F UGS→Forest metric.n = 325 values per box

Genetic differentiation
The mean genetic differentiation (D PS ) between populations was high overall (> 0.6) and increased from 0.62 to 0.72 and 0.85 in scenarios 1, 2 and 3, respectively (Figure 4).Genetic differentiation was lower between forest populations than between UGS populations and took intermediate values between forest and UGS populations.We first used a generalized linear mixed model (GLMM) assuming a beta distribution and using a logit link function to model the D PS values as a function of the cost scenario and type of population pair (i.e., Forest-Forest, Forest-UGS, UGS-UGS).We also used a LMM assuming a normal distribution.The GLMM provided a slightly better fit than the LMM.However, its residuals were slightly over-dispersed, whereas the LMM residuals had a satisfactory distribution.Both models provided similar results and we here present the LMM results.The random effects (city-level random intercepts) explained 76 % of the overall variance (ICC: 0.759) due to large differences among cities.The model fit was good (conditional R 2 = 0.90, marginal R 2 = 0.58).There were large differences among both cost scenarios and types of population pairs, after accounting for among-cities differences (Table 3).In that case also, the effect of the interaction between cost scenarios and types of population pairs was much lower than their main effects (main effect of cost scenario: χ 2 =11824.2,main effect of population pairs type: χ 2 =3714.1,interaction: χ 2 =850.0).Overall, the genetic differentiation among populations was lower in scenario 1 than in scenario 2, and the latter lower than in scenario 3, in accordance with Figure 4 (Table 3).The genetic differentiation among Forest populations was overall lower than between Forest and UGS populations, itself lower than among UGS populations.The increase in genetic differentiation from scenarios 1 to 2 and 3 was slightly lower for pairs of populations located in forests (see mean values and their CI in Table 3).

Distance of Maximum Correlation (DMC)
. IBLR relationships were very different from one dispersal cost scenario to another.In scenario 1, the DMC was overall equal to the maximum costdistance between populations (Figure 5A reflected by the Mantel correlation coefficient at the DMC (Figure 5B), was less steep in this scenario, indicating that the increase of genetic differentiation was somehow limited when costdistances increased.In contrast, in scenario 2, the mean DMC was equal to 0.47 and the corresponding correlation coefficients were high.This implies that for a subset of populations separated by cost-distances lower than a given threshold, there was a strong linear relationship between genetic distances and cost-distances reflecting significant gene flow between neighbor populations.Finally, in scenario 3, the low DMC values (Figure 5A) suggest that the IBLR relationships were weak and that genetic drift had a much stronger influence than gene flow on genetic differentiation.The mixed-effects models did not have a satisfactory fit for the DMC.In contrast, they confirmed that the differences in the correlation coefficients measured at the DMC across cost scenarios, visible in Figure 5B, were substantial even when taking among-cities differences into account.The values predicted by the LMM were equal to 0.55 (95% CI [0.54, 0.56]) in the scenario 1, 0.81 (95% CI [0.80, 0.83]) in the scenario 2, and 0.78 (95% CI [0.76, 0.79]) in the scenario 3 (see Supporting information 3).However, given that the distribution of the DMC values in the latter scenario reflects an absence of IBLR relationship, the corresponding correlation coefficients should not be interpreted for the scenario 3.
Population graph modularity.The modules identified in genetic graphs best coincided with the modules from similar graphs with links weighted by cost-distances in scenario 2 (Figure 6A).Indeed, mean ARI values were equal to 0.12, 0.27 and 0.00 in scenarios 1, 2 and 3, respectively (corresponding median values: 0.02, 0.25, 0.00).Similarly, in scenario 2, genetic modules best reflected the distinction between forest and UGS patches (Figure 6B).This means that, in this scenario, it was more likely for two populations belonging to the same genetic cluster to be  6 -Adjusted Rand Index (ARI) comparing the partitions into modules of the genetic graph with links weighted by genetic distances (D PS ) with (A) module partitions obtained from similar graphs with links weighted by cost-distances or with (B) the classifications of populations into forest or UGS populations according to the type of patch they occupy.An ARI value is computed for each city and each cost scenario, such that n = 325 values per box close when considering cost-distances, and to be located in the same habitat type.Thus, the genetic structure of populations and either their spatial structure when taking dispersal constraints into account, or their habitat type classification (UGS vs forests) matched in a stronger way in this intermediate cost scenario.This was not the case when gene flow was less constrained by unfavorable areas (i.e., scenario 1) or, on the contrary, highly constrained by these areas (i.e., scenario 3).We did not interpret the results of the mixed-effects models of ARI values, because the strong heteroscedasticity and atypical distribution of their values prevented us from obtaining satisfactory models.

Consistency of analysis results among cities having different sizes and habitat amount
The heterogeneity of city sizes, due to the varying radius of the disks containing 20 % of artificial areas, had negligible effects on the results.In the supplementary materials, we provide the same figures as included in this section plotted separately after splitting the urban areas in 4 quartiles based on their total area (Supporting information 2 -Figures S5, S6, S7, S8, S9, and S10).Similarly, the results were highly consistent when considering cities including varying total habitat areas (Supporting information 2 -Figures S11, S12, S13, and S14).Finally, in the few cities where the connectivity of UGS was higher than that of forests due to small forest areas, allelic richness contrasts were consistent with the ones observed for the whole set of cities in scenarios 2 and 3 (Supporting information 1 -Table S3).However, in these cities, under the scenario 1, UGS and UGS interface patches tended to maintain higher diversity levels than forest patches.
Paul Savary et al.

Discussion
We simulated the genetic structure of urban tolerant forest species occupying forests and urban green spaces in 325 European cities, while varying species abilities to cross the least favorable areas.We thereby found that in urban contexts, variations in dispersal movement behaviors alone can shape highly variable genetic diversity contrasts between habitat types and isolation by landscape resistance patterns.The substantial variations in simulated genetic responses between forests and urban green spaces could be due to their connectivity differences, reflecting their respective extent, spatial configuration and location within the urban matrix.These results were obtained without making any assumption regarding the respective quality of these habitats in our simulations.Urban ecologists should thus bear in mind the strong influence that dispersal between urban habitats exhibiting different spatial distributions can have on genetic patterns when assessing the relative influence of dispersal, adaptation, resource availability or biotic interactions on species responses to urban environments.Our results also provide insights into connectivity modeling and biodiversity conservation in these contexts.

Variations in dispersal constraints can shape highly contrasted genetic patterns
In our simulations, as dispersal across the least favorable land cover types became more costly, genetic diversity tended to decrease in both forest and UGS, and genetic differentiation tended to increase between populations.Although overall genetic responses, such as the mean allelic richness, differed among cities, variations among cost scenarios for a given city were relatively consistent.Contrasts were mainly due to differences of total habitat areas and number of patches, which determined the number of populations and individuals, and consequently the intensity of drift in each city.In urban contexts, differences in effective population sizes among species are known for being an important driver of the variability of their genetic responses (Schmidt et al., 2020), and are frequently invoked as the main driver of genetic diversity (Miles et al., 2019).However, these population size differences do not explain the strongly contrasted genetic responses among cost scenarios we obtained for a given city because the number of individuals in each population stayed constant whatever the scenario.Thus, only gene flow variations can explain these contrasts.
The resistance of artificial areas, roads and water bodies to species dispersal varied according to the cost scenario, but the total cost-distance that a fixed proportion of individuals could cover at each generation did not vary.In other words, the scenarios essentially modified the spatial pattern of dispersal movements, but not the effective distance they could cover.Consequently, for a given number of dispersing individuals, dispersal paths crossing unfavorable land cover types became less likely to be followed from scenario 1 to 3 because the cost of paths crossing other land cover types (e.g., grasslands, agricultural areas) remained the same whatever the scenario.Our results thus reflect the potential genetic responses of several species having the same individual density within patches, the same absolute dispersal abilities, and dispersal rates, but different dispersal behaviors in urban environments.It is noteworthy that simulated variations in the spatial pattern of dispersal can reproduce the inter-specific variability of patterns of genetic diversity and differentiation commonly observed in empirical landscape genetic studies (as reviewed by Fusco et al. (2021) and Miles et al. (2019)).For instance, the very subtle differences in genetic diversity simulated under the scenario 1 recall the empirical results of Khimoun et al. (2020).In contrast, the sharper differences simulated under the scenario 3 are akin to the significant genetic contrasts between urban and non-urban settings observed by Delaney et al. (2010) in several species.Although our results do not provide an explanation for these specific empirical findings, they show that differences in dispersal patterns could be sufficient to generate similar genetic patterns.
When drift intensity is constant but the spatial pattern of gene flow changes, the relative influence of drift and gene flow on genetic differentiation is logically modified in a different way according to the respective location of populations.Variations in the spatial range at which the relationship between genetic differentiation and cost-distances was the strongest revealed how contrasted genetic differentiation patterns were across scenarios.Indeed, when gene flow was not strongly restricted across artificial areas, roads and water bodies, genetic differentiation only increased progressively as cost-distance increased, producing a continuous pattern of genetic variations at the scale of the entire urban areas.This pattern is similar to the case I IBD pattern, according to the typology by Hutchison and Templeton (1999), and suggests that stepping-stone dispersal movements can prevent strong genetic differentiation in species with good abilities to disperse within cities.On the contrary, this pattern was only apparent between subsets of populations in scenario 2, and genetic differentiation increased more strongly with cost-distances.This corresponds to the case IV IBD pattern, according to Hutchison and Templeton (1999), and indicates that considering species having intermediate abilities to move across cities, the continuous and progressive increase of genetic differentiation resulting from dispersal limitation is only observed at small scales, within well connected subsets of populations.Finally, when dispersal movements in unfavorable areas were highly constrained, the relationship between genetic differentiation and cost-distances flattened out because drift became the main driver of the genetic response, and not gene flow anymore, which somehow recalls the case III IBD pattern described by Hutchison and Templeton (1999).In sum, changes in dispersal cost scenarios led to contrasted genetic differentiation patterns because they modified the relative frequency of dispersal events over each path connecting two populations.In other words, they rewired dispersal networks.The consequences of these changes in dispersal spatial patterns reinforce previous calls for a better consideration of population network topology in landscape genetics (Savary et al., 2021a;Van Strien, 2017).
In landscape genetics, unbridgeable barriers to dispersal are commonly inferred by identifying landscape features separating spatially structured genetic clusters of populations (Manel et al., 2007;Safner et al., 2011).Interestingly, our results show that genetic modules best reflected the spatial structure of genetic differentiation patterns in the second cost scenario, and not in the third one, which exerted the strongest constraints on dispersal movements.We could have expected the opposite and this suggests that there is an analytical limit in our ability to detect barriers to dispersal when gene flow is so constrained that drift is the main driver of genetic differentiation.
Several studies have evidenced that other processes than gene flow drive different population dynamics and individual fitness in urban populations as compared with non-urban ones.For example, urban bird populations can feed on more diverse food items, sometimes at the expense of their quality (Sinkovics et al., 2021), and can exhibit shifted and/or more variable phenotypes (Thompson et al., 2021), due to plasticity or genetic adaptations.If our results do not deny the Paul Savary et al. 19 Peer Community Journal, Vol. 4 (2024), article e40 https://doi.org/10.24072/pcjournal.407existence of these well-known processes, they nonetheless call for a better consideration of dispersal spatial patterns when inferring the respective influence of different drivers of population genetic structure in urban areas.

Influence of the spatial distribution of multiple urban habitat types
Forest populations tended to maintain a higher genetic diversity than UGS populations and to be less differentiated than pairs of UGS populations.Besides, genetic differentiation levels measured between these two types of habitats were intermediate, as compared with the high and low levels measured within UGS and forest patches, respectively.Similar genetic responses to urbanization have already been empirically observed in several species, from birds and reptiles (Delaney et al., 2010), to rodents (DeMarco et al., 2021;Gortat et al., 2015) and larger mammals (Stillfried et al., 2017;Wandeler et al., 2003).In our simulations, they mainly stem from the fact that the contribution of forest areas to the overall amount of reachable habitat was much larger than that of UGS in most cities.Besides, UGS patches are usually smaller and also harder to reach due to their location within the urban fabric.This explains why, except in scenario 1, even in cities where UGS contributed more to the amount of reachable habitat than forests, similar genetic contrasts were observed.These differences in terms of area and connectivity between UGS and surrounding natural areas provide a likely explanation to previous empirical observations in urban landscape genetics, as habitat amount and connectivity are often mentioned as key drivers of urban biodiversity (Beninde et al., 2015).
The stronger relative isolation of UGS was also apparent in the genetic clustering pattern.In the second cost scenario, we observed that forest and UGS populations tended to form separate genetic clusters.This sub-structuration of genetic differentiation patterns within the urban core areas had also been empirically evidenced in striped field mouse populations in Warsaw (Gortat et al., 2015) or in wild boar populations in Berlin (Stillfried et al., 2017).However, we also observed that the presence of habitat patches well connected to patches of another habitat type seemed to locally buffer these differences among habitat types by promoting gene flow at the interface between forest and UGS.Therefore, the connectivity analyzes and genetic simulations together suggest that peri-urban and less anthropized areas can be important sources of biodiversity in cities when they are connected to intra-urban habitat patches, in accordance with previous simulations (Snep et al., 2006) and empirical observations (Stillfried et al., 2017).The corollary of this is the potential sink role of UGS, as previously raised by Lepczyk et al. (2017) and Verrelli et al. (2022).This could have a negative influence on the long-term persistence and genetic adaptation of wild populations both within and outside cities, and remains to be investigated.

Implications for biodiversity management in urban areas
Making urban planning policies compatible with the conservation of biodiversity is crucial.These policies are commonly based on biodiversity surveys and on the conservation of so-called "green infrastructures", including several types of natural areas, UGS and the corridors connecting them.Our results stress several points that should deserve more attention in that context.First, genetic diversity and connectivity differences between forest and UGS were substantial regardless of the spatial extent of the city under consideration.In most cases, the largest habitat areas and biodiversity levels are to be found in natural areas surrounding city centers and not in UGS.This should encourage planners to consider large areas including the most biodiverse and favorable places to wildlife, often located at the periphery of cities.
Our connectivity analyzes and genetic simulations in 325 European cities also suggest that urban planners should identify habitat interface areas and consider them as "gateways" through which species can move from less to more anthropized habitats, as suggested by Gortat et al. (2017).We could also expect these areas to play a crucial role for maintaining species diversity within cities as long as gene flow and drift effects affect single species genetic structure in a comparable way as colonization and extinction processes affect species diversity (Vellend and Geber, 2005).
Finally, our results confirm that species which can hardly move across artificial areas and roads will not maintain high levels of genetic diversity within cities.This can explain why some species are very rarely observed in urban areas, but this could also mean that urban populations of some species are already engaged in a local extinction vortex.As such, considering the longterm effect of urbanization on genetic structure and its potential consequences for population persistence is key for biodiversity management (Sarrazin and Lecomte, 2016), and particularly in urban areas where these processes can be fast (Szulkin et al., 2020).

Limitations and perspectives
Our analyzes focused only on urban tolerant forest species occupying forest and UGS.These species only represent a small proportion of urban biodiversity.Reproducing these analyses by considering either specialist species using another type of habitat or more generalist species could help obtain a broader picture of biodiversity dynamics in urban settings.Besides, most forest patches were located at the periphery of cities, whereas UGS were more central, and this peculiar spatial distribution largely affected our results.Assessing how the spatial distribution of other types of urban and/or peri-urban habitats affects genetic patterns would also be needed.
For comparison purposes, we here assumed that UGS mapped by the Urban Atlas database were suitable for urban tolerant forest species.However, habitat patches may not all be suitable for these species and many other patches could fit their needs within cities (e.g., private wooded backyards in residential areas).This limitation should encourage the use of fine-grained remote sensing data to map urban habitats with more accuracy.Yet, acquiring such maps for a number of cities providing sufficient statistical power in a standardized way remains a challenge.
Although the genetic simulations we carried out can help prioritize habitat patches in urban planning context, they only reflect the potential genetic responses of a single species.These simulations could be advantageously completed by empirical surveys for assessing whether they closely reflect actual ecological processes.Our simulation framework could also be implemented in areas where genetic data are already available for testing whether dispersal limitations can explain the empirically observed genetic structure.Besides, adding to the simulation local habitat features, and how they locally influence several species based on their niche optimum, breadth, and competitive interactions (alike in the meta-community simulation framework of Thompson et al. (2020)) could provide insights into species diversity patterns in urban areas.This could help determine whether connectivity restoration measures are always positive for biodiversity conservation in urban areas; which may not always be guaranteed when invasive species also benefit from them.
Finally, although substantial differences among urban areas in terms of reachable amounts of habitat (i.e., connectivity) probably explain most of the variability in genetic responses across Paul Savary et al. 21 Peer Community Journal, Vol. 4 (2024), article e40 https://doi.org/10.24072/pcjournal.407cities, we did not investigate the structural causes of these differences.Indeed, the interplay of urban form with the spatial configuration of forest and UGS likely determines habitat connectivity patterns and species genetic structure.Further research is needed for understanding these relationships and providing broad guidance on urban planning at a time when increased urbanization and biodiversity conservation too often seem incompatible.

Figure 1
Figure 1 -(A) Land cover map of a city under consideration (city of Besançon, East of France), represented by a disc including 20% artificial areas.(B) The 325 European cities considered in our analyzes.(C) Example of least-cost paths (brown lines) connecting forest patches (dark green) and UGS patches (light green).(D) Example of genetic simulation output.The purple and orange dot sizes represent the simulated allelic richness in forest and UGS populations, respectively, according to the cost scenario 2.

10
Paul Savary et al.Peer Community Journal, Vol. 4 (2024), article e40 https://doi.org/10.24072/pcjournal.407random partitions.It takes its minimal value (-1) when partitions are totally distinct, i.e., when two nodes from the same module are in different modules in the other graph.

Figure 2 -
Figure 2 -Absolute and relative values of the different EC components computed in the 325 cities according to cost scenario 1. (A) EC Forest.Forest (green), EC Foret.UGS (purple) and EC UGS.UGS (orange) divided by the total study area, for each city.(B) Respective contributions of EC Forest.Forest (green), EC Forest.UGS (purple) and EC UGS.UGS (orange) to the connectivity of the habitat network.The total connectivity value of the network is the sum of the three EC components, which is slightly different from the global EC value because of square root properties.

Figure 5
Figure 5 -(A) Distance of Maximum Correlation (DMC), computed as the threshold distance used for selecting the subset of population pairs giving the maximum Mantel correlation coefficient between genetic distances (D PS ) and cost-distances, in the 325 cities and for the three dispersal cost scenarios.The DMC is divided by the maximum costdistances between populations in the corresponding city and cost scenario and therefore ranges from 0 to 1. (B) Mantel correlation coefficients measured at the DMC.n = 325 values per box

Table 1 -
Dispersal cost values associated with each land cover type according to the scenarios considered.Scenarios 1 to 3 represent an increasing aversion of forest species for dispersal movements across water, roads, and artificial areas.

Table 3 -
Results of the mixed-effects model of the simulated genetic differentiation.Predicted values and confidence intervals of the DPS among populations at the city level as a function of dispersal cost scenario, type of population pairs, and their interaction.