Defining Evolutionary Conservation Units in the Macedonian Crested Newt, Triturus macedonicus (Amphibia; Salamandridae), in a Biodiversity Hotspot

: In this study, we used genetic approaches to assess the conservation status of a protected amphibian species, the Macedonian crested newt, Triturus macedonicus , in Northern Pindos National Park (Epirus, Greece). Mitochondrial DNA sequences and multilocus genotypes of individuals from 38 breeding sites were used to infer their phylogenetic position and to detect and measure genetic variation patterns, population genetic structure, and levels of gene flow. The examined individuals fell within two major clades of the Macedonian crested newt phylogeny, being geographically sep ‐ arated by the Aoos River valley and Vikos Gorge. Both groups constitute separate gene pools, bear ‐ ing private haplotypes and alleles, and the groups were found to be highly differentiated in both their mitochondrial and microsatellite markers. Thus, they meet all of the criteria needed to be char ‐ acterized as evolutionary significant units (ESUs) that deserve a separate conservation status


Introduction
Conservation genetics encompasses the management of natural populations and aims to retain high levels of genetic diversity and minimize inbreeding, which reduces the risk of population and species extinction [1].The successful conservation of a species takes into consideration its evolutionary potential by resolving taxonomic uncertainties and delineating conservation units [2].Phylogeographic and population genetic analyses, apart from them significantly contributing to the resolution of taxonomic uncertainties, can also reveal areas of intraspecific diversity that are important for the conservation of a species [3].By recognizing and delineating conservation units, such as evolutionary significant units (ESUs) and management units (MUs) [4][5][6], the evolutionary context of a species can be effectively utilized in the implementation of targeted conservation actions [2].
Amphibians are often targeted for conservation and management programs because of worldwide population declines [7,8].Amphibians are particularly vulnerable to climate change and are likely to face serious survival problems; for example, 41% of species have been included in the IUCN Red List and are estimated to face some form of population decline in the future [9].Especially, pond-breeding amphibians face severe declines and extirpations due to extensive droughts and extreme temperatures that are caused by climate change, as water availability is linked to their life cycle and survival [10][11][12].In the Mediterranean region, due to the extended periods of high temperatures and droughts [13], amphibians suffer local population extinctions and geographical shifts [14].Functional connectivity, conducted through migration, is a process of major importance with respect to the maintenance of the structure, development, and longevity of amphibians' spatially structured population and metapopulation systems [15][16][17][18].Locally isolated populations are also at risk of further loss of genetic variation, genetic degradation due to random genetic drift, inbreeding, and reduced population fitness [19][20][21][22].Local isolation and loss of functional connectivity are crucial factors that affect the viability of pondbreeding amphibian populations and should be seriously considered in the design and implementation of effective conservation actions [23].
Among Mediterranean biodiversity hotspots, Northern Pindos National Park (Epirus prefecture, NW Greece) covers an area of 197,010 hectares of high landscape and geological complexity and hosts a large number of species of fauna and flora [24].Within the national park, numerous natural and artificial ponds constitute ideal breeding sites for several amphibian species [25].Among them, the Macedonian crested newt, Triturus macedonicus (Karaman, 1922), is an endemic salamandrid of the Balkan Peninsula [26][27][28][29] that shows a fragmented distribution; this is because it is usually associated with permanent ponds with rich aquatic vegetation [26] and prefers wetter conditions compared with its related species [28,30].Previous phylogenetic analyses recognized three distinct mitochondrial clades with relatively clear geographic distributions [31], all of which are represented in Greece: Clade 1, located in the northern distribution of the species, includes parts of northern Greece and Corfu Island; Clade 2 is located in a geographically restricted area in northern Greece and southern Albania; and Clade 3 is located in the eastern part of the species distribution in Greece [31].The distribution of species local breeding populations within the national park is highly fragmented and divided by an uninhabited area that is bounded by the Aoos River valley and Vikos Gorge [25].However, the phylogenetic position of these populations has not yet been assessed.
As a member of the Triturus cristatus complex, the Macedonian crested newt is included in Annex II of the Directive 92/43/EEC as a species of Community interest, whose conservation requires the establishment of special conservation areas (SCAs); meanwhile, a European action plan has been drawn up (Action Plan for the Conservation of the Crested Newt Species Complex in Europe) for its conservation and protection [32].Triturus macedonicus is currently designated as "Vulnerable" (VU), as its populations are decreasing due to the loss of its aquatic habitats through a combination of factors such as agricultural and livestock activities and the introduction of predatory fish [33].
Our study sought to assess the conservation status of the Macedonian crested newt inhabiting a protected area in the southern limit of its distribution.Our findings are expected to advance the current knowledge on the conservation requirements of the species in the area and may inform future conservation strategies.We employed two types of molecular markers, mtDNA haplotypes and multilocus microsatellite genotypes, to (a) complement the phylogeny of the species by resolving the phylogenetic position of the populations inhabiting Northern Pindos National Park; (b) infer the genetic conservation status of populations by assessing levels of genetic diversity, population structure, and connectivity through gene flow; and (c) address management priorities for the conservation of the species and its habitats in the area.We hypothesize that breeding individuals on the opposite sides of the major geographic barriers (Aoos River valley and Vikos Gorge) might constitute distinct and highly differentiated gene pools.

Laboratory Procedures
Tissue samples (toe clips and/or tail tips) from 309 individuals were collected from 38 water bodies in Northern Pindos National Park during April-July 2021.Additionally, seven samples collected in earlier samplings from sites outside of the national park were included in the subsequent genetic analyses (Table S1), resulting in a total of 316 samples.Total genomic DNA was extracted using the Macherey-Nagel NucleoSpin DNA Extraction kit following the manufacturer's protocol.
For the mtDNA analysis, a partial fragment of the mitochondrial subunit 4 of the NADH dehydrogenase gene complex gene (ND4) was amplified for 86 samples using the primer combination KARF4 (AGC-GCC-TGT-CGC-CGG-GTC-AAT-A) and DOBR2 (GTG-TTT-CAT-AAC-TCT-TCT-TGG-T) [34].Reactions were performed in a final volume of 20 μL containing 1X Taq buffer contained 1.5 mM of MgCl2, 0.4 μL of 10 mM of dNTPs, 2U Taq polymerase, 1 μL of each primer (10 μM), and 30 ng of DNA template, and the reactions occurred under the following conditions: an initial denaturation step at 94 °C for 3 min followed by 35 cycles of denaturation at 94 °C for 45 s, annealing at 58 °C for 45 s, extension at 72 °C for 60 s, and a final extension step at 72 °C for 5 min.The PCR products were purified according to the NucleoSpin PCR clean-up kit following the manufacturer's protocol.A Sanger sequencing of amplicons (~750 bp) was performed on an ABI 3730XL sequencer (Applied Biosystems) by CeMIA (Cellular and Molecular Immunological Applications, Larissa, Greece).
For the microsatellite analysis, 304 samples from 39 breeding sites were genotyped at 7 microsatellite loci, which were developed and tested specifically for the Macedonian crested newt in a previous study [42] (Table S2).Genetic loci were amplified by the PCR technique using forward primers labeled with fluorescent dyes (FAM, ROX, and HEX) in two multiplex reactions.Each reaction had a final volume of 12.5 μL, containing 15-20 ng of genomic DNA, 0.2 μM of each primer, 6.25 μL of KAPA2G Mix (Kapa Biosystems), and ddH2O.The PCR temperature profile consisted of an initial denaturation step at 95 °C for 3 min, 30 cycles of 15 s at 95 °C (denaturation-denaturation), 30 s at 60 °C (hybridizationannealing), 30 s at 72 °C (extension-extension), and a final extension step at 72 °C for 10 min.Genotyping of the seven microsatellite loci (Genescan 500 LIZ, Applied Biosystems) was performed using an ABI 3730XL sequencer (Applied Biosystems, CeMIA-Cellular and Molecular Immunological Applications, Larissa, Greece).Genotype scoring was performed using STRand ver.2.4.110[43].

Phylogeography
Mitochondrial haplotypes were identified in DnaSP ver.6.12 [44].We constructed an intraspecific phylogeny while following three different approaches: (a) maximum likelihood estimation (ML) [45], (b) Bayesian inference (BI) [46], and (c) neighbor joining (NJ) [47].The general time-reversible (GTR) model [48] with corrections for invariant sites (+I) was selected as the best-fit nucleotide substitution model based on the Bayesian information criterion (BIC) [49] for both the ML and BI approaches, which was incorporated using the Smart Model Selection software ver.1.8.4 (SMS) [50] implemented in the online version of the PhyML 3.0 program [51].The ML haplotype tree was constructed in the online version of PhyML 3.0, and the robustness of the topology was evaluated by 1000 bootstrap replications.The BI haplotype tree was constructed in MrBayes ver.3.2.2 [52], where four MCMC chains ran for 2 × 10 6 generations and were sampled every 100 generations; the first 25% of the generated trees were discarded as the burn-in.For tree visualization, we used FigTree ver.1.4.4 [53].The NJ haplotype tree was constructed using MEGA X based on the p-distance matrix using 10,000 bootstrap replications to assess branch statistical support.To visualize the relationships among the detected haplotypes, a median-joining (MJ) network [54] was constructed using PopArt ver.1.7 [55].

Historical Demography
To investigate the demographic history of the species in the national park, a mismatch distribution analysis [56] was performed in DnaSP software for the clades with a sufficient number of sequences.Three parameters were estimated for N0 and N1 (the initial population and the post-expansion population, respectively): τ, which refers to the time estimated from the beginning of the expansion in units of mutational time (τ = 2μt, where μ is the mutation rate of the genetic locus, and t is the number of generations since the beginning of the expansion), and θ0 and θ1, which refer to the size of the population before and after the expansion, respectively (θ0 = 2N0μ, and θ1 = 2N1μ) [57].Additionally, the goodness-of-fit and sudden expansion models were evaluated by 100,000 bootstrap replications.The sum of squared deviations (SSDs) between the observed and expected distributions was used as a statistical control, and for p > 0.05, the sudden expansion model can be accepted [56].Furthermore, Tajima's D [58] and Fu's Fs [59] tests for departure from neutrality were also estimated in Arlequin ver.3.5 [57].

Genetic Diversity
For each identified clade, the nucleotide (Pi) and haplotype (Hd) diversity were estimated in DnaSP.Pairwise genetic distances (p-distance) within and between clades, as well as between detected clades, were estimated using MEGA X using 1000 bootstraps.

Microsatellite Analysis
Individual genotypes were checked for selective amplification due to allele dropout and the presence of null alleles in Micro-Checker ver.2.2.3 (CI = 99%, 4000 iterations of MCMC-Monte Carlo Markov Chains); we only performed this check for sites with a sample size of greater than six individuals as Micro-Checker only considers the sites with a sample size of greater than six individuals in order to achieve an accurate dataset that is representative of the entire study area [60].

Population Structure and Gene Flow
The genetic population structure of all breeding populations in the study area was investigated using a Bayesian inference method in Structure ver.2.3.4.[63].Breeding sites with a small sample size (fewer than six individuals) were included to increase the representability of the area.Initially, to achieve unbiasedness in the model, parameter λ was calculated assuming that all sampled genotypes belong to one gene pool [64]; setting K = 1 with ten iterations; and using the "admixture model", "correlated allele frequencies", and LOCPRIOR (500,000 MCMC iterations, burn-in period of 250,000 iterations).The estimated λ was fixed for the following analysis using the same priors as before.To assess the optimal number of gene pools, ΔK was calculated using Structure Harvester [65] following the method of [66].An analysis of molecular variance (AMOVA) [67] was conducted using Arlequin to examine the distribution of genetic variation at three hierarchical levels: among the two geographic areas (eastern park-E.P.; western park-W.P.), among sites within areas, and within sites.The significance of the variance components was tested using 30,000 permutations.
To infer the gene flow rates among breeding sites, the rate of recent migration rates (i.e., during the last three generations) was calculated using BayesAss ver.1.3 [68] with 3 × 10 6 MCMC iterations, a burn-in period of 10 6 iterations, and a sampling frequency of 2000 iterations.The mixing parameters were adjusted to achieve an acceptable ratio ranging between 0.4 and 0.6 for the average number of changes per total iterations [69].The delta (Δ) values were set to 0.3, 0.2, and 0.65 for allele frequency (deltaP), migration rate (deltam), and inbreeding (deltaF), respectively.The 95% confidence intervals of the results obtained from the dataset were compared with the 95% confidence intervals for cases with insufficient data [70].

Genetic Diversity, Differentiation, and Isolation by Distance
Genetic diversity indices based on microsatellite data (average number of alleles per locus (A), number of private alleles (P), expected heterozygosity (He), and observed heterozygosity (Ho) were calculated in GenAlEx ver.6.5 [71].The Allelic richness (Ar) and Global FST or θ index (CI = 95%, 15,000 permutations) [72] were calculated using FSTAT.The pairwise FST was calculated in FSTAT, and its significance was tested by 10,000 permutations in Arlequin.The effective population size (Ne) of the gene pools was estimated using the program NeEstimator version 2.1 [73] using the molecular coancestry method [74].
For the entire park, and for each of the two geographic groups of breeding sites (W.P. and E.P.), the correlation between two matrices, a genetic distance matrix (FST/1 − FST), and a Euclidean geographic distance matrix (log km), generated by the dist() command in R software [75], was calculated using the Mantel test [76] (CI = 95%, 10,000 permutations) in GenAlEx.

Phylogeography and Historical Demography
Within the entire species distribution, the phylogenetic analysis identified 34 haplotypes.Four new haplotypes (Tmac31, Tmac32, Tmac33, and Tmac34) were added to the existing knowledge [27,31].Out of 658 nucleotide positions aligned, 583 were conserved and 75 were variable, all of which were coding for ND4.
Identical tree topologies were estimated by all of the different approaches (ML, BI, and NJ).However, some of the nodes are not highly supported for the ML and NJ methods (Figure 1b).Our phylogenetic analysis identified three main clades that differed by a substantial number of mutations (13-17 steps) while converging on an ancestral haplotype.
Clade 1 ("Northern Clade") consists of the individuals found mainly in the Northern distribution of the species, Corfu Island (Tmac16, 17,18,19) and Amvrakikos Gulf (Tmac13) (Figures 1 and A1).It includes 24 haplotypes that constitute an extensive complex mutational network, which originate from a central ancestral haplotype (Tmac04) (Figure A1).Clade 2 ("Zagori's Clade") includes six haplotypes and is further divided into two subclades.Subclade 2A includes individuals from Central Zagori and Ioannina Basin and shows a simpler typical star-like network, where some private haplotypes derive from a high-frequency ancestral haplotype (Tmac02) with one mutational step (Figures 1  and A1).Sub-Clade 2B includes individuals from the northwestern border of the national park (Tmac33) up to southern Albania (Tmac27) (Figures 1 and A1).Finally, Clade 3 ("Olympus' Clade") includes individuals from the eastern distribution of the species around Mount Olympus (Figures 1 and A1).The clade includes four haplotypes and shows a reticulation structure in which haplotypes are equidistant by one mutation step (Figure A1).The sudden expansion model appears to be supported for both Clades 1, 2, and Subclade 2A (Figure S1 and Table S3).The expansion time, represented by τ, was longer for Clade 1 than for the others.The θ0 value was approximately zero for all clades, indicating the possible existence of a bottleneck event prior to expansion.Regarding the neutrality tests, the D value was negative and statistically significant only for Clade 1 and Sub-Clade 2A (p = 0.047 and p = 0.010, respectively), indicating population growth after a bottleneck event [58].The Fs value was negative and statistically significant only for Sub-Clade 2A (p = 0.000), indicating an excess number of haplotypes due to recent growth [59] or genetic hitchhiking [77].Clade 2 showed a negative D value, whereas a positive Fs value reflected a deficiency in haplotypes due to a recent bottleneck event [59].Finally, the available data (nine sequences, two haplotypes) were insufficient for conducting a demographic analysis of Sub-Clade 2B.
The genetic distance (p-distance) between and within the clades and sub-clades of T. macedonicus and the related species T. carnifex and T. ivanbureschi is shown in Table 1.The greatest distance was observed between Subclade 2A and T. ivanbureschi.Within the T. macedonicus lineages, the largest distance was observed between Subclade 2A and Clade 3 (3.3%),whereas the smallest distance was observed between Subclades 2A and 2B (1%).The genetic distance between Clades 1 and 2 was 2.8%.S1.
Out of the new haplotypes found, Tmac31 was a private haplotype in site 5 and Tmac32 in site 30, whereas Tmac03 within the park was only detected once in site 6.The Tmac28 haplotype, which is a private haplotype of site 18 [27], was not recorded.

Microsatellite Analysis
All microsatellite loci were polymorphic in all studied breeding sites.The Micro-Checker results indicated the possible presence of null alleles at the TCM-83 locus in sites 27 and 30.For the remaining sites, this locus showed no evidence of null alleles, and the data were included in further analyses.No linkage disequilibrium (LD) was detected (p > 0.05) for all pairs of loci in all sites, while the inbreeding coefficient FIS ranged from −0.139 (site 24) to 0.198 (site 11) (Table S5); however, heterozygote deficiency and excess were not statistically significant (p < 0.00036).Specifically, the number of alleles per locus ranged from 2 (for TCM-83) to 18 (for TCM-517) (total of 47 alleles) (Table S2).

Population Structure and Gene Flow
The results obtained by the Evanno method showed the highest ΔK value (2672.3),which corresponded to K = 2; the second highest ΔK value (71.1) for K = 3 was considerably lower, suggesting that the optimal number of gene pools is two.The Bayesian clustering analysis of all individual genotypes indicates a clear population genetic structure, with extremely low admixture between the two distinct gene pools.Gene pool 1 (blue) consists of individuals from sites located east of the Aoos River valley (E.P.), while gene pool 2 (red) consists of individuals from sites located west of the valley (W.P.) (Figure 3).The AMOVA analysis revealed significant genetic differentiation between the two geographic areas (% variation = 17.82,FST = 0.23, p = 0.0000) (Table S6).The genetic differentiation among sites within areas (% variation= 5.16, FSC = 0.62, p = 0.0000) and within sites (% variation= 77.02, FCT = 0.18, p = 0.0000) was also significant (Table S6).
When simulating the effect of having insufficient data for migration rates estimation, we obtained values of 0.675-0.992(CI = 95%) for the non-migration rates (proportions of individuals originating from the "source" population of each generation) and 1.43 × 10 −19 -0.0933 (CI = 95%) for the immigration rates.For the non-migration rates, the confidence intervals that were recovered from the full dataset were significantly smaller than those obtained from the null hypothesis (insufficient data), suggesting that the dataset contained useful information to support the results [68].For the entire dataset, only seven sites exhibited significant migration rates (Table S7).Sites 27, 9, and 4 appear to be "source" populations, with site 27 being the most prominent, showing migration rates to four other sites (18,21,23,26).The highest migration rate (0.172) was found between "source" population 27 and site 26.Sites 4 and 9 showed significant migration rates to sites 5, 35, and 21, respectively (Table S7).

Genetic Diversity and Isolation by Distance
The mean number of alleles per locus (A) was estimated to be 4.9 in gene pool 1 and 4.4 in gene pool 2 (Table 2).In total, 25 private alleles were identified, 14 of which were private for gene pool 1 and 11 of which were private for gene pool 2 (Tables 2 and S8).The allelic richness (Ar) significantly differed, with a mean value of 3.0 for gene pool 1 and 1.9 for gene pool 2 (Table 2).The expected (He) and observed (Ho) heterozygosity was estimated to be 51.2% and 46.4% for gene pool 1 and 56.8% and 50.8% for gGene pool 2, respectively (Table 2).On the breeding site level, the mean number of alleles per locus ranged from 2.6 in site 25 to 4.0 in site 30 (Table S9).Two private alleles were identified: a 264 bp allele at the TCM-517 locus, private to site 4, and a 184 bp allele at the TCM-496 locus, private to site 38 (Table S9).The expected heterozygosity ranged from 40.4% to 56.3%, with a mean value of 48.35%.Accordingly, the observed heterozygosity ranged from 39% to 61.2%, with a mean value of 50.1%.The intra-population genetic diversity significantly differed between the breeding sites; the highest was observed in site 22 (Ar = 3.4,He = 0.537), and the lowest was observed in site 6 (Ar = 2.6,He = 0.404) (Table S9).Furthermore, the sites with the lowest genetic diversity were 6, 9, 25, 35, and 38 (Table S9).The overall (global) FST value in the national park was 0.131 (95% CI = 0.092-0.159),and the pairwise FST values ranged from 0 to 0.332 (76% of the comparisons were statistically significant).

Discussion
Our multimarker approach on the phylogeography and population genetic structure of the Macedonian crested newt within the biodiversity hot spot of Northern Pindos National Park adds new insights to the current phylogeny of the species as well as its conservation status in the area.

Identification and Characterization of ESUs
The phylogenetic reconstruction coincides with the one proposed by [31], which confirms the existence of the three major clades (Clade 1 or "Northern Clade", Clade 2 or "Zagori's Clade", and Clade 3 or "Olympus' Clade") but differs in the basal topology.Contrary to [31], our analysis seems to support the basal bifurcation between Clade 2 and the ancestor of the Clade 1 and Clade 3 branches (Figure 1b).However, the distinction between Clade 1 and Clade 3, although suggested by all three methods, is not as clear.The same basal node is also low supported in [31].Moreover, our analysis revealed the existence of two subclades within Clade 2 (Figure 1b).Our dense sampling in the study area and the four new haplotypes that were identified, three of which were recovered in the Zagori's Clade (Tmac32, 33, and 34; Figure 1b), might explain the finer resolution of species phylogeny.The degree of genetic differentiation between the three clades is considerably high (Table 1).Similar levels of genetic differentiation have been observed among recognized subspecies of the alpine newt (Ichthyosaura alpestris) [78], raising the question of whether the three Macedonian crested newt lineages might deserve a subspecific designation despite the cytonuclear discordance that has been observed for the species [31].In any case, sequencing of larger genome fragments and the analysis of additional genomic markers is expected to reconstruct the phylogenetic history of lineages with greater accuracy and fidelity (e.g., [79]).
The combination of high haplotype and low nucleotide diversity that was observed for both the entire species and the individual mitochondrial clades (Table S4) might reflect rapid demographic growth from an initial small source population [80].However, the recent expansion model is not supported for all the clades.Specifically, Tajima's D and Fu's Fs tests did not support the recent expansion model for Clade 2, whereas Tajima's D test did support the recent expansion of Clade 1. Clade 1 seems to have been expanded from an already identified refugium in the western coasts of Albania [27,28,31].In contrast, Clade 2 shows long demographic stability, and the populations probably survived in an independent refugium during the Last Glacial Maximum.Central Zagori was a suitable area that was absent of glaciers compared with the rest of the Northern Pindos area [81,82].However, the recent expansion model is supported for Sub-Clade 2A.This suggests that Sub-clade 2A individuals might have recently expanded across the Zagori Province and Ioannina Basin (Figures 1 and 2) after a bottleneck, which was probably caused by human intervention in the water table of the area [83].Such an expansion could have been facilitated by the existence of a dense network of permanent artificial ponds in the area, as has been the case in the close relative Triturus cristatus [84,85].
All Triturus macedonicus populations from the Northern Pindos National Park belong to two previously described phylogenetic clades.Sites from the eastern and northern parts of the national park constitute marginal populations at the southern periphery of the geographically extended northern lineage (Clade 1), while sites from the western parts of the national park belong to the geographically more localized Zagori lineage (Clade 2) (Figure 2).
These two population groups are geographically separated by physical features such as the Aoos River valley and Vikos Gorge (Figures 2 and 3), which appear to have been insurmountable barriers for the gene flow.Such a long-lasting physical separation could explain the observed high differentiation in both mitochondrial and nuclear markers (2.8% p-distance and 0.184 FST) and the formation of two independent gene pools with complete or extremely low admixture between them (Figures 2 and 3).Moreover, the detection of several private alleles in each gene pool indicates separate evolutionary histories in isolation [86].
Therefore, as these two population groups were found to be reciprocally monophyletic at the mtDNA level, differ considerably in their allele frequencies, and show unique genetic characteristics and diagnostic genetic diversity, they meet the criteria to be characterized as separate evolutionary significant units (ESUs), with particular importance for targeted conservation [4][5][6].

Genetic Variation and Population Connectivity
The overall genetic diversity of the Macedonian crested newt in the study area was found to be lower than that observed in other studies of closely related species at various spatial scales [87][88][89][90].
Although the studied breeding sites comprise two homogeneous gene pools, significant migration within each of them was detected on a few occasions.Only seven cases were observed: two in the eastern group and five in the western one.This can be indicative of reduced functional connectivity between breeding sites or might reflect a poor performance of the Bayesian approach that was used due to uneven sample sizes; thus, this observation should be treated with caution.However, 50% (18/36) and 59.2% (116/196) of the pairwise FSTs in the eastern and western genetic group, respectively, were found to be statistically significant and relatively high (0.05-0.23 and 0.02-0.19,respectively), having been higher than those found in similar studies [87,[89][90][91].Such levels of differentiation, independent of the geographic distance between sites (absence of IBD), implies local isolation and differentiation due to drift [1,20].The quite low Ne estimates in both geographic areas, being particularly low in the western one in spite of the larger number of sampled individuals, can be considered as indicative of a strong drift [92].
The significant differentiation among breeding sites of the eastern population group, the most prominent being between the sites residing in the northern and southern parts of the E.P., could be attributed to their geographic isolation, as the dispersal ability of Triturus newts is much lower than the between breeding sites geographic distance [85,93].
Although individuals from the western group occupy a smaller geographic area and breed in a dense network of artificial ponds, significant genetic differentiation and low dispersal between sites have been detected, which was probably due to the denser network of paved and dirt roads that might act as barriers to gene flow [94], promoting local differentiation even on a small scale (e.g., sites 15 and 34; Figures 3 and S2).The negative effects of road networks on amphibian dispersal have been addressed, such as the impact of road traffic on annual migration routes [95] and the high heat storage on paved surfaces, which poses a high risk of dehydration for amphibians [96].The negative impact of the road network on the population structure of Lissotriton graecus in the study area has been previously addressed [89].However, a more detailed landscape genetic analysis in the area could locate these barriers to gene flow more concretely.

Conservation Status and Management Priorities
Our results on the levels of differentiation (Figure S2), AMOVA (Table S6), and effective population size estimates (Table 2) indicate the local isolation of breeding sites and genetic differentiation due to drift in both geographic areas.Apart from the negative impact of drift on the genetic variation levels [1], the lack of functional connectivity makes local populations susceptible to alterations of their local habitat conditions [20].
The individuals of the two genetic groups exhibit contrasting habitat preferences (Table S1).Individuals from the eastern population mostly occupy natural water bodies, while individuals from the western population inhabit almost exclusively artificial cattle ponds (Table S1).Artificial ponds provide stable habitats, and their construction has been proved successful in boosting amphibian populations in the past [97][98][99][100][101]; their potential for conservation of crested newts has been demonstrated in situ [85].In the neighboring area of the Ioannina Basin, agricultural and livestock ponds have been identified as ideal habitats for newts to thrive [102].
Although artificial cattle ponds might seem stable enough (stable hydroperiod) and could potentially serve as effective breeding sites for numerous taxa, they simultaneously proved vulnerable to significant pressures locally.For example, fish introductions (releases) as well as high mechanical and chemical disturbances from cattle activity have been confirmed in the study area [25].Such pressures negatively affect artificial ponds by altering and downgrading habitat quality, with subsequent negative effects on aquatic communities [103].Introduced fish often prey on eggs or disrupt the spawning sites of newts and other aquatic breeders, reducing their evolutionary potential [103][104][105][106].In the neighboring Ioannina Basin, the presence of invasive fish in suburban ponds has been found to be responsible for the absence of the species from these sites [102].At the same time, the transition of traditional livestock practices (i.e., sheep farming) to the more "commercial" or financially profitable cattle farming has caused high mechanical and chemical disturbance of breeding sites in the study area [25].Other newt species, such as Triturus dobrogicus in Romania [107] and Lissotriton lantzi in Georgia [108], have been found to systematically avoid such problematic breeding environments.Within the study area, two robust and diverse breeding sites of the western genetic group, which served as a source to many others in a metapopulation context some 15 years ago [89], were extensively pressured to extinction due to the severe degradation of their habitat by both practices [25].
Certain management measures that would preserve and secure the quality and stability of the aquatic habitat of T. macedonicus, particularly in the western part of the national park, should implement fish removals, fencing of ponds, and pond restoration [109][110][111][112]. Additionally, pond construction could facilitate migration and dispersal movements, increasing the functional connectivity of local breeding sites [85].The maintenance of old artificial ponds and the establishment of new ones in the study area is expected to facilitate the metapopulation networks of Clade 2 and reinforce the presence of Clade 1 populations through natural colonization by individuals over a few years [85].At the same time, the aforementioned measures are anticipated to have a positive impact on agriculture and livestock farming [113].
Management actions should also be aimed at peripheral populations, such as sites 27 and 9 in the western group and site 4 in the eastern group, as these are found to be the most important sources of gene flow (Table S7).Peripheral populations are often characterized as genetic diversity hotspots [114][115][116][117] and have been the focus of numerous conservation programs [118][119][120].
Because T. macedonicus is considered as an "umbrella" species [121], the preservation and management of its habitats, apart from being beneficial for the breeding populations of the species, is expected to benefit a wide spectrum of aquatic breeders as well, which would thus preserve biodiversity in the hotspot that is Northern Pindos National Park.The coordination and cooperation of several public agencies are required to conserve each ESU and maintain the evolutionary continuum of the species [122,123].

Figure 1 .
Figure 1.(a) Geographic distribution of Triturus macedonicus (based on[31]).Location numbers correspond to TableS1.(b) Bayesian phylogeny of T. macedonicus haplotypes and five outgroups based on a 658 bp fragment of the mtDNA gene ND4.Numbers along internal branches represent posterior probability of Bayesian inference (left) and bootstrap support for maximum likelihood (middle) and neighbor joining (right).The coloration of clades follows[31].Newfound haplotypes (Tmac31, Tmac32, Tmac33, and Tmac34) are underlined.The examined T. macedonicus individuals in the national park fall within Clades 1 and 2. All individuals from the eastern part of the park belong to Clade 1 ("Northern Clade"), while all individuals from the western part belong to Clade 2 ("Zagori's Clade").The sudden expansion model appears to be supported for both Clades 1, 2, and Subclade 2A (FigureS1and TableS3).The expansion time, represented by τ, was longer for Clade 1 than for the others.The θ0 value was approximately zero for all clades, indicating the possible existence of a bottleneck event prior to expansion.Regarding the neutrality tests, the D value was negative and statistically significant only for Clade 1 and Sub-Clade 2A (p = 0.047 and p = 0.010, respectively), indicating population growth after a bottleneck event[58].The Fs value was negative and statistically significant only for Sub-Clade 2A (p = 0.000), indicating an excess number of haplotypes due to recent growth[59] or genetic hitchhiking[77].Clade 2 showed a negative D value, whereas a positive Fs value reflected a deficiency in haplotypes due to a recent bottleneck event[59].Finally, the available data (nine sequences, two haplotypes) were insufficient for conducting a demographic analysis of Sub-Clade 2B.

4 -
Within the national park, Clade 1 was represented by only three haplotypes.Similarly, Clade 2 was represented by three haplotypes.One dominant haplotype for each clade was identified: Tmac01 for Clade 1 and Tmac02 for Clade 2 (Figure2).

Figure 2 .
Figure 2. Spatial distribution of haplotype diversity in Northern Pindos National Park and the surrounding area.Location numbers correspond to TableS1.

Figure 3 .
Figure 3. Genetic population structure of Triturus macedonicus in Northern Pindos National Park: (a) STRUCTURE barplot showing two gene pools, comprising breeding sites from the eastern and western parts of the national park, respectively, with extremely low admixture between them.Vertical bars represent single individuals with their site of origin on the x-axis and the respective coefficients of membership on y-axis.(b) Geographic location and membership of each site in the two identified gene pools.Pie diagrams represent the average coefficient of membership of each site.E.P.: "Eastern Park" group; W.P.: "Western Park" group.

Table 1 .
p-distance (%) within and between clades and sub-clades of Triturus macedonicus and two sister Triturus species.Below diagonal: mean % p-distance between clades and sub-clades; Above diagonal: standard error estimates (after 1000 bootstrap replicates); Diagonal: % p-distance within clades and sub-clades with respective S.E. in parenthesis.

Table 2 .
Indices of genetic diversity per detected gene pool.Listed in order are: the sample size (n); mean number of alleles per locus (A); corrected allelic richness (Ar); number of private alleles (Ap); expected (He) and observed (Ho) heterozygosity, with corresponding ± SE in parentheses; and effective population size (Ne), with the 95% CI in parenthesis.