Genetic divergence and demography of pudu deer ( Pudu puda ) in five provinces of southern Chile, analyzed through latitudinal and longitudinal ranges

Pudu deer ( Pudu puda ) is endemic to the temperate rainforests of Chile. Genetic studies at different geographic scales for this species are required to better determine the genetic divergence within and among populations and their demography across the distribution range. These data can provide unique insights into the species or population status for conservation plans and decision-makers. We analyzed the mtDNA control region (CR) and cytochrome b (Cyt b) sequences of pudu deer in five provinces of southern Chile located at different latitudinal locations (Cautín, Valdivia, Osorno, Llanquihue and Chiloé Island) and three geographic areas within the studied provinces, representative of different longitudinal sites (Andes range, Central Valley and Coastal Range), to understand their genetic divergence and demography. The haplotype (H) and nucleotide (Π) diversities of CR and Cyt b ranged from 0.64286 to 0.98333 and from 0.00575 to 0.01022, respectively. CR diversity was significantly different among provinces, with Valdivia showing higher values than Llanquihue and Chiloé Island (H = 0.98333 vs. 0.64286–0.92727, P < 0.05). Cyt b variation also showed significant differences among provinces, particularly, among Cautín and Llanquihue (H = 1.000 vs. 0.222, P < 0.05). Genetic structuring among provinces was relatively high, as indicated by the F ST index (F ST = 0.41905). Clustering analysis indicated the presence of a distinctive cluster for Chiloé Island individuals. Fu’s F S and Tajima’s D based on CR revealed significant, negative deviations from equilibrium for Chiloé Island (D = -1.65898), Valdivia (Fs = -7.75335) and Llanquihue (Fs = -3.93267), suggesting population expansion in these provinces. Analysis at the longitudinal range showed significant differences among areas based on Π (P < 0.05), with the Andes range and Central Valley showing higher diversity than the Coastal Range. Neither population structuring (F ST = 0.01360, P > 0.05) nor distinctive clusters in the longitudinal range were observed. Fu’s Fs and Tajima’s D were negative and significant for the Coastal Range based on CR (Fs = -6.64752, P < 0.001) and Cyt b (D = -1.74110, P < 0.05), suggesting the existence of population expansion. Our results suggest that pudu deer in the analyzed provinces is a genetically structured species, which could be associated with reduced panmixia among populations. The genetic divergence pattern and the population expansion recorded are likely to be associated with past processes of recolonization after Pleistocene glaciation events.


Introduction
Pudu deer, Pudu puda (Molina, 1782), is endemic to southern South America and is characterized as being one of the smallest deer in the world due to its short shoulder height (30-40 cm) and small body weight (< 15 kg) (Jiménez 2010). This species is distributed in Chile from 35°10'S to 46°45'S, i.e., from the Maule to Aysen regions (Jiménez 2010), occupying an area of 128,278 km 2 according to the International Union for Conservation of Nature (IUCN), and is mainly located in Chile. Pudu deer characteristically inhabit pristine temperate rainforests, particularly in areas of dense understory growth and native bamboo thickets (Eldridge et al. 1987;Meier and Merino 2007), but they can also be found in disturbed and secondary forest habitats (Jiménez 2010). The current conservation status of pudu deer is near threatened according to the IUCN or vulnerable based on the threatened species list of the Chilean Ministry of Environment. This conservation status is related to different threats that appear to have affected the viability of the species. Among them, local threats have been identified to be linked to the expansion of human activities, such as forest loss and fragmentation, predation by domestic dogs, competition with exotic species, and poaching activities (Miller et al. 1973;Wemmer et al. 1998;Silva-Rodríguez et al. 2010;Silva-Rodríguez and Sieving 2012;Jiménez and Ramilo 2013). For example, loss or fragmentation of forest in Osorno and Llanquihue provinces from southern Chile is a matter of particular concern for pudu deer conservation, given that between 2006 and 2013, the area planted with exotic tree species increased significantly (+20.6% and +61.3%, respectively), and a large part of this growth occurred at the expense of native forest (CONAF-UACh (Corporación Nacional Forestal-Universidad Austral de Chile) 2014). This process is more accentuated in the central valleys located between the Andean and Coastal mountain ranges (Miranda et al. 2017). Thus, the conservation of pudu deer in southern Chile is a matter of public interest that requires combined efforts in several research areas, including genetic studies to assess the effect of different threats affecting their survival in their natural environment (Silva-Rodríguez et al. 2011).
Cautín, Valdivia, Osorno, Llanquihue and Chiloé Island are provinces in Southern Chile (39-44°S latitude) that contain a significant remnant of native temperate rainforest that covers a large proportion of each province. For example, Osorno Province contains approximately 42.9% of native forest, and in Llanquihue and Chiloé provinces they cover approximately 54.5% and 68.3%, respectively (CONAF-UACh (Corporación Nacional Forestal-Universidad Austral de Chile) 2014), mostly distributed in the Coastal and Andes mountain ranges in the first two provinces (Miranda et al. 2017). Both mountain ranges present a predominance of vegetational formations comprising temperate laurifoliar rainforests that include the Valdivian, North Patagonian and Subantarctic types (Villagrán and Hinojosa 2005). The temperate rainforest of Chile encompasses the Valdivian rainforest ecoregion, which has been listed among the most endangered ecoregions of the world and has a critical conservation status (Dinerstein et al. 1995;Olson and Dinerstein 1998;Miranda et al. 2017). In addition, the Valdivian rainforest ecosystem is considered a biodiversity hotspot and, therefore, a region of high conservation priority (Ormazábal et al. 1993;Myers et al. 2000;Smith-Ramírez 2004). Previous (Vanoli 1967) and recent (Pavez-Fox and Estay 2016;Colihueque et al. 2020) records support the presence of pudu deer in the aforementioned provinces, with higher probabilities of occurrence in the Coastal mountain range, where extensive coverage of native forest persists (Colihueque et al. 2020).
Inferences based on analyses of the geographic distribution of genetic diversity indicate that the Coastal mountain range in Osorno Province and other neighboring provinces, such as Cautín, Valdivia and Llanquihue, was a lowland refuge for vertebrates during the Pleistocene glaciation in the Last Glacial Maximum (LGM), which occurred approximately 23,000-18,000 years ago (Hulton et al. 2002), since this site is part of a refugial area identified in the Chilean Coastal Range north of 42°S (Sérsic et al. 2011). Paleobotanical evidence regarding the Quaternary history of the Chilean Coastal mountain range also supports this statement (Villagrán et al. 2019). In this area, the vertebrates persisted during this geological period to later recolonize the areas where the ice sheet had retreated (Rodríguez-Serrano et al. 2006;Himes et al. 2008;Victoriano et al. 2008). Genetic analyses of different vertebrate species distributed in this area (e.g., 39°S-42°S) have shown that they present private haplotypes and moderate to high levels of genetic diversity in comparison to those located along the Andean Mountain north-south axis of Patagonia, which comprises previously glaciated areas (Sérsic et al. 2011). Moreover, it should be noted that provinces in southern Chile contain several rivers that originate in the western Andean Mountains, which flow in a relatively straight line until reaching the Pacific Ocean, such as the Cautín, San Pedro and Bueno rivers (Errazuriz et al. 2000). These rivers constitute important geographic barriers due to their constant and large flow; therefore, it is likely that they may have promoted the genetic differentiation among populations of pudu deer. Moreover, it is likely that an important genetic divergence exists between pudu deer from Osorno and those from Chiloé Island because available studies support that continental populations of this cervid present a strong genetic divergence with respect to pudu deer inhabiting this island (ϕst = 0.75, Fuentes-Hurtado et al. 2011). This level of genetic divergence has been attributed to the effect of the Pleistocene glaciations that gave rise to a particular refugia on this island, whose populations were maintained in isolation from the continent by the Chacao channel . In fact, the estimated divergence time among Chiloé Island and continental populations of pudu deer would have occurred at least 217,500 years ago .
Phylogenetic studies performed on pudu deer populations from the northernmost distribution of southern Chile (36-42°S) using mtDNA control region and cytochrome b markers ) have indicated a significant genetic structure, with two main divergent clades corresponding to continental and Chiloé Island groups, and with lower genetic divergence among continental populations than between continental and island populations. Thus, these data suggested the existence of two different subspecies across the distribution range of the species; these subspecies have been attributed to ancient geographic isolation of the original population. In addition, a recent study that analyzed several populations that inhabit most areas of distribution of the species, including samples from the southernmost distribution (northern and southern Patagonia, 42-45°S) and Chiloé Island, confirmed this genetic structuring but revealed that both Patagonian groups are more related to the Chiloé Island group than to other continental groups based on the consistent subclusters that they formed (Cabello 2019). Despite the research progress on the phylogenetic relationships of pudu deer populations, until now, no detailed studies have addressed the genetic divergence of this cervid at a narrower geographic scale, for example, within and among provinces in southern Chile. This issue is interesting to address given that it is likely that the genetic structure of pudu deer populations of these sites could be affected by ancient paleoclimatic events or recent habitat modifications. For example, it is unknown whether the consequence of paleoclimatic events, such as the LGM, on the genetic structure of pudu deer populations differs across provinces or geographic areas. It is likely that presumably less affected populations from the Pacific coast in the Coastal Range should show a contrasting pattern of genetic variation. Supporting evidence indicates that these areas were relatively stable during this period (Villagrán and Armesto 2005;Villagrán et al. 2019); therefore, stronger signals of demographic equilibrium with increasing northern latitudes should be expected. In contrast, in areas with a pattern of stronger influences of glaciations (i.e., southernmost), which contain glaciated and recolonized areas, a predominance of derived haplotypes should be expected. Thus, by assessing the extent of gain or loss of genetic variation, the level of genetic divergence within and among populations and the identification of unique haplotypes could provide evidence to clarify the effect of paleoclimatic events on the modification of the genetic structure of pudu deer in southern Chile.
The objective of this study was to assess the level of genetic divergence and the demography of pudu deer in five provinces in southern Chile that are located at different latitudinal sites (from 39°S to 42°S) (Cautín, Valdivia, Osorno, Llanquihue and Chiloé Island) and three characteristic geographic areas within these provinces representative of distinctive longitudinal sites (Andes range, Central Valley (i.e., Longitudinal Valley) and Coastal Range), based on new and publicly available mtDNA control region (CR) and cytochrome b (Cyt b) sequences. We hypothesized that among provinces located at the extreme of the distribution, including those that inhabit an island, pudu deer populations must exhibit marked genetic divergence and reduced gene flow, given the existence of geographic barriers or the effect of Pleistocene glaciations that have modified the genetic structure of the populations. Moreover, pudu deer populations from the Coastal Range, are likely to have experienced an expansion process, since this site was a lowland refuge for vertebrates during the Pleistocene glaciations, to later recolonize the areas where the ice sheet retreated. In addition, populations located in this geographic area should exhibit a higher level of genetic diversity than populations located in other sites, for example, to those located in the Andean mountain range, due to its refugial nature. It is also possible that pudu deer population from some provinces may conform to a distinctive genetic cluster since these populations can present characteristic haplotype frequencies due to local selection pressure or population-specific haplotypes due to its relictual nature. Clarifying this issue is relevant for the sustainability of pudu deer in southern Chile, given that some populations may constitute an evolutionarily significant unit, whose identification is of major importance when conservation plans for threatened species are considered (Höglund 2009;Torres-Florez et al. 2018).

Sampling and DNA extraction
Pudu deer samples collected from Osorno (n = 6), Llanquihue (n = 8) and Chiloé (n = 6) provinces were analyzed (Table 1). Muscle tissue were obtained from individuals found dead along the highways that were donated by different institutions to our laboratory (Museo de Historia Natural, Purranque; Parque Nacional Puyehue, Osorno; Sitio Paleontológico de Pilauco, Osorno; Centro de Reproducción del Pudu, Osorno). Tissue samples were then fixed in 80% ethanol. Moreover, preserved blood samples taken from live individuals at the Wildlife Rehabilitation Center (Centro de Rehabilitación de Fauna Silvestre) of the Universidad San Sebastián, Puerto Montt, Chile, were also used. DNA was extracted from the fixed muscle or blood tissue using the phenol-chloroform method, as described in Taggart et al. (1992). Extracts were standardized at 100 ng/μL using Tris-EDTA buffer at pH 8.0. To increase the sample size aimed at robustness of the genetic analysis, CR (n = 27) and Cyt b (n = 11) sequences previously published by Fuentes-Hurtado et al. (2011) were included (Table 1). These sequences represent different haplotypes observed for pudu deer populations in the analyzed areas, particularly, the Y02-Y25 haplotypes for CR and the A-C and E-I haplotypes for Cyt b.

Population differentiation and genetic diversity
We used DNASP version 5.1 software to estimate genetic variation within populations based on the number of haplotypes (h), polymorphic sites (s), haplotype diversity (H), nucleotide diversity (Π) and average number of nucleotide differences (k) (Librado and Rozas 2009). Coalescent analysis based on a neutral infinite-site model and assuming a large constant population size (Hudson 1990), with theta per gene being estimated from the data, no recombination and 1,000 replicate settings, was also performed using the same software. This approach allowed us to estimate the average number and variance (95% confidence interval) of expected haplotypes (Neh) to assess whether sample sizes were sufficient to capture the genetic diversity. We estimated pairwise F ST values based on the Kimura 2P distance method (Nei and Kumar 2000) and haplotype frequencies that were implemented in ARLEQUIN version 3.1 (Excoffier and Lischer 2010) to quantify the extent of genetic differentiation among populations. For pairwise F ST , the significance was tested using 1,023 permutations between populations with a level of significance of α = 0.05. The heterogeneity of haplotype frequencies was assessed using the chi-square homogeneity test (Sokal and Rohlf 2009). In this test the null hypothesis of equality in haplotype frequencies within each geographic area was tested. We performed the Welch t-test for unequal variances to test for statistical significance of comparisons among study areas based on the H and Π indices (Ruxton 2006). Analysis of molecular variance (AMOVA) (Excoffier et al. 1992) implemented in the program ARLEQUIN version 3.1 (Excoffier and Lischer 2010) was used to assess the geographic pattern of population subdivision. In AMOVA the correlation of haplotype frequency was used as an F-statistic analogue at various hierarchical levels. To carry out this analysis, samples were separated into five groups according to their sampling province, which have different latitudinal locations from north to south as follows: Cautín Province (CAU, 39°0'0"S, 72°30'0"W), Valdivia Province (VAL, 39°45'0"S, 72°30'0"W), Osorno Province (OSO, 40°35'0"S, 73°10'0"W), Llanquihue Province (LLA, 41°20'0"S, 72°50'0"W) and Chiloé Island Province (CHI, 42°30'0"S, 74°0'0"W) (Fig. 1a). We selected these provinces because there was a sufficient sample size to carry out the genetic analyses in these geographic areas (from 5 to 16). In addition, to test the possible existence of genetic differentiation between populations located at different longitudinal locations, samples from different provinces, except those of Chiloé Island Province, were grouped as belonging to the Coastal Range (COA), Central Valley or Longitudinal Valley (CEN) and Andes Mountain range (AND) based on their geographic coordinates (Fig. 1b).  Table 1.

Gene flow among provinces
The number of migrants per generation (Nm) among provinces was estimated from F ST values using ARLEQUIN version 3.1 (Excoffier and Lischer 2010). This estimation assumes that two populations of size N drawn from a large pool of populations exchange a fraction m of migrants each generation and that the mutation rate u is negligible compared to the migration rate m.

Population clustering
We used a clustering analysis through kMeans version 1.1 (Meirmans and Van Tienderen 2004) to split objects into groups by minimizing the within-group genetic diversity while maximizing it among-groups. For this analysis, the input file comprised a matrix of allele frequencies coded as binary data. An a priori assigned number of groups (k) from 1 to 10 was tested with 50,000 permutations. We used the Caliński-Harabasz pseudo-F statistic (Caliński and Harabasz 1974) to determine the optimal number of clusters. In the pseudo-F statistic, the optimal clustering is the one with the highest value. Simulated annealing as a clustering algorithm was used based on the Monte Carlo Markov chain (MCMC), which prevents local optima becoming stuck.

Population dynamics
We tested for patterns of demographic history with Tajima's D (Tajima 1989) and Fu's Fs (Fu 1997) tests of neutrality. Tajima's D test is based on the infinite-site model without recombination and is appropriate for short DNA sequences. The significance of the D statistic is tested by generating random samples under the hypothesis of selective neutrality and population equilibrium using a coalescent simulation algorithm adapted from Hudson (1990). Significant D values can be due to factors other than selective effects, such as population expansion, bottleneck, or heterogeneity of mutation rates (Tajima 1996). Fu's Fs test is based on an infinite-site model without recombination. This test evaluates the probability of observing a random neutral sample with a number of alleles that is equal to or smaller than the observed number of alleles. When the Fs statistic is large and negative, one can infer demographic expansion. This test has been shown to be especially sensitive to departure from population equilibrium, as in the case of population expansion. Tajima's D and Fu's Fs tests were run in ARLEQUIN version 3.1 (Excoffier and Lischer 2010).

Analysis at the latitudinal range CR and Cyt b haplotypes and nucleotide polymorphism
Thirty-three CR haplotypes were detected in 47 pudu deer individuals across the five provinces analyzed ( Table 2). The coalescent analysis indicated that populations approached the estimated total haplotype number relatively well since, with the exception of Valdivia Province, all provinces fell within a 95% confidence interval (Table 3). Thus, this analysis indicated that most haplotypes were captured. Thirty-six of 626 nucleotide positions were variable (5.75%), including 27 parsimony informative sites and nine singleton variable sites. The number of haplotypes per province ranged from 4 to 14, with a mean of 7.2. Haplotype diversity (H) ranged from 0.64286 to 0.98333 among provinces, and nucleotide diversity (Π) ranged from 0.00575 to 0.01022 (Table 3). Haplotype diversity (H) and nucleotide diver-   Columns and bar errors represent the means and standard deviation of the means, respectively. Significant differences between means were calculated according to Welch´s t-test. *P < 0.05, NS = Non significant difference. sity (Π) were significantly distinct among provinces especially among Valdivia and Llanquihue plus Chiloé provinces in the former (0.98333 versus 0.64286-0.92727, Welch´s t-test, P < 0.05) (Fig. 2a), and between Cautín and the other four provinces in the latter (0.01022 versus 0.00575-0.00797, Welch´s t-test, P < 0.05) (Fig. 2b).
Twelve Cyt b haplotypes were detected in 30 pudu deer individuals across the five provinces analyzed ( Table 2). The coalescent analysis indicated that populations approached the estimated total haplotype number relatively well, since all provinces fell within a 95% confidence interval (Table 3). Fourteen of 743 nucleotide positions were variable (1.88%), including six parsimony informative sites and eight singleton variable sites. Thus, the Cyt b sequence was approximately 70% less variable than the CR sequence. The number of haplotypes per province ranged from 2 to 4, with a mean of 3.4. Haplotype diversity (H) ranged from 0.222 to 1.000 among provinces, and nucleotide diversity (Π) ranged from 0.00030 to 0.00336 (Table 3). Haplotype diversity (H) and nucleotide diversity (Π) were significantly distinct among some provinces, especially among Cautín and Llanquihue in the former (1.000 versus 0.222, Welch´s t-test, P < 0.05) and between Valdivia and Llanquihue in the latter (0.00336 versus 0.00030, Welch´s t-test, P < 0.05).

Population genetic structure
Global AMOVA results as a weighted average over loci, conducted by partitioning variation among and within populations, revealed that CR haplotype variation was attributed mainly to within-population differences (58.1%) but also to differences  among populations (41.9%). In fact, the F ST index reached a value of 0.41905 (P < 0.001). Population pairwise F ST values ranged from 0.03617 to 0.68651, with seven out of ten pairwise comparisons showing significant differences (P < 0.05) ( Table 6). This pattern of genetic differentiation was clearly observed between Cautín and Valdivia and Llanquihue (F ST = 0.15385-0.16810) and among Chiloé Island compared to the other four continental provinces (F ST = 0.59004-0.68651). The clustering analysis based on allelic frequencies of CR revealed the existence of either two (K = 2, pseudo-F = 29.620) or three (K = 3, pseudo-F = 24.695) groups as the optimal number of clusters. In both simulations, individuals from Cautín, Valdivia, Osorno and Llanquihue provinces overlapped into one or two clusters, in contrast to those of Chiloé Island Province, which mostly formed a distinctive single group (Fig. 7a).

Gene flow among provinces
Pairwise estimate of Nm for CR haploytpe data ranged from 0.22832 to 13.32296. The number of migrants per generation was high among Osorno, Valdivia and Llanquihue provinces (Nm = 5.18197-13.32296), medium between Cautín and Osorno plus Valdivia and Llanquihue (Nm = 2.47437-3.66488), and low for Chiloé Island compared to the other four continental provinces (Nm = 0.2283-0.34740) ( Table 6). Thus, patterns in pairwise Nm values were consistent with population structure based on AMOVA and clustering analysis.

Evolutionary neutrality and population expansion
Tajima's D estimate for CR was negative and significant in one province that corresponds to Chiloé (D = -1.65898, P < 0.05). In contrast, Tajima's D statistic of other provinces was consistent with population equilibrium (Table 3). Fu's Fs was negative and significant for Valdivia (Fs = -7.75335, P < 0.001) and Llanquihue (Fs = -3.93267, P < 0.05) provinces, which occurs when an excess of rare haplotypes is present and suggests that either population expansion or genetic hitchhiking has occurred. These results are consistent with deviation from neutrality due to either selection or population expansion in some provinces. Tajima's D and Fu's Fs estimations based on Cyt b were negative and not significant for all provinces (Table 3). These results are consistent with population equilibrium.

Analysis at the longitudinal range
Analysis at the longitudinal range based on CR revealed a number of haplotypes per area ranging from 6 to 11 in 32 analyzed individuals (Table 4), with a mean of 7.3. Haplotype diversity (H) ranged from 0.955 to 1.000 among areas and nucleotide diversity (Π) ranged from 0.006513 to 0.008307 (Table 5). Only nucleotide diversity (Π) was significantly different among areas, with the Andes range and Central Valley showing higher diversity than the Coastal Range (Welch´s t-test, P < 0.05) (Fig.  2d). Only four out of 22 haplotypes were shared among the Andes range, Central Valley and Coastal Range. In fact, the Y27 haplotype was shared among the Andes range and Central Valley, as were the Y17 haplotype between the Andes range and Coastal Range, and the Y05 and Y10 haplotypes among the Central Valley and Coastal range (Fig. 3b). Haplotype frequency were heterogeneous only for Coastal range (chi-square test: χ 2 = 23.49, df = 10, P < 0.01). In fact, the Y17 and Y32 haplotypes showed high frequencies in this area (21.4 and 14.3%, respectively) (Fig. 4b). In the case of Cyt b, haplotypes per area ranging from 4 to 6 in 22 analyzed individuals were observed (Table 4). Haplotype diversity was significantly higher in the Andes range than in the Coastal Range (0.803 vs. 0.533, Welch´s t-test, P < 0.05) (Table 5). Moreover, the haplotype frequencies were heterogeneous for Andes (chi-square test: χ 2 = 58.37, df = 5, P < 0.001) and Coastal (chi-square test: χ 2 = 108, df = 3, P < 0.001) ranges. In fact, the A haplotype showed a high frequency in both provinces (41.7 and 70%, respectively) (Fig. 6b). The clustering analysis based on allelic frequencies of CR revealed the existence of either two (K = 2, pseudo-F = 8.802) or three (K = 3, pseudo-F = 8.011) groups as the optimal number of clusters. In both simulations, individuals from the Andes range, Central Valley and Coastal Range overlapped into two or three clusters (Fig. 7b).
Global AMOVA results as a weighted average over loci revealed that CR haplotype variation was attributed mostly to within area differences (98.6%) and scarcely to differences between areas (1.3%). In fact, the F ST index was not significant and reached a value of only 0.01360 (P > 0.05). Pairwise F ST values based on CR haplotype frequencies ranged from -0.01741 to 0.03434, with no significant differences   -4.230 observed in any comparison (P > 0.05). Pairwise estimate of Nm for CR haplotype data ranged from 14.05932 to infinity, revealing the occurrence of an important number of migrants per generation among the Andes range, Central Valley and Coastal Range. Tajima's D was negative in all areas, but it was not significant (D = from -0.02394 to -0.97905, P > 0.05) (Table 5). Thus, these results were consistent with population equilibrium. However, Fu's Fs was negative and highly significant for the Coastal Range (Fs = -6.64752, P < 0.001), which suggests that in this area population expansion may have occurred. Tajima's D estimated based on Cyt b showed a similar result as CR given that the D value was negative and significant for the Coastal Range (D = -1.74110, P < 0.05) ( Table 5).

Discussion
Analysis of the CR and Cyt b mitochondrial sequences from pudu deer populations of five provinces from southern Chile revealed different levels of genetic variation, haplotype heterogeneity, a large proportion of unique haplotypes, and strong genetic structuring and clustering of individuals in some provinces. The most likely explanation for our results is that pudu deer from the study area have a low to moderate level of gene flow among populations, which suggests the occurrence of reduced panmixia across the collection range. This genetic pattern is in accordance with our predictions of marked genetic divergence among pudu deer populations, especially those located at the extremes of the distribution or those that inhabit islands, due to the existence of geographic barriers that may have promoted genetic divergence. In fact, individuals from Chiloé Island that are separated from the continental populations by the Chacao channel formed clear single groups in the clustering analysis in concordance with the existence of such a geographic barrier, which should have reduced gene flow with the other continental analyzed populations. In this regard, the separation among continental and island populations of pudu deer, estimated at least 217,500 year ago , would have been key event that has promoted genetic divergence of pudu deer that inhabit Chiloé Island. Pudu deer of other sites, such as those in Cautín Province, also showed a marked divergence with respect to other analyzed provinces (F ST = 0.15385-0.59004) and a reduced number of migrants per generation (Nm = 0.34740-3.66488). Since this province is located in the northernmost part of the range of the species distribution, it is likely that geographic barriers, such as rivers, and/or distance constrain gene flow and, therefore, facilitate an increase in genetic divergence. All haplotypes of the CR marker found in Cautín Province were private, a finding that provides further support for a divergent evolutionary process for pudu deer located at this site. Although further analysis is required, genetic differentiation between provinces recorded for pudu deer in our study could be also related to their habitat loss and fragmentation in southern Chile as a consequence of the expansion of human activities (Echeverría et al. 2006(Echeverría et al. , 2007Miranda et al. 2017). This could be so given that habitat fragmentation contributes to genetic differentiation between populations, since fragmented populations are expected to see more differentiation because migration between them is impaired. Of note is that our previous interpretations should be taken with caution because, as discussed below, the modifications of the genetic structure of pudu deer could also be a consequence of Pleistocene glaciations that affected the landscape in southern Chile in the past. Thus, the modifications of the genetic structure of pudu deer are likely to be a more complex process in which different factors have intervened.
In addition, our results also suggest that some provinces, particularly Valdivia, Llanquihue and Chiloé Island, underwent population expansion given that they showed significant negative Tajima's D and Fu's Fs values. This expansion process could explain the shared haplotypes among neighboring provinces, such as Valdivia, Osorno and Llanquihue, where the Y05 and Y17 CR haplotypes were highly frequent. In the case of the Cyt b marker, although it did not present any significant population disequilibrium, the occurrence of shared haplotypes for this marker among all provinces analyzed, particularly the A and I haplotypes, suggests that population expansion for pudu deer could be a more extensive process. Of note is that the analysis at the longitudinal range supports that populations located in the Coastal Range would be an important source of individuals for colonization, given that both markers showed significant population disequilibrium in this geographic area. These results are in accordance with the statement that this area was a lowland refuge for vertebrates during the Pleistocene glaciations (Sérsic et al. 2011), where the vertebrates persisted during this geological period to later recolonize the areas where the ice sheet had retreated (Rodríguez-Serrano et al. 2006;Himes et al. 2008;Victoriano et al. 2008). In addition, we also observed that pudu deer of Coastal Range presented lower or similar genetic diversity, measured as haplotype diversity and nucleotide diversity, than those located in the Andes range. It has been observed that vertebrate populations distributed in the Andes range (e.g., 39°S-42°S) tend to exhibit lower genetic diversity, while the opposite is observed for populations distributed in the Coastal Range (Sérsic et al. 2011). This pattern of genetic variation seems to be related to the recolonization process that occurred in the Andes Mountains after the Pleistocene glaciations, where founder effects or selection probably changed the genetic structure of populations. The modifications of the genetic structure as a result of Pleistocene glaciations in the Andes range have been well documented in other vertebrates that inhabit this geographic area, such as the lizard Liolaemus tenuis from 38°S to the southern limit of its distribution (Victoriano et al. 2008), the marsupial Dromiciops gliroides located between 40 and 43°S (Himes et al. 2008) and the lizard Liolaemus pictus (Vera-Escalona et al. 2012). In fact, for L. tenuis, there is a general pattern, whereby populations located in the mesomorphic zone have a higher genetic diversity than those inhabiting in the periglacial and glacial zones. In the case of D. gliroides, southern Clade C, strongly affected by glacial recession, showed lower nucleotide diversity than northern Clade A. Similarly, L. pictus in the northern group (37°S-39°S) was heterogeneous, with well differentiated and deeply divergent internal subgroups, whereas the southern group (39°S-41°S) was more homogeneous. Since Andes range populations of pudu deer analyzed in this work presented high genetic variation, an alternative explanation is required. Among them a mixing origin of this population cannot be ruled out, which implies the occurrence of various recolonization processes with different source populations. This hypothesis is not unlikely because the haplotype composition of CR observed in the Andes range included haplotypes proper to Cautín (Y02, Y08), Valdivia (Y09) and Osorno (Y15, Y26, Y28) provinces. The multiple sources origin hypothesis is not new for vertebrates that inhabit the western Andes range in southern Chile. For example, Vidal et al. (2012) by using the Cyt b marker, found evidence that the lizard Liolemus pictus that inhabits this area may have experienced such a process during Pleistocene glaciations, given their geographically differentiated gene pools. For this lizard, postglacial colonization could have occurred from sources located in ice-free sectors of the Longitudinal Valley adjacent to the Antillanca locality and not from a single western refuge.
A number of studies have indicated that Pleistocene glaciations modified the landscape in southern Chile, which affected the species distribution and created habitat fragmentation of areas that species used as refuges (Lessa et al. 2010;Vera-Escalona et al. 2012;Vidal et al. 2012;Marín et al. 2013). From a population genetic point of view, this glaciation process may have promoted the genetic divergence among populations, with some characteristics of this process being the existence of populations with private haplotypes, haplotype heterogeneity or a reduction in genetic variation as a consequence of a decrease in population effective sizes. Thus, the genetic differentiation between pudu deer from Chilóe Island and surrounding provinces in continental areas fit well with this divergence pattern. In fact, we determined that pudu deer from Chilóe Island do not share CR haplotypes with any continental populations. Moreover, they showed significantly lower haplotype diversity than other analyzed provinces (e.g., with Valdivia Province), which may reflect past processes of population size reduction. Chilóe Island Province also showed significant haplotype heterogeneity, a process that may be related to reduction in population size, since this provokes genetic drift, causing genetic differentiation in populations. In other vertebrates that inhabit this island a similar divergence pattern has been observed. For example, the lizard L. pictus conforms to a well-supported subclade and is composed exclusively of haplotypes proper from Chiloé Island, although the genetic distances between haplotypes from this island and the surround-ing mainland populations are low compared to those from the northern clade (Vera-Escalona et al. 2012). The absence of shared haplotypes between populations from the mainland and Chiloé Island, which give rise to divergent haplogroups in L. pictus, might reflect ancient isolation, much older than the LGM (Vidal et al. 2012). Of note is that further niche simulation analyses support that this pattern of divergence in L. pictus may be related mostly to the LGM (Vera-Escalona et al. 2012), given that at that age a greater suitability of the habitat in the extreme north of the island of Chiloé may have existed, which probably acted as a refuge for different vertebrates.
Previous reports indicate that pudu deer from Chiloé Island have marked genetic divergence with respect to continental populations of this species, especially those from the northernmost distribution Cabello 2019). For example, Fuentes-Hurtado et al. (2011) using the CR marker, found that there were higher genetic distances among continental populations and Chiloé Island populations (2.2%) than between continental pudu deer populations (0.9%). Our clustering analysis results are in accordance with this divergence pattern, since almost all pudu deer individuals from Chiloé Island formed a consistent single cluster when using a clustering algorithm based on a Monte Carlo Markov chain. Moreover, pudu deer from this site showed a marked divergence with respect to those inhabiting the other analyzed provinces (F ST = 0.59004-0.68651). Therefore, these results support that they may constitute an evolutionarily significant unit and could be treated as a separate management unit. This finding has conservation implications because pudu deer in Chile is classified as a threatened species according to the Chilean Ministry of Environment, and the identification of an evolutionarily significant unit is of major importance when conservation plans for threatened species are considered (Torres-Florez et al. 2018). Currently, this issue is of public interest because in Chiloé Island Province, pudu deer faces several local threats linked to the expansion of human activities, such as forest loss and fragmentation, predation by domestic dogs and competition with exotic species (Miller et al. 1973;Wemmer et al. 1998;Silva-Rodríguez et al. 2010Silva-Rodríguez and Sieving 2012;Jiménez and Ramilo 2013). Pudu deer from the other analyzed provinces, could be also treated as different management units because they presented a large ratio of unique haplotypes (e.g., from 0.66 to 1.0 for CR), which may reflect the existence of local adaptation, a process that may arise in response to localized natural selection pressures. In this regard, we hope that the genetic divergence data presented in this study will help to improve the conservation status of pudu deer in southern Chile.

Conclusion
Genetic and demographic analyses of the mtDNA control region and cytochrome b sequences carried out for pudu deer of five provinces located at different latitudinal sites in southern Chile, indicate the existence of different levels of genetic variation, haplotype heterogeneity, a large proportion of unique haplotypes, strong genetic structuring, clustering of individuals in a province (i.e., Chiloé Island) and a signature of population expansion in some provinces (i.e., Valdivia, Llanquihue and Chiloé Island). These results indicate that pudu deer have a low to moderate level of gene flow, especially among distant continental provinces or between continental provinces and Chiloé Island, which suggests the occurrence of reduced panmixia. This genetic pattern could be related to the existence of geographic barriers and/ or paleoclimatic events that may have promoted genetic divergence, such as the Pleistocene glaciations that modified the landscape in southern Chile. Analysis at the longitudinal range within the studied provinces revealed haplotype variation attributed mostly to within areas differences and the existence of a population expansion process in the Coastal mountain range. The last finding supports the statement that this area was a lowland refuge for vertebrates during the Pleistocene glaciations, where the vertebrates persisted during this geological period to later recolonize the areas where the ice sheet retreated. Since pudu deer in Chile is a threatened species, our results have conservation implications because populations of some analyzed geographic areas may constitute significant evolutionary units, which is of major importance when conservation plans are considered.