Inferring fine-scale spatial structure of the brown bear (Ursus arctos) population in the Carpathians prior to infrastructure development

Landscape genetics is increasingly being used in landscape planning for biodiversity conservation by assessing habitat connectivity and identifying landscape barriers, using intraspecific genetic data and quantification of landscape heterogeneity to statistically test the link between genetic variation and landscape variability. In this study we used genetic data to understand how landscape features and environmental factors influence demographic connectedness in Europe’s largest brown bear population and to assist in mitigating planned infrastructure development in Romania. Model-based clustering inferred one large and continuous bear population across the Carpathians suggesting that suitable bear habitat has not become sufficiently fragmented to restrict movement of individuals. However, at a finer scale, large rivers, often located alongside large roads with heavy traffic, were found to restrict gene flow significantly, while eastern facing slopes promoted genetic exchange. Since the proposed highway infrastructure development threatens to fragment regions of the Carpathians where brown bears occur, we develop a decision support tool based on models that assess the landscape configuration needed for brown bear conservation using wildlife corridor parameters. Critical brown bear corridors were identified through spatial mapping and connectivity models, which may be negatively influenced by infrastructure development and which therefore require mitigation. We recommend that current and proposed infrastructure developments incorporate these findings into their design and where possible avoid construction measures that may further fragment Romania’s brown bear population or include mitigation measures where alternative routes are not feasible.

Landscape genetics is increasingly being used in landscape planning for biodiversity conservation by assessing habitat connectivity and identifying landscape barriers, using intraspecific genetic data and quantification of landscape heterogeneity to statistically test the link between genetic variation and landscape variability. In this study we used genetic data to understand how landscape features and environmental factors influence demographic connectedness in Europe's largest brown bear population and to assist in mitigating planned infrastructure development in Romania. Model-based clustering inferred one large and continuous bear population across the Carpathians suggesting that suitable bear habitat has not become sufficiently fragmented to restrict movement of individuals. However, at a finer scale, large rivers, often located alongside large roads with heavy traffic, were found to restrict gene flow significantly, while eastern facing slopes promoted genetic exchange. Since the proposed highway infrastructure development threatens to fragment regions of the Carpathians where brown bears occur, we develop a decision support tool based on models that assess the landscape configuration needed for brown bear conservation using wildlife corridor parameters. Critical brown bear corridors were identified through spatial mapping and connectivity models, which may be negatively influenced by infrastructure development and which therefore require mitigation. We recommend that current and proposed infrastructure developments incorporate these findings into their design and where possible avoid construction measures that may further fragment Romania's brown bear population or include mitigation measures where alternative routes are not feasible.
Europe is characterised by a fragmented natural landscape, interspersed with a high density of human settlements and associated infrastructure. Habitat loss, fragmentation and infrastructure development are commonly regarded as being among the greatest threats to biodiversity. The construction of highways, has become an issue of increasing concern to wildlife populations 1 . Connectivity models are progressively being used as a tool to address the effects of fragmentation induced by barriers such as highways and other human infrastructure, and landscape resistance models provide a useful method for mapping potential mitigation corridors 2,3 . In a recent study it has been shown that eastern Europe has some of the highest road densities on the continent 4 . However, while 14.7% of the world's land mass is covered by Protected Areas (PA), globally only 7.5% is covered by connected PAs 5 .
One of Europe's largest brown bear populations (Ursus arctos Linaeus, 1758) is found in the Romanian Carpathians. While the population declined significantly after WWII, when only approximately 800 individuals remained, it recovered to almost 8,000 individuals by 1988 6 . Within the last two decades the population size has been estimated at 6,000 individuals 7 , and the population in the Carpathians has been shown to be one of the most Spatial autocorrelation. Results showed a higher and significant genetic correlation (class 5 km (P = 0.001, r = 0.064) and class 15 km (P = 0.003, r = 0.017)) among individuals than expected within specific classes of distance (Fig. 1a). Significant genetic correlation among individuals was detected for females within the first distance class (0-15 km) (class 5 km (P = 0.001, r = 0.10) and class 15 km (P = 0.001, r = 0.051)) ( Fig. 1b). For males we detected no significant spatial autocorrelation at any of the distance classes (Fig. 1c).
Landscape resistance modeling. In order for landscape models to be accepted, models had to pass the following tests: causal modelling with IBR (isolation by resistance) and the reduced model test. In addition, candidate models were also assessed against each other rather than simply against IBR to evaluate models based on relative support (RS). First, we will discuss the results for the best univariate models followed by the best multivariate models.
The best univariate models of effective resistance (Table 1) based on partial Mantel correlations were: rivers (r = 0.175; P = 0.003), roads (r = 0.091; P = 0.016), aspect (r = 0.076; P = 0.042) and slope (r = 0.051; P = 0.003). Large rivers explained about 18% of the variation in genetic distance, followed by roads, eastern aspects and by relatively small slope gradients. Elevation and land use were not significant and were therefore excluded from any further analyses.
Causal modelling with IBR for rivers was defined as such: partial Mantel test between genetic distance and rivers after removing the effect of IBR was significant (P = 0.002) and partial Mantel test between genetic distance  www.nature.com/scientificreports www.nature.com/scientificreports/ and IBR, removing the effect of the rivers was not significant (P = 0.997), suggesting that landscape resistance as a function of rivers was supported ( Table 2).
Model optimisation based on RS, showed that the best univariate model was rivers when compared to the IBR model, followed by aspect and roads ( Table 2). Only rivers, aspect and roads have met the IBR causal modelling criteria. R max values were greater for rivers and roads than in previous univariate models using casual modelling, while model parameters for aspect remained constant.
Causal modelling after removing the effect of IBR calculated for Ri + Ro + A (Table 3; Ri -rivers, Ro -roads, A -aspect) shown that: (1) partial Mantel test between the genetic distance and Ri + Ro + A after removing the effect of IBR (GD ~ (Ri + Ro + A)|IBR was significant (P = 0.037, Table 3; column 1); (2) partial Mantel test between genetic distance and IBR partialling out the landscape variable (GD ~ IBR|(Ri + Ro + A) was not significant (P = 0.671, Table 3; column 2); (3) the partial Mantel correlation was significant (P = 0.003, Table 3; column 3) when partialling out the effect of the reduced model (GD ~ (Ri + Ro + A)|(Ri + Ro); 4) while the opposite (GD~(Ri + Ro)|(Ri + Ro + A) was not significant (P = 0.959, Table 3; column 4). Model optimization based on RS, indicated that the best multivariate model comprised of aspect and rivers (  Table 3. The best multivariate models based on relative support (RS), causal modelling after removing the effect of the isolation-by-resistance (IBR) model (1, 2) and causal modelling criteria with the reduced model (3,4). We reported optimized parameter values, RS as compared to IBR, partial Mantel r and significance of support. Optimized values include equation parameters for x or SD (contrast; shape of the relationship) and R max (magnitude of the relationship). Causal modeling after removing the effect of the IBR consisted in: partial Mantel test between genetic distance and the landscape variable, partialling out the effect of IBR (GD ~ LV|IBR) and partial Mantel test between genetic distance and IBR, removing the effect of the landscape variable (GD ~ IBR|LV). Causal modeling criteria with the reduced model was: (3) partial Mantel test between genetic distance and landscape model after removing the effect of the reduced model GD ~ LM|) and partial Mantel test between genetic distance and the reduced model, partialling out the effects of the landscape model (G ~ |LM). Mantel r-value is reported in the first column for each test while the second column is reported the P-value. Abbreviation: Ri -River, Ro -Road, A -Aspect.
www.nature.com/scientificreports www.nature.com/scientificreports/ where (A) was the partial Mantel test between genetic distance and landscape variable 1 (LV 1 ), partialling out the effect of second landscape variable (LV 2 ) and (B) was the partial Mantel test between genetic distance and LV 2 partialling out the effect of LV 1 . In order to be supported, RS of the first model compared to the second model should be positive in every comparison. Including aspect in multivariate model 1 did not improve the RS value (0.097). The RS for model 2 compared to model 1 was slightly lower, but this model still passed the causal modeling criteria. This model included the same parameters as in the univariate analysis (Table 1). In addition, model Ri + Ro + A performed significantly better (r = 0.079, P = 0.037) than resistance distance alone (r = −0.018, The best multivariate model based on the partial Mantel correlation after partialling out the IBR model ( Table 3, column 1: P = 0.019; column 2: P = 0.964) included rivers and roads (model 1). Roads, rivers and aspect (model 2) were supported based on the partial Mantel correlation after removing the effect of IBR (Table 3).
The best-supported model using the AICc statistic consisted of a combination of aspect, rivers, roads and slope (C: R 2 = 0.021, w i = 1). When VIF values were calculated we observed collinearity between roads and rivers. We therefore calculated MRM excluding rivers (model C1; R 2 = 0.020, w i = 1; Table 4). When we excluded roads, model C 2 showed a lower AICc value (Table 4). Rivers, aspect and slope were identified as the landscape predictors influencing gene flow.

Current density maps and connectivity.
Using the best multivariate model, a map highlighting core areas with the highest current density (areas of connectedness) for brown bears along the Carpathians was generated. The current map of the most probable routes of gene flow indicated the presence of wildlife corridors which plays a major role in maintaining gene flow within this large brown bear population (Fig. 2). We therefore identified several wildlife corridors in areas were development has been planned ( Fig. 3a): Prahova Valley (Fig. 3b, 3 corridors), Olt Valley (Fig. 3c, 3 corridors), Apuseni (Fig. 3d, 3 large areas with small corridors) and Targu Mures -Iasi (Fig. 3e, 5 corridors). Some of the corridors are large, while others are very small. However, they are not the only ones, there are others that have not been the subject of this study. Moreover, we generated a framework of decision (Fig. 4) for new highways to be built based on our landscape genetic model outputs: (a) identifying if the infrastructure passing through an area where bears are known to occur; (b) identifying if the proposed  Table 4. Multiple regressions on distance matrices (MRM) indicating the relationship between pairwise genetic distances and the resistance distances for different landscape variables. In the model C1 and C2 we alternately removed roads and rivers in order to minimize colinearity among predictors. We reported β = intercept only for the best-supported model while P = P-value and VIF = Variance Inflation Factor were reported for each model. We present the results of matrix regressions (model R 2 ) and Akaike's Information Criterion (AICc, ∆AICc, w i ).
Models with the highest AIC support are in bold (∆AICc ≤ 2).
www.nature.com/scientificreports www.nature.com/scientificreports/ infrastructure bisect bear dispersal routes; (c) inferring if the proposed infrastructure pass through regions of high bear connectedness as a result of the landscape analysis; (d) analysing if the road can be re-routed without major disruption to the environment and at acceptable cost; (e) identifying mitigation measures (bear tunnels, green bridges, other solutions) which can be included into highway design; (f) given the landscape and results of the landscape genetics analysis, identifying the most suitable and cost-effective mitigation measure (bear tunnel at valley bottom, green bridge on valley slope with forest, etc.).

Discussion
Model-based clustering inferred one large population; the overall lack of genetic sub-structure suggests high levels of dispersal and gene flow. This implies that bear habitat has yet to become fragmented to restrict movement of individuals. This is also the case in nearby Croatia where a single population for brown bears was detected and highlighted the importance of habitat connectivity 36 . Spatial autocorrelation revealed a significant genetic correlation among females at relatively fine spatial scales (up to 14 km) supporting observations of female philopatric behaviour 14 . A previous study 37 showed that females can disperse up to 10 km while males disperse on average over distances of 30 km. These observations are supported by telemetry research carried out at the National Institute for Research and Development in Forestry Marin Dracea (INCDS Marin Dracea) in order to establish home ranges of brown bears in Romania (10 to 25 km 2 for females). The pattern detected of fine-scale spatial correlation in females could be explained by young females establishing their home ranges close to their mothers 14 , although we did not have information on the age of the approximately 50 females analysed here. As expected, males did not show significant spatial correlation, confirming ecological evidence for long-distance juvenile male dispersal 38 .
The optimal landscape model predicting genetic connectivity among individuals comprised a combination of rivers and eastern facing slopes. Large rivers, often located alongside villages and large roads featuring heavy traffic, were found to restrict gene flow, while eastern facing slopes promoted dispersal. These results are surprisingly similar as to those recently found for giant panda in the Qinling mountains in China, and could be linked to abundant vegetation and food resources 39 . In our study, land cover did not significantly correlate with gene flow. This is likely to be due to a lack of fragmentation of forest regions at a national level in Romania 40 . Similarly, other studies have reported a lack of correlation with land cover 34 , thus we infer that bear preferences for slope and specifically eastern facing slopes could be determined by forest type and food availability (as with the giant panda) 41 . A recent study carried out in the eastern Romanian Carpathians using telemetry showed that during hyperphagia bears seeking food were capable of moving along altitudinal gradients of high variability 42 , moreover, when analysing an extensive telemetry dataset of 70964 GPS locations collected from 32 brown bears (from which 9 were males) across Romania (INCDS Marin Dracea research projects, unpublished results) a high preference for www.nature.com/scientificreports www.nature.com/scientificreports/ the forest (especially mixed Norway spruce-silver-fir-beech and beech dominant forest) was revealed (Figs S5 and S6). In Romania, European beech and Norway spruce are found in large stands and high intrapatch/interpatch connectivity has been observed 40 suggesting that bears benefit from these favorable conditions in terms of food availability, shelter and movement. Slope aspect in beech forests seems to be an extremely important for herbaceous composition: a lower number of species but with more mature vegetation are present on north-eastern aspects from April to August, suggesting that north-eastern slopes could have higher biomass production due to more humid soil and competition for light 43 , facilitating thus bear movement in accordance with food availability. Moreover, Ciucci et al. 44 described the preference of the brown bears for consumption of herbaceous vegetation in spring and early summer, berries in latte summer and hard mast in the autumn (periodically available from beech (Fagus sylvatica), oak (Quercus) and common hazel (Corylus avellana)).
In our study, bears were inferred to use the gentlest slopes available and this likely reflects the energetic cost of dispersal 45 . However, topographic variables were not significant in a previous, broad-scale habitat model for brown bears across Poland 46 , while studies in other parts of Europe 47-49 implied a preference for rugged terrain. This observation might be explained by the geography of the Romanian Carpathian Mountain range, which is of low to medium altitude and no wider than 100 km, deeply fragmented by longitudinal and transverse valleys crossed by several major rivers; surrounded by large ground cavities. These features are expected to facilitate the movement of brown bear individuals by selecting the easiest paths due to energetic constraints.
When superimposed on the future infrastructure map, our landscape model suggested the presence of wildlife corridors in important areas for connectivity (Fig. 3a): (1) the Prahova Valley (Fig. 3b) also features a European-level road with extremely heavy traffic and a railway, crossing one of the most popular tourist destinations in the country (Sinaia, Busteni, Predeal, Azuga); (2) the gorge of the Olt River which includes many artificial lakes, and features extremely heavy traffic and a railway line (Olt Valley, Fig. 3c); while 3) Apuseni (Fig. 3d) and (4) Targu Mures -Iasi (Fig. 3e) are included in the future development plan. The influence of European roads/ highways is likely to be considerable, especially due to their collinearity with large rivers, and they are likely to have become more and more restrictive due to traffic intensification during the last two decades. However, highways with properly designed mitigation measures can ensure connectivity despite roads/railways/river systems.
Romanian bear habitat remains among the least fragmented in Europe (Joint EEA-FOEN report). If current traffic density plays an important role in restricting dispersal and gene flow in brown bears and road/railway mortalities confirms this (see Fig. S7), the proposed extension of the developing transport network could have a substantial impact on brown bears, which are not currently affected by dispersal constraints. Over the next few www.nature.com/scientificreports www.nature.com/scientificreports/ years, governmental road infrastructure development will include the construction of highways crossing core brown bear habitat and this is likely to cause further fragmentation if connectivity is not ensured. Since anthropogenic fragmentation acts rapidly and can affect the adaptive potential of populations 50 , investigating its effects post hoc can be obscured by time lags in response to disturbance 51 . However, landscape resistance models are a useful tool in mapping potential corridors and predicting connectivity prior to the construction process 52 , as we have attempted to apply here.
While economic development is seen as critical for any country, we advocate the use of landscape analyses to drive 'smart' development, by suggesting optimal mitigation measures. Restoration of wildlife habitats and ecosystem functionality is far more expensive than applying measures in advance to reduce fragmentation risks.
At a national level, landscape modelling can act as a decision-support tool for stakeholders to inform connectivity management and to encourage evidence-based conservation. Fine-scale maps of landscape resistance can be used to counteract the potential effects of habitat fragmentation and anthropogenic disturbance. Our approach can be seen as a means to inform authorities about ongoing ecological changes and could encourage a proactive approach to the design of infrastructure development through the proposed decision tool described in Fig. 4.
Areas where wildlife corridors might have the highest benefits for bear connectivity should be prioritised (Prahova Valley (Fig. 3b), Olt Valley (Fig. 3c), Apuseni (Fig. 3d) and Targu Mures-Iasi (Fig. 3d)) and barrier mitigation measures should be considered for improving connectivity e.g. considering these particular cases: (1) when European/national roads have very high volumes of traffic (while no highway is planned to be developed in the area) and wildlife corridors are identified and still used by brown bears, local measures should be applied: e.g. signs for reducing speed, warnings to drivers, improving driver visibility, installing an audible system to alarm wildlife when they are approaching a road; and (2) when a highway/railway is due to be constructed, mitigation measures such as overpasses and underpasses should be considered from the early planning stages (Fig. 4). Such an approach will reduce costs compared to when measures are adopted at later stages. While our results are robust with respect to recent changes in the landscape, fragmentation effects can be biased by time lags in species response, and corridor mitigation should be considered when implementing future developments. Brown bear conservation will depend on the management of corridors, monitoring of wildlife populations, evidence-based science outputs, collaborations between scientists, stakeholders and policy makers, and clear policies on landscape connectivity for wildlife species.
Furthermore, ensuring connectivity would have multiple benefits not only for biodiversity but also for other environmental objectives, while avoiding unnecessary or potentially bureaucratic mitigation planning. With these connectivity measures in place that will promote gene flow in the Carpathians, brown bear has a unique chance to co-exist with economic development. Statistical analysis. Microsatellite data was tested for the presence of null alleles using MICROCHECKER 58 .
Observed and expected heterozygosity was estimated for each locus using GENETIX 4.05 59 , while Hardy-Weinberg equilibrium was analysed for each locus using GENEPOP v.4.2 60 . Statistical significance was assessed by P-values using a sequential Bonferroni correction for multiple comparisons 61 .
We used a Bayesian clustering method in STRUCTURE v 2.3.4 62,63 to infer population structure using 1,000,000 simulations preceded by 100,000 burn-in replicates assuming an admixture model and correlated alleles frequencies. K values varied from one to ten and we performed ten repetitions for each K value. Both the posterior probability of the data for the given value of K (Ln Pr(X|K)) and its rate of change (ΔK) were used to identify the number of population clusters 64 . We generated genetic distances between individuals based on the proportion of shared alleles (Dps) in MSA 65 . We estimated pairwise Euclidean geographic distances between all sampling points allowing to test for the presence of a simple (isolation-by-distance) IBD 66 model, which assumes genetic differentiation is a by-product of geographic distance without taking any landscape features into account. IBD was calculated in GenAlEx v 6.5 67 using 10,000 random permutations. We investigated the genetic autocorrelation of multilocus genotypes at multiple spatial scales using GenAlEx v 6.5 67 . We then performed spatial autocorrelation, in addition, spatial autocorrelation was examined within males and females to test for sex-specific differences in autocorrelation. The 95% confidence interval around the point estimate of correlation (r) was obtained by bootstrapping with replacement. A value of r was considered statistically significant if its 95% confidence interval was above the null hypothesis of r = 0. We also performed multiple regressions on distance matrices MRM (R; "ecodist" package) 68 and this method was used to identify the contribution of each landscape variable to the overall variance in the dependent variable (genetic distance), since analysing data using only one method could result in false, method-dependent outcomes 69,70 . AICc weights (wi) were calculated by using the "MuMIn" R package 71 . Second-order Akaike Information Criterion values (AICc) were calculated among all competing models based on resistance distances, in order to select the best model (smallest AICc) 72 . Models that registered the lowest change in the AICc score (∆AICc = 0) and the highest Akaike weight were considered optimal 73 . To test for collinearity, the Variance Inflation Factor (VIF) was calculated for each predictor 74 in the model using the "VIF" R package. The presence of multi-collinearity, VIF values >10 75 revealed a linear association between two explanatory variables. Beta values (β) was calculated for the best model, in order to detect which of the independent variables had a greater effect on the dependent variable 76 using the R package "QuantPsyc 77 ". Raster maps for aspect, elevation, and slope were extracted from a Digital Elevation Model (DEM) at a resolution of 1 km square pixels. In addition, land cover (CLC 2012), rivers and roads maps with the same resolution were used. In order to assign resistance values to each cell the "Reclass" function 78 in ArcGIS v 10.1 (ESRI 2012) was used and rasters were converted to ASCII files to be used in CIRCUITSCAPE v.3.5 79 . The landscape parameters considered in this study are listed in Table S2 (Supplementary Materials).
Landscape resistance. An isolation-by-resistance (IBR) raster was generated in CIRCUITSCAPE v.3.5 by assigning a value of 1 to all pixels. Resistance distances between individuals were generated in CIRCUITSCAPE v.3.5, using the pairwise modelling mode with focal points and connected to eight neighbours 80 . We determined the best univariate models of effective landscape resistances based on partial Mantel correlation after removing the effect of the IBR model using the "vegan" R package 81 . In addition, a causal modelling approaches were used 82-84 that consisted of two steps: the first step suggests that if a resistance hypothesis model is supported independently of the null model then partial Mantel tests between the genetic distance and the environmental variable would be significant after removing the effect of IBR; partial Mantel tests between genetic distance and IBR distance would not be significant, partialling out the environmental variable. The second step allows for the comparison of causal modelling with a reduced model. If the true model is supported independently of the other candidate models then partial Mantel tests between genetic distance and the true model would be significant, removing the effect of the reduced model and partial Mantel tests between genetic distance and the reduced model would not be significant, partialling out the effect of the true model. Scaled transformations 85 for each landscape variable were ranked based on partial Mantel correlation coefficients after removing the effect of IBR using 9,999 permutations 19 . The function with the highest partial Mantel r-value (significant P-value) was included in the next step of the analyses. Further, the best univariate models based on relative support (RS) and causal modelling were evaluated after removing the effect of IBR, thus we tested if landscape resistance was supported independently by the null model 86 .
Models were ranked by model support values. Optimized values included parameter equations for x or SD (contrast and the shape of the relationship, respectively) and R max (magnitude of the relationship). The first step included a partial Mantel test between genetic distance and the landscape variable (LV), partialling out the effect of IBR (A = GD ~ LV|IBR); while the second step included a partial Mantel test between genetic distance and the IBR model, removing the effect of the landscape variable (B = GD ~ IBR|LV). In order for a landscape variable to (2019) 9:9494 | https://doi.org/10.1038/s41598-019-45999-y www.nature.com/scientificreports www.nature.com/scientificreports/ be included in the final multivariate model, it had to pass the causal modeling criteria with the reduced model. We report optimized parameter values, RS compared to IBR, partial Mantel r-values and significance of support.
Multivariate resistance surfaces were built by generating rasters equal to the sum of univariate model rasters for each landscape feature. First models were built using the best two landscape variables (rivers and roads, highest partial Mantel r) by maintaining parameters of the first variable constant, while the second variable parameters were varying. We identified the best-supported model for the second variable by analysing the partial Mantel correlation (after removing IBR). The optimum parameter of the first variable was obtained by keeping the parameter of the second variable constant. Additional landscape variables (aspect and slope) were added one at a time by keeping the first two constant until the best-supported model did not change. We also evaluated multivariate models by their RS rather than just removing the effect of IBR. Both procedures were repeated until all parameters were stable.
Cumulative current maps were generated in CIRCUITSCAPE v.3.5, best-supported maps of the landscape parameters that influenced gene flow have been generated by combining the cumulative currents maps in ArcGIS, displaying the road infrastructure layer and the current infrastructure development plan.

Data Availability
Distance (genetic and resistance) matrices; raw microsatellite data and GIS layers will be archived on DRYAD.