Beyond the landscape: Resistance modelling infers physical and behavioural gene flow barriers to a mobile carnivore across a metropolitan area

Urbanization affects key aspects of wildlife ecology. Dispersal in urban wildlife species may be impacted by geographical barriers but also by a species’ inherent behavioural variability. There are no functional connectivity analyses using continuous individual‐based sampling across an urban‐rural continuum that would allow a thorough assessment of the relative importance of physical and behavioural dispersal barriers. We used 16 microsatellite loci to genotype 374 red foxes (Vulpes vulpes) from the city of Berlin and surrounding rural regions in Brandenburg in order to study genetic structure and dispersal behaviour of a mobile carnivore across the urban‐rural landscape. We assessed functional connectivity by applying an individual‐based landscape genetic optimization procedure. Three commonly used genetic distance measures yielded different model selection results, with only the results of an eigenvector‐based multivariate analysis reasonably explaining genetic differentiation patterns. Genetic clustering methods and landscape resistance modelling supported the presence of an urban population with reduced dispersal across the city border. Artificial structures (railways, motorways) served as main dispersal corridors within the cityscape, yet urban foxes avoided densely built‐up areas. We show that despite their ubiquitous presence in urban areas, their mobility and behavioural plasticity, foxes were affected in their dispersal by anthropogenic presence. Distinguishing between man‐made structures and sites of human activity, rather than between natural and artificial structures, is thus essential for better understanding urban fox dispersal. This differentiation may also help to understand dispersal of other urban wildlife and to predict how behaviour can shape population genetic structure beyond physical barriers.

However, the urban environment imposes much more on wildlife than the need to navigate a highly altered landscape. Animals often perceive humans as predators and avoid areas of human activity (Samia, Nakagawa, Nomura, Rangel, & Blumstein, 2015). Individuals from the rural surroundings of an urban area might thus face a behavioural barrier to enter urbanized areas. Within the city, species with the physical capability of crossing the urban matrix may face behavioural barriers if they avoid man-made objects (with their artificial structures or scents) or human presence per se. Different scenarios are thus conceivable for population structure and drivers of gene flow across the urban-rural continuum and the perception of human-induced risks may differentiate urban and rural populations beyond physical barriers. The role of behavioural limitations to movement has been frequently overlooked. Examining the functional connectivity-the connectivity of the landscape from the species' perspective (Tischendorf & Fahrig, 2000)-of the urban landscape would thus help to assess the relative importance of physical and behavioural dispersal barriers and thereby make an important contribution to understanding the ecology and evolution of wildlife in novel environments.
Molecular genetic methods permit inferences about wildlife dispersal without the need to collect extensive field data on individual movements (Frantz, Do Linh San, Pope, & Burke, 2010;Guillot, Leblois, Coulon, & Frantz, 2009). Recently, numerous studies of gene flow in urban areas have been published, but many of those focus on smaller and less mobile species that are thought to face major barriers in urban areas (Beninde, Feldmeier, Veith, & Hochkirch, 2018;Combs, Puckett, Richardson, Mims, & Munshi-South, 2018;Munshi-South, 2012). Studies on larger and more vagile species, in contrast, analysed the population genetic structure of animals from peripheral suburban populations or from isolated sampling sites within urban and rural areas (Blanchong, Sorin, & Scribner, 2013;Santonastaso, Dubach, Hauver, Graser, & Gehrt, 2012;Stillfried, Fickel, et al., 2017;Wandeler, Funk, Largiadèr, Gloor, & Breitenmoser, 2003). There is currently no thorough analysis of the population and landscape genetic structure of a vagile species in an urban-rural continuum available, using continuous individual-based sampling. This would permit to identify drivers of urban gene flow, including those unrelated to the physical properties of the landscape.
Landscape genetic methods are particularly suited to assess functional connectivity. Specifically, hypotheses on how the environment influences gene flow in a target species can be evaluated by statistically relating the distribution of genetic similarities among individuals to landscape characteristics (Cushman, McKelvey, Hayden, & Schwartz, 2006;Schwartz et al., 2009). Several statistical problems have been recently solved, such as the nonindependence among ecological distances and the subjective assignment of resistance values to environmental features (Peterman, 2018;Prunier, Colyn, Legendre, Nimon, & Flamand, 2015;Sawyer, Epps, & Brashares, 2011). Landscape genetic approaches are still evolving (Balkenhol, Waits, & Dezzani, 2009;Manel & Holderegger, 2013;Richardson, Brady, Wang, & Spear, 2016) and some methodological aspects remain relatively underexplored. For example, while a simulation study by Shirk, Landguth, and Cushman (2017) has suggested that not all genetic distance measures perform equally well in model selection, different genetic distance measures have not been tested with the same empirical data set.
Aiming to gain a more fundamental understanding of the impact of urbanization on wildlife populations at a large spatial scale, we here focus on a mobile mesopredator, the red fox (Vulpes vulpes). Red foxes are ecologically flexible (Voigt & Macdonald, 1984) and occur in various habitat types. Their populations prosper even in highly urbanized habitats. In Berlin, our focal city, the first reports of foxes date from the 1950s (Saar, 1957) and by the 1990s the entire city was populated (Börner, Wittstatt, Schneider, 2009). Their ubiquitous distribution in highly artificial and fragmented areas as well as their movement ecology make foxes an ideal model for this study. On the one hand, foxes are very mobile. Urban animals have been reported to routinely cross streets and even rivers (Adkins & Stott, 1998) and gene flow may be unhampered by the urban landscape. On the other hand, anthropogenic infrastructure could represent significant gene flow barriers for mobile carnivores (Riley et al., 2006) and both telemetry and genetic studies point towards the existence of distinct urban and rural fox populations (Janko et al., 2011;Wandeler et al., 2003).
Here, we used continuous sampling of individuals both within Berlin as well as the adjoining rural countryside to evaluate three hypotheses. (a) The null hypothesis was that, due to their niche breadth and mobility, foxes disperse unhampered throughout the city and urban and adjoining rural populations are panmictic. This predicts that the urban fabric has no influence on gene flow, resulting in the absence of population and landscape genetic structure. (b) The fragmentation hypothesis posits that fox dispersal was (solely) affected by physical barriers such as rivers, built-up areas and highways. Under this hypothesis, multiple physical barriers limit gene flow, resulting in several scattered genetic populations. (c) The urban island hypothesis (Gloor, Bontadina, Hegglin, Deplazes & Breitenmoser, 2001) expects that dispersal may (also) be affected by behavioural barriers, which are most likely to occur at the border of the city where the rural landscape changes into the urban environment. Accordingly, individuals within the city are habituated to manmade structures and human presence, while individuals from the rural surroundings are not and thus face a behavioural barrier to enter the urban area. This predicts two genetic populations resulting from limited gene flow across the city border. We further expect that urban foxes disperse along artificial structures and through built-up areas when crossing the urban matrix.
In order to examine these predictions, we used assignment-based population genetic approaches to identify the location of abrupt genetic discontinuities and resistance-modelling-based landscape genetic approaches to assess the functional connectivity of the landscape. We tested three genetic distance measures to address the performance of different genetic distance measures in model selection and to generate robust results.

| Study area, sampling and laboratory procedures
The Berlin metropolitan area ( Figure 1a) has 3.5 million inhabitants and covers ~900 km 2 . It has been steadily changing during the last century and independent villages and satellite agglomerations were incorporated into the city. Thus, the urban landscape structure is quite heterogeneous, ranging from extremely urbanized areas of dense housing and high proportions of impervious surfaces to districts where forests and lakes represent up to 75% of land cover. The city area includes around 2,500 city parks, some areas of agricultural cultivation, 160 km 2 of forest and several lakes. The countryside around Berlin is characterised by sparse urban agglomerations, agriculture and forest. The landscape transition from the countryside to (sub-)urban areas does not fully correspond to the administrative boundaries between Berlin and Brandenburg as there are several forests, lakes and green areas that reach into the city (Figure 1b).
These green spaces and lakes are commonly used as recreational areas.
Detailed information on laboratory procedures is given in Appendix S1.

| Population genetic analysis
To assess the suitability of the microsatellites for population genetic analyses, we tested each locus for deviations from Hardy-Weinberg and linkage equilibrium using genepop v.4.7.0 (Rousset, 2008). We also used genepop to calculate F IS values (Weir & Cockerham, 1984).
To avoid deviations resulting from Wahlund effects and isolationby-distance (Frantz, Cellina, Krier, Schley, & Burke, 2009), when analysing the full data set, we subsampled the complete data set to generate 10 data sets consisting of 24 spatially clustered individuals (details in Appendix S2). We tested each set for significant deficiency or excess of heterozygotes and linkage disequilibrium (LD) among loci using the Markov chain method in genepop with 10,000 dememorization steps, 20 batches and 5,000 subsequent iterations.
We used the false discovery rate (FDR) to account for multiple tests (Verhoeven, Simonsen, & McIntyre, 2005).
We used two Bayesian-based clustering methods to estimate the number of genetic subpopulations (K), structure v. 2.3.4 (Pritchard, Stephens, & Donnelly, 2000) and geneland v. 3.3.2 (Guillot, Estoup, Mortier, & Cosson, 2005). Running conditions and specifications are described in Appendix S3. We used Microsatellite analyser MSA 4.05 (Dieringer & Schlötterer, 2003) to calculate observed and expected heterozygosities and the number of alleles across all loci. The level of genetic differentiation between genetic clusters inferred by structure and geneland was assessed using F ST (Weir & Cockerham, 1984) in spagedi 1.4 (Hardy & Vekemans, 2002); significance was tested with 10,000 permutations of individual genotypes between populations. We analysed the complete data set using the Estimated Effective Migration Surfaces (EEMS) method (Petkova, Novembre, & Stephens, 2016).
It uses georeferenced genetic data and can identify locations of abrupt genetic discontinuities (i.e., gene flow barriers) in data sets characterised by isolation-by-distance patterns (details in Appendix S4). We plotted the results for the run with the highest log-likelihood, using the reemsplots package in r (Petkova et al., 2016).

| Landscape resistance modelling
Functional connectivity was assessed using resistancega 3.1-3 (Peterman, 2018). It calculates pairwise resistance distances between individuals and uses a linear mixed effects model based on genetic algorithms to maximize the fit of resistance surfaces to the data. The process is based on stochastic search algorithms that solve optimisation problems by mimicking processes of natural selection (Scrucca, 2013). The optimisation process uses log-likelihood as the objective function. Mixed models were fitted using the maximum likelihood population effects (MPLE) parameterization (Clarke, Rothery, & Raybould, 2002) implemented in the R package lme4 (Bates, Mächler, Bolker, & Walker, 2014). A simulation study by Shirk, Landguth, and Cushman, (2018) has shown that this linear-mixed-effects-model-based method had a high accuracy in model selection.
The extent of the study areas was therefore obtained by plotting a minimum convex polygon around the sampling locations. Then, based on dispersal distances obtained by capture-mark-recapture methods (Harris & Trewhella, 1988;Trewhella & Harris, 1990) we added a 5 km buffer around this feature.
Landscape classification was based on the German topographic cartographic information system ATKIS (Figure 1b, Gruenreich, 1992). Seven landscape categories, hereafter called environmental predictors, were considered to potentially influence gene flow ( Figure 1c-d, Figure S1). In addition to ATKIS, we used data from the 2012 Copernicus Urban Atlas (https ://land.coper nicus.eu/) that uses high-resolution remote sensing data to provide detailed land cover information of larger European urban areas and their hinterland. It comprised 27 different land cover types (Table S1) We converted all layers into grids. Not considering cells without data, each cell in the initial grid had a value of zero or one, depending on whether it contained a feature under consideration. For linear predictors and water bodies we used a priority rule, meaning that every grid cell containing a linear predictor was coded as belonging to that predictor, independently of the proportion of the cell it covered. For shape predictors, we used a majority rule, with the cell being attributed to the single predictor with the largest area within the cell. Grid cell size was set to 250 × 250 m, giving rise to 233,798 grid cells without 'no data' cells, when ignoring three geographic outliers ( Figure 1a). Since we only considered one animal per grid cell, this resulted in 286 individuals being included in the analysis. Given the ecology of the species and the occasional lack of precision of the location information, we considered this to be an adequate compromise between computation time and spatial resolution. When focusing on the individuals in Berlin only, grid cell size was set to 100 × 100 m, given the higher accuracy of the sampling location and the smaller size of the study area. This resulted in 125,728 grid cells without "no data" cells, and 184 individuals being analysed.

| Genetic distance measures
Interindividual genetic distance measures are not equally accurate in model selection, especially when faced with weaker genetic structure (see Shirk et al., 2017). We therefore compared the perfor- We used genetix v. 4.05.2 (Belkhir, 2004) to perform an FCA on a multiple contingency table of the genetic data and used the first 10 FCA axes as a compromise between model accuracy and noise generation (Shirk et al., 2017). We calculated an Euclidean distance matrix for all individuals from their values on each FCA axis using the r package Ecodist (Goslee & Urban, 2007) and refer to this distance measure as "FCA". We used the r packages Alleles in Space (Miller, 2005) to calculate D Nei and Adegenet (Jombart, 2008) to calculate D PS .

| Optimisation of resistance surfaces: Single categorical environmental predictors
We first used the ss _ optim() command in resistancega to optimise the resistance of single categorical environmental predictors and test model selection performance of the genetic distance measures. In order to complete these analyses within a reasonable time, we limited initial tests to the seven ATKIS predictors (five for Berlin only).
We performed a (pseudo-)bootstrap procedure using the resist.boot() command, which subsamples individuals and resistance matrices without replacement at each iteration, refits the MLPE model to different resistance distance matrices and recalculates AIC c scores. We sampled 75% of the observations at each iteration. This was done in order to assess the relative support of each optimised resistance surface and the robustness of the model selection results given different sample combinations. For each genetic distance measure, we assessed model fit based on the differences between corrected Akaike information criterion (ΔAIC c ) values. When comparing genetic distance measures, the measure that gave rise to the highest marginal R 2 values (while generating biologically meaningful results) was considered the most adequate.

| Multiple resistance surfaces
After optimising individual categorical resistance surfaces, the relevant variables must be combined into a composite resistance surface. This is necessary to test whether models with several landscape features are better supported than models with single landscape features and, ultimately, to gain an understanding of the functional connectivity of the entire landscape (Khimoun et al., 2017;Ruiz-Lopez et al., 2016).

| Multiple resistance surfaces: Automatically combining categorical predictors
We used the samples from Berlin, the best genetic distance measure and the five ATKIS categories to compare two approaches that combine categorical predictors into a single surface. Firstly, we used resistancega's all _ comb() wrapper function which automatically combines and optimises all possible combinations of the five categorical predictors and runs the resist.boot() command to conduct a bootstrap analysis. However, the ms _optim() command gives different resistance values to a linear feature depending on which other feature it overlaps with, which may lead to erroneous conclusions (Section 3).

| Multiple resistance surfaces: Single-surface optimisation for combining categorical predictors
We therefore also tested a second approach for combining categorical predictors into a single surface: Rather than letting resistancega automatically combine different surfaces, we applied the singlesurface optimisation (ss _optim()) procedure to resistance grids containing multiple environmental predictors, i.e. each grid contained N categorical predictors and each cell in the grid had a value ranging from zero to N, depending on whether it was classified as one of the N predictors or whether it was classified as matrix, i.e., the remaining uniform study area not containing the features under investigation.
We will refer to these grids as "multicategorical" surfaces (to differentiate them from composite surfaces obtained using all _ comb()).
The principle underlying the multicategorical models is to add individual predictors based on model support (AIC c values) but to only retain a new predictor if its addition improved support (ΔAIC c > 2 after a resist.boot() bootstrap analysis). Individual predictors whose model support was ΔAIC c < 2 with distance were not considered.
The optimisation for each feature or combination of features was performed twice and only included the distance matrix from the optimisation run with the lowest AIC c value in bootstrap analysis. We will refer to this as the "stepwise optimisation" procedure.

| Multiple resistance surfaces: Dealing with overlapping linear features
In order to assess the effect of the overlap of linear predictors, we considered all possible priority combinations of predictors. We tested, for example, individual surfaces where linear predictor 1 took precedence over linear predictor 2 at points of overlap and vice versa. We also tested the support of a surface where all points of overlap between linear features were classified as a distinct feature.
In each case, the combination with the highest model support after bootstrapping was retained for further analysis.

| Multiple resistance surfaces: Effect of initial cell values of multicategorical surfaces
Preliminary exploratory analyses suggested that in the stepwise optimisation, the initial cell values of a predictor influenced the optimised resistance value for the predictors (and hence model support). We therefore coded individual predictors relative to their resistance/ permeability inferred in the initial individual analysis. For example, in order to obtain the highest model support when manually combining two different categorical predictors in a single grid, a predictor inferred to be permeable had to be given a grid value of zero, a predictor resisting gene flow a grid value of two and all other cells a value of one. In order to test more formally whether the optimised resistance values were sensitive to the starting values of the input surface, we took multicategorical surfaces with different combinations of predictors that were retained in the stepwise optimisation procedure and inverted the values of the input surface. We performed a total of four independent optimisation runs for each initial and inverted surface.

| Multiple resistance surfaces: Optimising multicategorical surfaces of the whole study area
We also used a stepwise procedure to generate a multicategorical surface for the whole Berlin/Brandenburg study area. To gain a more detailed assessment of the interface between the city and the surrounding countryside as a possible gene flow barrier, we created a concave hull of the administrative city border using the ConcaveHull plug-in for QGIS (QGIS Development Team, 2018).
We then drew 1, 2, 3, 4 and 5 km buffers around the concave hull and used the ss _ optim() command in resistancega to separately optimise the resistance of all five inner and outer borders ( Figure S4).
We then used the boundary model with the highest support together with the six remaining ATKIS predictors to identify the bestsupported multicategorical surface. Again, we performed the resist.
boot() bootstrap analysis for each optimisation run, to circumvent potential problems with imprecise locations of individuals sampled in Brandenburg.
The best-supported multicategorical resistance surfaces for the Berlin and Berlin/Brandenburg data sets were used to predict movement/gene flow patterns across both study areas using circuitscape v.4.0.5 (McRae, 2006). Animal movement paths were inferred between all pairs of sample location as well as between 200 random locations generated for both data sets using arcmap v.10.3 and located along the border of the study areas.

| RE SULTS
After correction for multiple tests, locus V502 deviated from Hardy-Weinberg equilibrium (HWE) in five out of 10 subsampled data sets (Table S2), its F IS values ranging between 0.40-0.53 in these five data sets. The locus was thus excluded from further analyses. No other locus showed systematic deviation from HWE. Some loci were in linkage disequilibrium in some subsampled data sets, but no pair deviated more than once (Table S3). We therefore performed further population genetic analyses with all loci except V502. identified a band of (slightly) reduced long-distance migration rates that covered most of the city, but also extended to the south of the study area (Figure 3a). In the east of Berlin, migration rates were significantly lower than the overall average rate (Figure 3b).
Independently of the clustering method, the more urban cluster had reduced genetic diversity compared to the rural cluster (Table   S4) and the clusters were significantly differentiated from each other. Differentiation between both structure clusters (F ST = 0.026; p < .0001) was higher than between the two geneland clusters (F ST = 0.011; p < .0001). The EEMS contour plot of effective diversity illustrated that in southwest Berlin effective diversity rates were significantly lower than the overall average rate (Figure 3c, d).

| Optimisation of resistance surfaces: Single categorical environmental predictors
When considering the five ATKIS environmental predictors, the results obtained after bootstrapping (Table 1) were qualitatively similar to initial model results (Table S5). The three genetic distance measures did not converge on the same results in the model selection process (Table 1). In the analyses using D Nei and D PS , motorways was always identified as the most significant factor facilitating gene flow, with railways ranked as second best model (also facilitating gene flow) and all other models (except one: water bodies in one run using the D PS genetic measure) having a difference in ΔAIC c < 2 with the distance model. In the FCA, the difference between the distance model and all five predictors was large (ΔAIC c > 6.2), with the ranking of the five models remaining identical between the two independent optimisation runs (  Figure S7).

| Multiple resistance surfaces: Single-surface optimisation for combining categorical predictors & dealing with overlapping linear features
The two best-supported models in the FCA-based analysis of individual features were water bodies and railways. When combining these two predictors in a single-surface optimisation, the highest model support was obtained (after bootstrapping) when giving water bodies precedence over railways in the resistance grid (when water bodies overlap with railways, the cell is classified as water body; Table S8). When adding the next best-supported motorways model to the single-surface analysis, the highest model support was obtained when water bodies took precedence over motorways and motorways took precedence over railways (water bodies > motorways > railways; Table S9). After bootstrapping, the three best "overlap" models had almost identical model support (ΔAIC c < 2; Table S9). Water bodies always strongly impeded gene flow, while railways and motorways conducted gene flow. When adding all vegetation (and hence built-up areas) to each of these three overlap models in a single-surface analysis, the model with water bodies > motorways > railways was again the best-supported model after bootstrapping (Table S10). In summary, when only considering the ATKIS data, the best permeability model for Berlin included all five tested features (Table 2; Figure S8). Water bodies strongly resisted gene flow (resistance: 1574), railways (resistance: 1) and motorways (resistance: 4) enhanced gene flow. Built-up areas (resistance: 291) were more permeable than all vegetation (resistance: 494). The following analyses were based on the water bodies > motorways > railways overlap model.

| Multiple resistance surfaces: Effect of initial cell values of multicategorical surfaces
Multicategorical models whose starting cell values had been inverted gave rise to different optimised resistance values for the predictors and had a lower model support than the noninverted original multicategorical surfaces (Table S11).

| Multiple resistance surfaces: Optimising multicategorical surfaces of Berlin
When repeating single-feature optimisations but splitting the all vegetation predictor into the two predictors forest and arable/ green and the built-up areas into the five Urban Atlas categories, water bodies, railways and motorways were still the individual features with the highest model support (Table 2), with arable/green generating a better model support than the all vegetation model (Table 2). Similarly, two Urban Atlas land cover categories (sealing levels [S.L.] of 30%-50% and >80%) had higher model support than the predictor including all built-up areas ( Table 2). The fourth-best (arable/green) and the fifth-best (S.L. 50%-80%) individual models had similar model support (ΔAIC c = 0.5; Table 2).
The resistance value inferred for each single predictor is given in Table S12. A better-supported model was obtained when adding S.L. 30%-50% to the water bodies, railways, motorways ("mrw") model than when adding the arable/green predictor or both arable/green and S.L. 30%-50% to the mrw model (Table 2). Adding further single-feature predictors to the single-feature optimisation procedure in order of decreasing support (and testing all possible combinations when ΔAIC c < 2 between two individual predictors) resulted in three multicategorical models having comparable support ( Table 2). The overall best model ( Figure 4) included railways (inferred resistance value: 1), motorways (resistance: 2), S.L.
30%-50% (resistance 8), S.L. 50%-80% (resistance: 103), S.L.>80% (resistance: 469), water bodies (resistance: 784) as well as the remaining matrix (resistance: 282). Despite differences in the resistance surface values between the models, the circuitscape current maps for the best supported model and the model with the fewest predictors were very similar, both suggesting that gene flow within the city of Berlin mostly occurred along linear landscape elements (railways and motorways, Figure 4).

| Multiple resistance surfaces: Optimising multicategorical surfaces of the whole study area
All city border models and obtained better support than the distance only model and inferred the city border to resist gene flow (Table  S13). The best-supported model (the city border converted into a concave hull) had a marginal R 2 of 0.42. The city border concave was also the most significant single predictor influencing gene flow when considering all other predictors (Table 3; Figure S9). Considering the bootstrapping results of the FCA-based genetic distance only, five of the six remaining single-feature models better explained gene flow than the distance only model (the exception being motorways; Table 3). Forest and arable/green were the only environmental features inferred to facilitate gene flow (Table S14). Again, the three genetic distance measures did not converge on the same results in the model selection process (Table 3), with the support of the city border model in particular changing with the genetic distance measure considered. Also, the marginal R 2 values obtained with D Nei and and D PS were considerably lower than those obtained with the FCAbased measure (Table 3).
When performing a stepwise procedure to create a multicategorical surface, two multicategorical models had almost identical support. The overall best-supported model (after bootstrapping) contained city border concave, built-up areas, railways and water bodies (Table 4), where water bodies took precedence over railways (Table S15). While railways (resistance: 1) and the remaining habitat matrix (resistance: 2) enhanced gene flow, city border concave (resistance: 498) and water bodies (resistance: 70) provided a greater resistance than built-up areas (resistance: 6). The second-best model (with almost identical support) had the same predictors (with similar resistance values) but did not include railways.
Considering both models, the circuitscape maps did not show a clearly-defined corridor network in the Brandenburg countryside ( Figure S9).

| D ISCUSS I ON
In the present work, we aimed to assess the importance of physical and behavioural dispersal barriers to drive population and landscape genetic structure of the red fox across the Berlin metropolitan area. We found support for the fragmentation hypothesis with major water bodies and densely built-up areas resisting gene flow. Contrary to our prediction, however, these barriers did not create several scattered populations across the city, possibly TA B L E 1 Boostrap results of the single-predictor resistancega analysis for the city of Berlin We also found support for the urban island hypothesis and inferred limited gene flow across the city border, indicating an effect of behavioural barriers. Urban foxes further made use of artificial structures when dispersing through the urban matrix. Our results may thus suggest a hierarchy of drivers of genetic structure with a general behavioural effect and impediment through physical barriers underneath. However, the specifics of our results also suggest that genetic structure was relatively weak and, therefore, dispersal rates still high.

| Population genetic structure and gene flow
The genetic clustering algorithms both inferred K = 2 as the mostly likely number of subpopulations, yet they differed in the spatial distribution of the clusters. For geneland, the cluster boundary closely coincided with the administrative city border, whereas for structure the urban cluster mostly excluded the north and northeast of the city. The location of the structure-inferred genetic discontinuity approximately corresponded to the course of the rivers Spree and Havel ( Figure S5). EEMS also identified reduced TA B L E 2 Results of the multicategorical functional connectivity analysis for the city of Berlin

Notes:
The best-supported multicategorical surfaces combining different environmental predictors were obtained using a stepwise procedure: Individual predictors were added based on model support (corrected Akaike information criterion values, AIC c ), but only retained if their addition improved support of the multicategorical model (ΔAIC c > 2). Presented here are the bootstrapping results based on two optimisation runs (Table  S12) that were performed for each (combination of) landscape features. avg. AIC c , average of the AIC c values obtained for each model in 1,000 bootstrap iterations; k, number of parameters; ΔAIC c , difference in the avg; AIC c , values between the best supported model (lowest avg.AIC c ) and each subsequent model; avg.weight, average of the AIC c weights obtained for each model in 1,000 bootstrap iterations; avg.mR 2 , average marginal R 2 of 1,000 bootstrap iterations. Predictors are sorted according to increasing avg.AIC c values. S.L., sealing level migration (broadly) around the city of Berlin, but especially in East Berlin. Despite discrepancies, all three population genetic methods inferred the presence of a cluster located within the confines of the city. Furthermore, the landscape resistance modelling identified (a concave hull of) the administrative city border as the main barrier to fox dispersal in the study area ( Figure S9). Finally, the F ST -based approach and the EEMS method confirmed reduced genetic diversity within (parts of) the city compared to the surrounding countryside. Our results therefore provided general support for a genetic differentiation between urban and rural areas, i.e., the urban island hypothesis.
While the three population genetic methods inferred the presence of an urban island, they differed in its proposed location and composition. Different solutions for the partitioning of a data set may result from differences in the assumptions and algorithms underlying the statistical methods (Guillot et al., 2009) and the way they deal with weak or hierarchical genetic structure (Frantz et al., 2006;Puechmaille, 2016;Rowe & Beebee, 2007) as well as with deviations from random mating that are not due to physical barriers (e.g., isolation-by-distance, presence of relatives, Rodríguez-Ramilo & Wang, 2012). As all three methods inferred a 'circular' cluster in the centre of the sampling distribution and the diversity within the city was reduced, it appears unlikely that the partitioning was an artefact of an isolation-by-distance pattern ).
Perhaps the most likely explanation for the observed outcome is that population genetic structure is weak because of high dispersal rates in our vagile study species. A simulation study suggested that structure was efficient at inferring the correct number of genetic clusters even at lower levels of genetic differentiation (i.e., F ST = 0.02-0.03), but this was not necessarily the case for its accuracy in assigning individuals to these clusters (Latch, Dharmarajan, Glaubitz, & Rhodes, 2006). While, by definition, geneland infers abrupt genetic discontinuities, the deviation from IBD inferred by EEMS also appeared to be relatively slight (Figure 3). We therefore conclude that our results provided evidence for genetic differentiation between urban and rural foxes, but that dispersal between urban and rural areas was ongoing.

| Performance of genetic distance measures
While resistancega offers high potential to gain a species-specific understanding of the functional connectivity of the landscape, careful consideration of some technical aspects seems necessary.
In the present study, the fit of a model testing single categorical environmental predictors and its rank relative to other predictors clearly differed between genetic distance measures. In the simula-  (Dytham, 2011), whereas Factorial Correspondence Analysis (FCA) was designed for multistate categorical variables (She, Autemm, Kotulas, Pasteur, & Bonhomme, 1987) and is thus more suitable for the analysis of allele states. Analogous to Shirk et al. (2017), our 10-axes FCA metric led to a better model fit (in terms of marginal R 2 ) than the other two measures and generated biologically meaningful results. Future research will show whether this is a general feature of FCA and how much this depends on the number of axes included. With a modest strength of the genetic signal, a few large eigenvectors may have insufficient diagnostic power to infer more subtle processes. The geographical distribution of the target species may also matter (Shirk et al., 2017). We considered 10 axes to be a good compromise between accuracy and noise and (almost)

| Pitfalls in landscape resistance modelling
Our results show that a subtle understanding of gene flow requires the simultaneous consideration of multiple landscape features.

| The urban island
Our results provided general support for a genetic differentiation between urban and rural areas, i.e. the urban island hypothesis.
The observed genetic structure was relatively weak, indicating that some individuals from the surrounding areas do disperse into Berlin. With abundant high-quality food and a lack of hunting pressure, the city is possibly a better-quality habitat for foxes, despite an increased mortality. Urban foxes could therefore be expected to stay within the city and individuals from the surrounding areas to disperse into the urban area. However, there was no support for a constant influx of foxes from the countryside and the (genetic) exchange between the urban agglomeration and the rural countryside was sufficiently reduced to maintain genetic structure.
In line with this, a radio-tracking study of foxes in Zurich showed limited movement across an the urban-rural boundary (Gloor, 2002). Colonising urban areas may thus require behavioural shifts such as an improved tolerance of the presence of humans (Gloor et al., 2001). Such behavioural changes have often been interpreted as resulting from phenotypic plasticity, allowing habituation to humans (Bateman & Fleming, 2012;Kauhala, Talvitie, & Vuorisalo, 2016;Vuorisalo, Talvitie, Kauhala, Bläuer, & Lahtinen, 2014). However, work on urban birds suggested that avoidance of humans may have a genetic basis and urban colonisation may result from selection for fearless individuals (Carrete et al., 2016;Carrete & Tella, 2009;Møller et al., 2015). The presence of a genetically distinct urban population may thus result from a founder effect followed by limited urban-rural exchange due to differences in avoidance behaviour (see also below).

| Gene flow within the cityscape
Gene flow in Berlin foxes was hampered by physical barriers. The landscape resistance models identified major water bodies as the most significant predictor resisting gene flow in the urban area.
Contrary to our predictions, foxes did not freely move through the urban landscape. The best-supported multicategorical model(s) inferred highly urbanised areas (sealing levels >80%) to represent an important impediment to gene flow. On the other hand, urban fox dispersal did not depend on corridors of natural vegetation as it was described for other species (Goldingay et al., 2013;Munshi-South, 2012) either. While suburban areas with low degrees of imperviousness were inferred to be more permeable for dispersers than densely built-up areas, our results suggest that railways and motorways served as the main dispersal corridors. This last result is in line with results from radio-tracking studies in Edinburgh where railway lines were the main conduit for long-distance dispersal of male foxes (Kolb 1984).
Railway lines and motorways are highly artificial structures.
On the circular railway around the city centre, trains pass continuously day and night. Similarly, the multilane motorways connecting the districts of Berlin are extremely busy with high-speed traffic.
While railway-tracks are usually embedded within vegetated verges, motorways are not, and generally, dispersal along such transport infrastructure carries a high mortality risk (200-250 road-killed foxes are found in Berlin each year: Börner, 2014). Yet, what both landscape elements have in common (besides their linearity), is the absence of human activity, in terms of pedestrians and cyclists. The green spaces of Berlin, in contrast, are usually crowded. Although the actual mortality risk in green spaces and sparse built-up areas is low, they show less conductance to gene flow than motorways and railways (Figure 3), despite the latter's inherent mortality risk.
Consequently, foxes may use artificial structures as corridors but avoid areas of human activity (see also Table 5). Adkins and Stott (1998) reported that city foxes stayed shy and preferably used sites when human activity was low. The authors concluded that foxes do not avoid human constructions-but humans themselves. Beyond physical barriers, human activity may thus represent a significant impediment to dispersal in urban foxes.

| The landscape of fear
Over centuries, foxes have been intensively hunted by humans -and still are. Although no hunting is conducted within the city, foxes should thus maintain a certain level of "background fear" (see Laundré, Hernández, & Ripple, 2010). The concept of a "landscape of fear" (Laundré, Hernández, & Altendorf, 2001;Laundré et al., 2010) is frequently applied to foraging behaviour and predator-prey relationships, but the authors promote its consideration for various life history traits. It describes how fear (or predator-induced stress) affects how animals use landscapes. It is not the actual predation risk but the anticipation of risks that limits movement in a landscape of fear (Laundré et al., 2010;Lima, 1998). In the context of our study, this could indicate that human activity drives urban foxes into costly trade-offs as they primarily disperse along structures with little human activity (hence low perceived risk) but high inherent mortality risks. This result conflicts with a model of fearless individuals entering and roaming through the city. Rather, behavioural plasticity may have allowed some foxes to enter the city and facilitate habituation to human presence to some extent, modifying but not obliterating their landscape of fear.
Movement constraints imposed by human activity could be even more relevant for rural foxes that are less accustomed to human presence (see also Stillfried, Gras, et al., 2017). Our results show that rural foxes, unlike their city relatives, do not use artificial structures as dispersal corridors and that dispersal was limited by the city border ( Figure 3). It may thus not be the rural foxes' physical capacity to move but the fear to do so that hinders rural foxes from entering the urban island and prevents admixture.
No matter how the genetic differentiation arose, the urban island could persist due to additional behavioural movement limitations. Human presence may thus be the key driver of red fox dispersal behaviour and impact both the separation into rural and urban clusters as well as the dispersal processes within the urban area.

ACK N OWLED G EM ENTS
We are indebted to the Stiftung Naturschutz Berlin as well as the National Natural History Museum of Luxembourg for providing funding. We would like to thank Bill Peterman for his help with the resistancega analysis and Gerald Kerth for his support. We further acknowledge the technical staff of LLBB for sampling and four anonymous referees for helping to improve this manuscript.

AUTH O R CO NTR I B UTI O N S
A.C.F., K.B., and S.E.K. contributed to the design of this research.