Detecting a hierarchical genetic population structure: the case study of the Fire Salamander (Salamandra salamandra) in Northern Italy

The multistep method here applied in studying the genetic structure of a low dispersal and philopatric species, such as the Fire Salamander Salamandra salamandra, was proved to be effective in identifying the hierarchical structure of populations living in broad-leaved forest ecosystems in Northern Italy. In this study, 477 salamander larvae, collected in 28 sampling populations (SPs) in the Prealpine and in the foothill areas of Northern Italy, were genotyped at 16 specie-specific microsatellites. SPs showed a significant overall genetic variation (Global FST = 0.032, P < 0.001). The genetic population structure was assessed by using STRUCTURE 2.3.4. We found two main genetic groups, one represented by SPs inhabiting the Prealpine belt, which maintain connections with those of the Eastern foothill lowland (PEF), and a second group with the SPs of the Western foothill lowland (WF). The two groups were significantly distinct with a Global FST of 0.010 (P < 0.001). While the first group showed a moderate structure, with only one divergent SP (Global FST = 0.006, P < 0.001), the second group proved more structured being divided in four clusters (Global FST = 0.017, P = 0.058). This genetic population structure should be due to the large conurbations and main roads that separate the WF group from the Prealpine belt and the Eastern foothill lowland. The adopted methods allowed the analysis of the genetic population structure of Fire Salamander from wide to local scale, identifying different degrees of genetic divergence of their populations derived from forest fragmentation induced by urban and infrastructure sprawl.


Introduction
Habitat destruction and degradation are considered the most important threats for wild populations, but the effects produced by the overall habitat loss are often complex to understand because many impacting factors do not act separately, but rather they play cumulatively or interactively (Giplin and Soul e 1986;Lindenmayer 1995;Young et al. 1996). These factors can shape the landscape, triggering a fragmentation process that affects the spatial distribution of animal populations by confining them to residual habitat fragments with a consequent reduction of population size. Small populations are usually more vulnerable to intrinsic demographic and genetic threatening factors, such as a higher variance of birth and death rates which leads to a higher probability of extinction and a less effective demographic response to environmental stochasticity. Moreover, small populations suffer from a higher genetic drift and inbreeding, leading to a loss of heterozygosity and genetic variability (Keller and Waller 2002;Frankham 2003Frankham , 2005. The mating of closely related individuals causes the inbreeding depression that has negative effects on demography (e.g., juvenile fitness and mortality rate among offspring; see Ralls et al. 1988;Lacy 1993;Lacy and Lindenmayer 1995), reducing population growth rates and thus the population size. In addition, the decrease of genetic diversity makes populations less adaptable to environmental variability (Frankham 2005). The smaller the population is, the more important are the effects of intrinsic threatening factors (Giplin and Soul e 1986).
Fragmentation processes generate metapopulations that are defined as a network of spatially discrete populations (i.e., subpopulations) linked by dispersal (Hanski and Simberloff 1997). The amount of dispersal between subpopulations represents the degree of their ecological connectivity. Several species naturally live in metapopulations, but their long-term persistence could be threatened by anthropogenic habitat fragmentation. Indeed, this process could lead to isolation, which is the disruption of dispersal movements and, consequently, the halting of the gene flow between subpopulations, which further emphasizes the negative effects produced by the habitat loss. Knowledge of the ecology of fragmented populations is essential in order to prevent their isolation or even restoring the ecological connectivity between them (Saunders et al. 1987;Burgman and Lindenmayer 1998). We studied the Fire Salamander Salamandra salamandra (AMPHIBIA, URODELA) genetic population structure because amphibians are usually model candidates for studies of fragmentation effects on connectivity (Moore et al. 2011). As they generally have low dispersal capabilities (Caldonazzi et al. 2007; Allentoft and O'Brien 2010) and are rather philopatric to breeding sites (Blaustein et al. 1994), amphibians are particularly susceptible to isolation. These characteristics often determine a high genetic differentiation among populations, even at restricted spatial scales (Allentoft and O'Brien 2010).
In the study area, located in Northern Italy, the species distribution still covers a wide area, although some subpopulations have declined or become locally extinct, where anthropic pressure is strong. Here, broad-leaved forest ecosystems which represent its habitat are mainly threatened by degradation, reduction, and fragmentation. Thus, the species ecology and distribution look particularly suitable to perform a general screening of forest ecosystems connectivity.
Molecular markers, such as polymorphic proteins or DNA sequences, are widely used for evaluating the effective genetic connectivity between populations, because they can distinguish breeding events, measure migration rates between generations and estimate gene flow (Avise 1994;Frankham et al. 2002;Frankham 2006). Moreover, molecular techniques require a lower sampling effort, as they usually rely on biological samples collected in a single period (Neville et al. 2006).
We chose microsatellites as molecular markers, because they are widely used in genetic population structure analyses and they have already been identified for the Fire Salamander in discrete numbers (Steinfartz et al. 2004;Hendrix et al. 2010). Microsatellites pertain to a noncoding DNA part of the genome, with no known function. This "neutral" region of DNA is thus particularly useful as it could change over time without bias induced by selection pressures. However, we stress that microsatellites have sometimes been suspected to be not completely neutral, meaning that at least some of the variation observed within and among populations may be attributed to selection (Kauer et al. 2003). Microsatellite markers generally have high mutation rates resulting in high standing allelic diversity (Selkoe and Toonen 2006). For this reason, when used for evaluating the ecological connection between populations, they should be identified in sequences with mutation rates which are low relative to the migration rates of individuals (Beebee and Rowe 2004).
The traditional Bayesian approach for the evaluation of genetic population structure from microsatellite data consists in the assignment of a genotyped representative sample of individuals to genetically homogenous groups. Usually, the sample is analyzed in one stage, and results can only highlight the main genetic structure of studied populations at wide scale. Nevertheless, some studies have stressed the usefulness of a multistep approach (V€ ah€ a et al. 2007;Harris et al. 2014), which entail a separate re-analysis of the groups identified in the previous steps, in order to identify hierarchical and more detailed genetic structure (Evanno et al. 2005;Pritchard et al. 2010). We adopted this method, as it can be very useful in large study areas with potentially fragmented habitats where local highly divergent populations can be hidden by the overall wide-scale genetic population structure.

Study area and sampling design
The study area is located in the Prealpine belt and in the foothill lowland of the Lombardy region (Northern Italy, Fig. 1). These areas were originally covered by extensive broad-leaved forests, which have been progressively removed and fragmented, especially during the last century. Forests have been replaced by a conspicuous urban sprawl particularly in the foothill lowland, while in the Prealps, they appear to be more continuous (Fig. 1).
Sampling design was realized in order to collect genetic data over most of the species range in the study area. We identified 28 sampling populations (SPs; Fig. 1), inhabiting suitable areas separated from each other by anthropogenic (i.e., continuous urban areas, main roads, etc.) or natural barriers (i.e., lakes, main rivers) (min distance: 2900 m). Within each SP, we selected between 1 and 5 sampling locations (making a total of 57), according to Fire Salamander ecology, breeding sites availability, and area accessibility (distance between sampling locations within SPs ranges from 1300 m to 9700 m). In each sampling location, we identified between 1 and 15 breeding sites (making a total of 174) corresponding mainly to slow-flowing streams and their meanders (ponds, in some cases) located in forest areas. The sampling location contained all close breeding sites within the same continuous forest patch in the same watershed (the distance between breeding sites within the sampling locations ranges from 12 m to 4100 m) (Fig. 2). This sampling design was real-ized in order to reduce the probability of having only siblings within the same sample unit (sampling location).
We collected 477 tissue samples by cutting the tip (about 3-4 mm) of 1-4 salamander larvae tail in each breeding site, all year round, from 2010 to 2013. We chose not to sample more than four individuals per site, in order to avoid possible bias deriving from full-sibling individuals when inferring population genetic structure by sampling larvae in breeding sites (Goldberg and Waits 2010). Tissue samples were stored in 95% ethanol in the field and subsequently kept in the laboratory at À20°C.
Salamander larvae were captured and handled with permission of the Lombardy regional administration:

DNA extraction and analyses of microsatellite markers
DNA was extracted using the Quick-g DNA TM MiniPrep kit (Zymo Research, USA) according to the manufacturer's instructions, eluted in 180 lL of elution buffer (10 mmol/L Tris-HCl, pH 8; 0.1 mmol/L EDTA), and stored at À20°C until subsequent handlings.
PCR products were analyzed in an Applied Biosystems 3130XL DNA sequencer (Life Technology, Inc.), and allele sizes were estimated using the software GENEMAPPER 4.0 (Life Technology, Inc.). Positive (known genotypes) and negative (no DNA) controls were used to check for laboratory contaminations, which never occurred. A 10% randomly selected subset of the other samples was PCRreplicated two times to check for allelic dropout and false alleles. Four of 20 microsatellite loci (Sal29, SST-F10, SST-G6, and SST-G9) were excluded from the analysis because we were not able to obtain PCR products that could be clearly interpreted.
We used Colony 2.0 (Jones and Wang 2010) to identify possible full-siblings, in order to avoid allele frequencies bias due to sampling larvae (Allendorf and Phelps 1981). We randomly selected one individual from each full-siblings dyad (Goldberg and Waits 2010) and removed the others from all further analyses. We used GenAlEx v. 6.501 Smouse 2006, 2012) to calculate the probability of identity (that is the probability of two independent samples having by chance the same identical genotype) in order to check whether the number of loci was suitable to identify univocally the individuals. The program LOSITAN was used to test for loci out of neutrality, which can influence the estimation of most population genetic parameters (Antao et al. 2008). The measures of the degree to which alleles at two loci were associated (linkage disequilibrium) was calculated using GENEPOP (Rousset 2008). We estimated the degree of genetic diversity among the 28 SPs, calculating global and pairwise F ST (Wright 1931;Weir and Cockerham 1984), and we tested their significance by a randomization test based on 9999 permutations. These estimated probabilities were corrected using the sequentially rejective Bonferroni method (Holm 1979), that is more powerful than regular Bonferroni correction (Allendorf et al. 2013). Finally, we calculated population genetic parameters, such as the number of genotyped individuals (N), number of different alleles (Na), number of effective alleles (Ne), allelic range (AR), observed (Ho) and expected heterozygosity (He), and fixation index (F IS ), for each locus on all individuals. Before performing genetic population structure analyses, we tested each locus for Hardy-Weinberg Equilibrium (HWE), checking its F IS , in order to detect an excess of homozygosity, which can bias population structure analyses. All these analyses were performed using GenAlEx v. 6.501. Nevertheless, as suggested by Allendorf et al. (2013), a deficit of heterozygosity may be caused by the presence of multiple demes (the Wahund effect, typically occurring in separate populations), sex linkage (occurring when some alleles are linked to one of the two sexes), and null alleles, but only the latter implies the removal of loci out of the HWE. We checked all loci for null alleles using Micro-Checker (Van Oosterhout et al. 2004) which can identify genotyping errors that can cause deviations from Hardy-Weinberg equilibrium and may bias population genetic analyses (Chapuis and Estoup 2007). Finally, we discharged only those loci not at equilibrium after sequentially rejective Bonferroni correction (see Allendorf et al. 2013), which also showed a probability to have null alleles.

Fire Salamander genetic population structure
Genetic population structure was performed analyzing the biparental multilocus genotypes by means of the Bayesian clustering procedure implemented in STRUCTURE 2.3.4 (Pritchard et al. 2000;Falush et al. 2003Falush et al. , 2007Hubisz et al. 2009), which was designed to identify genetically distinct clusters (populations of origin, K) of the sampled individuals. This analysis gave an assignment probability of each individual (Q i ) pertaining to the identified populations of origin. Populations were constructed by minimizing the departures from HWE and linkage equilibrium, which could result from recent admixtures, migration, or hybridization.
When genetic data have a weak signal, the analysis of genetic population structure by means of a Bayesian approach (as the one used by this software) can be improved with the knowledge of the sampling location of each individual. In fact, including the sampling location as a prior in the model, the clustering algorithm assumes that the probability of assignment of each individual to a population of origin varies among locations, without finding population structure where it does not exist (Hubisz et al. 2009). For this reason, we used the sampling location as LOCPRIOR in all performed analyses (see Tucker et al. 2012;Walsh et al. 2012;Kovach et al. 2013;Parsons et al. 2013).
As we did not know whether studied populations were discrete or had an admixed ancestry, we chose to run STRUCTURE using two ancestry models: "no-admixture" and "admixture". Actually, in some cases, it was unlikely that individuals from different SPs shared recent common ancestors (Falush et al. 2003), as the distances between SPs were higher than dispersal distance for several orders of magnitude. For this reason, we performed the analysis using the "no-admixture" model, which is also recognized to be more powerful in detecting subtle population structure (Pritchard et al. 2010). Nevertheless, in other cases, SPs could be relatively close and probably not completely isolated by barriers and might actually share an admixed ancestry. In these cases, the "admixture" model seemed to be more appropriate (Porras-Hurtado et al. 2013).
The effectiveness of the identification of populations of origin (K) may be further affected by the choice of the allele frequency model. Allele frequencies may show marked differences between distinct populations of origin and ancestral relationships between them is not expected (Rosenberg et al. 2005;Porras-Hurtado et al. 2013). In this case, the independent allele frequencies model seems to be more powerful in identifying populations highly divergent from each other (Pritchard et al. 2010) and reduces the likelihood of overestimating K (Hale et al. 2013). On the other hand, populations may show similar allele frequencies due to a shared recent ancestry (Pritchard et al. 2000;Rosenberg et al. 2005), and in this case, the correlated allele frequency model proves more powerful in detecting distinct populations of origin (Falush et al. 2003;Rosenberg et al. 2005;Porras-Hurtado et al. 2013).
The possible presence of different types of relationships between our SPs can hide a hierarchical structure that cannot be detected by a traditional one-step analysis. Indeed, in the presence of highly divergent populations, their separate analysis may improve the effectiveness of the correlated allele frequencies model in finding substructures (Evanno et al. 2005;Pritchard et al. 2010). For this reason, we adopted a multistep approach that consists of running STRUCTURE a first time in order to identify a main structure and then re-running the software separately for clusters identified in the previous step. Individuals were assigned to a cluster according to the best assignment probability (Evanno et al. 2005) of the sampling location (see LOCPRIOR above) in which they were collected. The process can be repeated until substructures of divergent populations are detected (see V€ ah€ a et al. 2007). At each step, we genetically characterized the identified clusters, testing all loci for HWE, and calculating F IS for each of them. Thus, we removed all loci out of equilibrium after sequential Bonferroni correction.
At each step, we ran STRUCTURE according to four parameter combinations using the two possible settings for "ancestry models" and "allele frequency models", because we did not have any information on the origins and the relationships of the sampled populations. All the runs for the "no-admixture" models were performed using a burnin period of 20,000 followed by 200,000 MCMC (Hastings 1970;Green 1995), while we chose to increment burnin period to 50,000 followed by 500,000 MCMC for the "admixture" models in order to reach an equilibrium in the estimated values of their key statistical parameters (see Porras-Hurtado et al. 2013). All the analyses were replicated 20 times for K ranging from 1 to 10. The optimal K values were selected by means of STRUCTURE HARVESTER (Earl and vonHoldt 2012) following the Evanno method (Evanno et al. 2005), based on the second order rate of change in the log probability of data between successive K values. In order to identify more detailed or subtle population structures, we also considered the K values with the highest average likelihood examining the plateau of the ln Pr(X|K) (Pritchard et al. 2000). The average population assignment probability Q p and the average individual assignment probability Q i were calculated using CLUMPP 1.1.2 with the Greedy algorithm (Jakobsson and Rosenberg 2007). The graphical displays (histograms) of genetic population structures were performed using the software DISTRUCT 1.1 (Rosenberg 2004), which plots individuals as colored bars according to their average Q i s.
Finally, we estimated the main genetic parameters of the groups identified in subsequent steps and calculated pairwise genetic distances among them in order to quantify the strength of genetic divergences and evaluate the ability of the multistep approach in discerning these genetic groups.

Genetic variability
Among the 477 sampled individuals, the Colony software identified 26 full-sibling dyads, corresponding to 23 pairs and three triplets of full-siblings. After the removal of full-siblings, the whole dataset consisted of 448 individuals. All loci resulted neutral to selection from LOSITAN analysis after Holm-Bonferroni correction. The probability of identity for increasing locus combinations (16 loci) resulted in values as low as 7.8 9 10 À12 . Values lower than 0.01 are believed to adequate for population studies (Waits et al. 2001). The panel of microsatellites thus supported reliable individual genotype identification. All loci were polymorphic, with the number of different alleles (Na) ranging from 5 to 24, and the effective number of alleles (Ne) varying between 1.03 and 5.23. All these parameters, with allelic range (AR), observed (Ho) and expected (He) heterozygosity, and fixation index (F IS ), for each locus, are shown in Table 1. We calculated F ST statistics on a subsample (436 individuals) to allow convergence in GeneAlex algorithm, deleting locus Sal23 (which had a high number of no data) from the initial pool of 16 loci and four SPs (SP1, SP10, SP16, and SP28) with fewer than five individuals. Significant overall genetic variation was observed (15 loci: global F ST = 0.032, P < 0.001). Pairwise F ST between the 24 SPs ranged from 0.010 to 0.101, with 66 significant comparisons after Holm-Bonferroni correction. Sampling populations SP5, SP6, SP7, and SP12 showed the highest number of significant pairwise tests (Table 2). In the linkage disequilibrium test, no association between alleles at two loci within SPs was found. Two loci of 16 resulted not at HWE after Bonferroni correction and were removed from the further genetic population structure analyses (SalE12 and SSTB11; Table 1).

Fire Salamander population structure
The first step of population structure analysis, performed on all samples (28 SPs, 448 individuals), produced consistent results among the four parameter combinations. According to the Evanno method, we identified two distinct clusters (K = 2; see Fig. 3A for the "no-admixture" model after CLUMPP analysis). CLUMPP Q p values, esti-mated by the four models, were always higher than 80% for both cluster 1 and 2, except for only one sampling population (SP12) for which the two "independent allele frequencies" models produced a Q p value of 58% for cluster 1. For the same SP, the two "correlated allele frequencies" models estimated a Q p value of 82% pertaining to cluster 1. For this reason, we considered the SP12 pertaining to cluster 1. The r LOCPRIOR parameter was always less than one, indicating that prior groups were informative in detecting population structure (Hubisz et al. 2009;Pritchard et al. 2010). The first cluster included all the Prealps area and the Eastern foothill lowland (and was named "PEF" group), with 24 SPs (323 individuals), while the second one corresponded to the Western foothill lowland (named "WF" group) with 4 SPs (from SP5 to SP8, 125 individuals; see Fig. 3A for the "no-admixture" model after CLUMPP analysis). Global F ST for the whole sample divided in the two groups was 0.010 (14 loci: P < 0.001). See Table 3 for all population genetic parameters and the HWE test.
Mean population parameters did not significantly differ between the two groups (Table 3). Nonetheless, the PEF population had a significantly larger number of sampled individuals than the WF and a higher number of private alleles (37 in the PEF and 9 in the WF). In particular, the SalE5 locus became monomorphic in the WF subsample (Table 3).    In the second step, we re-ran a population structure analysis separately for each of the two clusters previously identified (PEF and WF groups). We first tested for HWE using GeneAlex and deleted from the analyses all loci not at equilibrium after Holm-Bonferroni correction and with null alleles. In the PEF group, only Sal23 was out of HWE, while in the WF group, SalE5 resulted monomorphic and SST-A6-I not at HWE. We therefore ran population structure without locus Sal23 for the PEF and without loci SalE5 and SST-A6-I for the WF. Neither of the "independent allele frequencies" models, ran in this step for both the PEF and the WF, converged, even increasing the length of burnin and MCMCs, respectively to 100,000 and 1,000,000. We therefore considered only the "correlated allele frequencies" models hereafter.
The results of the PEF group were similar in terms of the number of clusters identified using the "admixture" and "no-admixture" models. Both models identified two clusters, according to both the Evanno and Pritchard methods. The "admixture model" showed a distinct cluster corresponding to SP12 with a Q p of 67%. All others SPs were assigned to a single cluster, with Q p s always higher than 69%. The "no-admixture" model supported the presence of 2 clusters (K = 2), and the results were consistent with the "admixture" model, identifying SP12 as a single cluster (Q p = 94%) and all other SPs pertaining to a second cluster (Q p always higher than 85%), except for SP10 which had a Q p of 67% (see Fig. 3B for the "no-admixture" model after CLUMPP analysis). Although the mean population genetic parameters did not significantly differ between the two clusters of the PEF (Table 4a), their Global F ST was significantly different from zero (13 loci, F ST = 0.006, P < 0.001). The two clusters showed a different number of private alleles, 33 in cluster 1 and 1 in cluster 2 (SP12). The latter also showed the locus SalE5 as monomorphic (Table 4a).
The analysis of the WF group produced different results using the "admixture" and "no-admixture" models, because the second one identified a more subtle structure. In the "admixture" model, the Evanno method allowed us to identify two clusters, with Q p s always higher than 73%. According to the Pritchard method, three clusters could be identified, although in this case, Q p s were never higher than 63% and it was showed the presence of two clusters corresponding to SP12 (blue) and all the others SPs (dark red); (C) the three histograms of the WF group (125 individuals, four SPs) showed population structures with 2, 3, and 4 clusters, respectively (see Results). For K = 4, cluster 1 in green, cluster 2 in red, cluster 3 in yellow, and cluster 4 in blue. These colors correspond to those in Fig. 4. Each bar represents an individual and is colored (using DISTRUCT) according to CLUMPP individual probability of assignment (Q i ) obtained from "No-admixture" models with "correlated allele frequencies" and LOCPRIOR information.
not possible to assign each SP to a defined population of origin. Conversely, the results of the "no-admixture" model suggested two clusters using the Evanno method and three to four clusters using the Pritchard one. In this model, Q p s were generally rather higher than the previous model allowing us to ascertain a stronger genetic structure. Assuming K = 2, Q p s were always higher than 86%, assigning SP5 and SP6 to cluster 1 and SP7 and SP8 to cluster 2. Considering K = 3 or K = 4, the general pattern of population structure did not change (Q p s > 84% [K = 3] and Q p s > 69% [K = 4] pertaining to the clusters 1 or 2). Nonetheless, the Q p s pertaining to cluster three or four were relatively low for all SPs. These two clusters were represented by two groups of genetically distinct individuals, each of which collected in a unique sampling location. The group assigned to cluster 3 was split from SP5 (Q i s > 93% with K = 3 and Q i s > 99% with K = 4), and the group assigned to cluster 4 was split from SP7 (Q i s > 66%) (see Fig. 3C for the "no-admixture" model after CLUMPP analysis). The mean population genetic parameters were not significantly different among the four cluster of the WF group, except for the number of different alleles (Na ; Table 4b). In cluster 2, locus SalE11 and SST-C3 resulted monomorphic (Table 4b). Global F ST for the WF group was 0.017 (12 loci: P = 0.058). Pairwise F ST among the four clusters of the WF group ranged between 0.007 and 0.083, but none of them were significantly different from zero after sequential Bonferroni correction (Table 5). Table 3. Population genetic parameters for the PEF (323 individuals) and Western foothill lowland (WF) (125 individuals) groups: number of genotyped individuals (N), number of different alleles (Na), number of effective alleles (Ne), allelic range (AR), observed (Ho) and expected heterozygosity (He), and fixation index (F IS ). Hardy-Weinberg Equilibrium test: degree of freedom (Df), chi-squared value (Chi-sq), P-value (p), significance after Holm-Bonferroni correction (HB sig).

Discussion
This research allowed us to detect a hierarchical population structure of the Fire Salamander in a wide area of the species range in Northern Italy. The 28 sampled populations were distributed in the Prealpine area, where forest cover is more continuous, and in the foothill lowland, where forests are more fragmented (see Fig. 1). Some SPs from the Western foothill lowland appeared genetically distant from the others, according to pairwise F ST (Table 2). This result was confirmed by genetic structure analyses which highlighted how these populations were strongly structured. We identified two groups, one represented by the Western Foothill populations (WF) and another one by all other populations laying in the Prealpine and Eastern Foothill areas (PEF). This was confirmed by the overall high probability of assignment of individuals to the two groups which resulted genetically distant as shown by significant pairwise F ST . Even if locus SalE5 became monomorphic in the WF group, the number of effective alleles (Ne) and the observed heterozygosity (Ho) were very similar between PEF and WF, suggesting that currently there was not a loss of genetic diversity on the whole. After the identification of these two genetically separated groups, the use of a multistep approach allowed us to identify substructured populations.
In the PEF groups, two more clusters were identified: A first one included all populations except SP12, which formed a second cluster whose separation was confirmed by both population structure analysis and pairwise F ST . SP12 cluster is located in a foothill forest separated from Prealpine areas by both wide conurbation and a main road. In addition, this area was also found to be isolated by the lack of a stream network connection that may support the dispersal process and seems to be important in connecting other Eastern foothill populations to the Prealps.
The WF group could be divided from two to four groups adopting the Evanno or Pritchard method. Considering K = 2, the group is divided into a Western (SP5 and SP6) and Eastern (SP7 and SP8) cluster, divided by large  conurbations and a south-north main road (Fig. 4). In addition, according to the K = 4 result, we observed another subdivision occurring in the two previous clusters (cluster 1: SP5 and SP6; cluster 3: SP7 and SP8), identifying two new clusters corresponding to two sampling locations, sampling location 5.5 (cluster 2), and sampling location 7.3 (cluster 4). Although these two new clusters are very close (with a distance of about three kilometers) to the other sampling locations of the same SP, they seem to represent separated demes. Cluster 2 appeared to be almost completely isolated by unsuitable areas for the Fire Salamander, largely represented by urban areas without stream network connections, while cluster 4 is isolated mainly by agricultural areas but with a stream connection, which could make the isolation less sharp (Fig. 4). Nevertheless, none of the pairwise F ST tests for the WF group resulted significant after Holm-Bonferroni correction. This result may be conditioned by the small sample size of clusters 2 and 4, each represented by a single sampling location, with respect to the SPs they were split from. On the whole, the identified population structure may be explained by the hypothesis of a relatively short-term fragmentation in the WF group: In this case, we cannot currently detect a significant loss of alleles, which could arise in the long term due to an inadequate gene flow. Conversely, it is likely that the PEF group (except SP12) will maintain a high number of alleles in the long term, because it is represented by a large and mainly continuous population (Fig. 3B).
Despite the fragmentation process currently only led to an evident population structure without significant differences of genetic parameters among identified clusters, in the medium and long term, it is expected it could even produce a loss of genetic diversity (e.g., Dixo et al. 2009). Indeed, the observed strong genetic population structure within the WF group is particularly worrying, because this will probably lead to a reduction in the capacity of adaptation of these populations, resulting in an increase in their probability of extinction, in relatively short time (Young et al. 1996;Reed and Frankham 2003;Arens et al. 2007). In addition, these populations will face other important threats such as further habitat loss and degradation caused by infrastructural, as well as commercial/ residential development projects in their area, which has one of the highest rate of urbanization in Italy. The same fate of the WF is expected for SP12, split from the PEF group, although its isolation seems to be less strong than the one occurring between the WF and the PEF groups and within the WF group.
From a methodological point of view, the multistep approach we adopted for the genetic population structure analysis at large spatial scale allowed us to identify different cases of fragmentation in the overall population, likely corresponding to a different degree of separation, or even isolation, of populations (see sampling location 5.5, i.e., cluster 2 of WF group). The effectiveness of this approach was proved by the use of the species here investigated, the Fire Salamander, an amphibian with rather strict ecological requirements (Di Cerbo and Razzetti 2004;Schmidt et al. 2005). Nonetheless, the genetic results suggest that at wide scale, most of the salamander populations are rather continuous or large enough to not suffer the typical problem of small population; they also highlight how locally the species may suffer habitat fragmentation. This is confirmed by the presence of strong structured populations, Figure 4. Population structure of the Western foothill lowland (WF) group. Numbered dots indicate sampling locations, grouped according to their sampling population (dot colors). Pie charts represent the sampling location probability of assignment to the four clusters identified by running STRUCTURE. Cluster 1: green, cluster 2: red, cluster 3: yellow, and cluster 4: blue. These colors correspond to those in Fig. 3C with K = 4. In the background main land uses: forests in green, urban areas in dark-gray, farmland, and other open habitats in gray, lakes, wetlands, and rivers in blue. Black lines represent main roads. with a generalized genetic depletion, evidenced by the presence of some monomorphic loci and many private alleles.
In conclusion, the multistep approach can be useful to detect substructured populations at local scale which cannot be observed analyzing wide-scale data in one-step only. Thus, the screening of genetic population structure from a regional to local scale can be considered effective in detecting local populations of conservation concern which resulted imperilled by threats acting locally and also for a preliminary assessment of the ecological connectivity in fragmented landscapes.