Introduction

The Pampas deer (Ozotoceros bezoarticus) was once a widespread and abundant species occupying a wide range of open habitats, including grasslands, pampas in Argentina and the Brazilian savanna known as the Cerrado (^{Cabrera 1943}; ^{Jackson 1987}; ^{Merino et al. 1997}; ^{González et al. 2002}; ^{Weber and González 2003}; ^{González et al. 2010}). Historical records showed that this species had a wide range distribution in southeastern South America (from -5° to -41° S), reported by several naturalists mentioned in their records as Charles ^{Darwin (Darwin 1860}). The area encompassed by these habitats has been dramatically reduced to less than 2 % by human activities such as agriculture, urbanization and poaching of that present in 1900 (^{González et al. 2023}). Currently, Pampas deer populations are generally small and isolated (^{Jackson and Langguth 1987}; ^{Pinder 1994}; ^{González et al. 2010}, Figure 1).

Mitochondrial DNA has been the most widely used tool for reconstructing population and species histories, presumably because it is relatively easy to amplify, typically non-recombining, supposedly nearly neutral, and highly variable between and within species (^{Taberlet 1996}; ^{Avise 1998}). Furthermore, this molecular marker of maternal inheritance is useful for the genetic analysis of populations in fragmented habitats, such as those of the Neotropical deer (^{Avise 1992}, ¹⁹⁹⁵). Most deer species usually have: 1) asymmetric genetic flow and dispersal rate, females frequently being phylopatric. 2) Females and fawns are spatial associated. 3) A strong maternal lineages structuration resulting in demographic autonomy among populations in an ecological scale (^{González et al. 1998}; ^{Márquez et al. 2006}; ^{González et al. 2010}). Historical population sizes based on mtDNA control region sequences were estimated to be several magnitudes larger than present day estimates. Gene flow patterns also showed high levels of genetic differentiation among isolated populations using representative samples of Pampas deer from throughout its geographic range. In that study, ^{González et al. (1998)} identified five conservation genetic units for the six localities surveyed: two in Brazil (Emas National Park and Pantanal da Nhecolandia), two in Uruguay (Salto and Rocha), and one in Argentina (composed of two populations: Samborombón Bay and San Luis).

The comparative craniometric analysis of the different Pampas deer populations revealed that differentiation was concordant with the levels of genetic differentiation found with mitochondrial markers. Furthermore, two new subspecies were recognized in the Uruguayan northwestern (O. b. arerunguaensis, Salto Department) and eastern (O. b. uruguayensis, Rocha Department) grasslands (^{González et al. 2002}).

This species is currently considered Near Threatened (NT) in the global IUCN Red List (^{González et al. 2016}). In their southern range of South America, the most threatened populations occur in Argentina, Bolivia, Paraguay, and Uruguay with fewer than 2,500 mature individuals. This is reflected by the Red List categories of the subspecies: O. b. celer, -Argentina-: Endangered [EN B1ab(iii)]; O. b. arerunguaensis -Uruguay- [CR B1ab(iii); O. b. uruguayensis -Uruguay- [CR B1ab(iii)]; O. b. bezoarticus -Brazil- (DD); O. b. leucogaster - Argentina, Bolivia, Brazil; Paraguay- Near Threatened (NT) (^{González et al. 2016}). The species is also included in CITES Appendix I (^{Giménez-Dixon 1987}).

Three decades ago, we conducted the first molecular genetic study of the Pampas deer based on representative samples from throughout their geographic range (^{González et al. 1998}). We aimed to deduce genetic units for conservation (^{Moritz 1995}<DB>) and to better understand the effect of habitat fragmentation on gene flow and genetic variation. In this study, we evaluated the effects of habitat fragmentation on gene flow among eight wild Pampas deer populations and one from the captive breeding centre Estación de Cria de Fauna Autóctona (ECFA) in Uruguay. We examined DNA sequences from three mitochondrial markers: the control region (D-loop), Cytochrome b (Cytb), and Cytochrome Oxidase I (COI) to determine levels of genetic differentiation among isolated populations. Additionally, we compared the resolution of the different mitochondrial markers to elucidate the phylogenetic and phylogeographic patterns of the species and to define Evolutionary Significant Units (ESU`s).

Almost thirty years later, we increased the sample size from the previous locations and incorporated new sites from its wide geographic range to revisit the genetic characterization of Pampas deer by implementing additional mitochondrial markers. Our results will be providing a comprehensive approach for understanding the current genetic status and the future viability of the species.

Materials and methods

Sample collection. We analyzed 164 Pampas deer samples from eight geographic localities across the species range from Argentina, Bolivia, Brazil, and Uruguay and one captive population in Uruguay, ECFA (Figure 1; Supplementary Material Table 1). This captive stock was founded in 1980 with 10 individuals from the wild population from Salto.

DNA extraction and PCR amplification. Genomic DNA was extracted from tissue samples (50 mg, see details on Supplementary material) following ^{González et al. (1998}) protocol. DNA from fresh feces that were stored in ethanol and refrigerated was extracted using the commercial QIAamp® Fast DNA Stool Mini Kit (QIAGEN, Hilden, Germany) following manufacturer’s instructions. All PCRs reactions were carried out in an automatic TC 9639 Thermal Cycler (Benchmark Scientific) in a mixture of final volume of 15 uL containing 3 ng/ul of sample genomic DNA, 7.5 μL of Immomix^TM mastermix (Bioline), 0.5 μL of each primer (10 μM, Table 1) and ultra-pure water.

The profile consisted of an initial denaturation at 95°C for 10 min, followed by 35 cycles of denaturation at 95°C for 1 min, an annealing step for 2 min (Table 1), an extension at 72°C for 1.5 min, and a final extension at 72°C for 7 min. Positive and negative controls were included in every PCR to check for contamination in different experiments.

PCR products were purified with DNA Clean and Concentrator™ (Zymo Research™) kit and diluted to an equal final concentration using a Nanodrop 1000™ Spectrophotometer (Thermo Fisher Scientific). The amplicons were sequenced by the Sanger method on an automatic sequencer ABI 3130 (Applied Biosystems) at the Pasteur Institute (Montevideo, Uruguay) and on an automated ABI 3730xl System Sequencer (MACROGEN Inc., Korea).

Bioinformatic Analysis. Sequences were aligned and edited in MEGA11 (Molecular Evolutionary Genetics Analysis software version 11,^{Tamura et al. 2021}) and compared with the nucleotide database available in the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/) using BLAST (Basic Local Alignment Search Tool) utility. We calculated the number of paired differences for each mitochondrial marker. Diversity indexes as haplotype, nucleotide diversities and the number of polymorphic sites were calculated using DnaSP v. 6.12.3 (^{Rozas et al. 2017}<DB>). The evolutionary genetic distances among haplotypes were determined by the ^{Kimura 2-Parameter distance (Kimura 1980}). Nucleotide sequence data were analyzed using maximum parsimony (MP) and neighbor joining (NJ) in the software MEGA11 (^{Tamura et al. 2021}), as well the cluster support trees were evaluated with 1000 bootstrap pseudo-replicates.

Haplotype networks. To evaluate the evolutionary relationships among haplotypes and their geographic distribution, we constructed two haplotype networks (from Cytb and D-loop sequences) using the median-joining network approach (^{Bandelt et al. 1999}) implemented on PopART 1.7 (^{Leigh and Bryant 2015}).

Patterns of geographic subdivision and gene flow. We used AMOVA (Analysis of Molecular Variance) to deduce the significance of geographic divisions among local and regional population groupings (^{Excoffier et al. 1992}; ^{Schneider et al. 2000}). We calculated the fixation index within populations (ΦST) and among populations (ΦSC), as well as within groups and between groups (ΦCT). We also estimated the average number of migrants per generation (using ΦST estimates between populations) to measure the degree of isolation of populations or the degree of subdivision among populations. The significance of F-statistic analogues was evaluated by 1,000 random permutations of sequences among populations. We experimented with various grouping of populations as suggested by the analysis of DNA sequences and population trees. The groupings which maximized values of ΦCT and were significantly different from random distributions of individuals were assumed to be the most probable geographical subdivisions (^{Excoffier et al. 1992}; ^{Schneider et al. 2000}).

Gene flow within and among regions was approximated as Nm, the number of female migrants occurring between populations units per generation and was estimated using the expression Fst = 1/ (1+2Nm) where N is the female effective population size and m is the female migration rate (^{Slatkin 1987}; ¹⁹⁹³; ^{Baker et al. 1994}). We used pairwise estimates of ΦST as surrogates for FST among regional groupings of populations (e. g.^{Stanly et al. 1996}).

Phylogenetics and evolutionary rate estimation of three mtDNA regions. The phylogenetic relationships and evolutionary rates based on three markers (D-loop, Cytb and COI) were estimated using the Neighbour Joining algorithm. Trees were drawn to scale with the length of the branches in the same units as the evolutionary distance used to infer the phylogenetic tree. The genetic distances were calculated using the Kimura 2-Parameters algorithm and the units are the number of base substitutions per site.

We used the divergence time between the Pampas deer O. bezoarticus and the grey brocket deer Subulo gouazoubira based on ^{Duarte et al. (2008}), to estimate the evolutionary rate of the three mitochondrial markers Cytb, COI, and D-loop.

Figure 1 Distribution map of Pampas deer populations. Grey shadow indicated the presumed past distribution of the species. Yellow shaded areas represent the current distribution. Numbers indicated the 8 localities sampled: 1. Rocha; 2. Salto; 3. Samborombón Bay; 4. Santa Fé; 5. San Luis; 6. Emas National Park; 7. Pantanal da Nhecolandia; and 8. Paraná. Asterisk (*) shows ECFA (Estación de Cría y Fauna Autóctona) locales in Maldonado, Uruguay. 

The evolutionary rate of Cytb was estimated using sequences deposited in GenBank by ^{Duarte et al. (2008}; accession numbers: DQ789173-DQ789231) together with additional sequences retrieved from GenBank and employed by the authors in the phylogenetic analysis. To calculate the evolutionary rate of the D-loop mitochondrial region, we used Pampas deer sequences and grey brocket deer available in GenBank (Accession numbers: AF012556-AFO12602).

Results

We sequenced 138 individuals of the mitochondrial control region (an amplicon of 603 bp). We analyzed a 423 bp fragment and found 87 different haplotypes in the eight localities defined by base-pair substitutions (Table 2). As was reported previously, this species has a polymorphic dinucleotide TA repeat sequence within this amplicon that had four to eight tandem repeats beginning at nucleotide position 186, with position 1 as the first nucleotide of our control region sequence (^{González et al. 1998}). Because the same allele sizes were found in divergent sequences from geographically distant populations (e. g. Emas and Rocha) this region has a high degree of homoplasy and we excluded the tandem repeat region from the analysis, leaving 423 base pairs of DNA sequence to be analyzed.

Patterns of geographical subdivision and gene flow revealed by D-loop region. We analyzed samples from nine geographic locations that showed high levels of polymorphism (Table 3). The highest value obtained for the average number of pairwise differences within population (PiX) was in the Paraná samples and the lowest in Samborombón Bay.

Geographic distribution of control region sequences. The control region haplotypes were perfectly segregated as no locality shared haplotypes except for the Salto population that shared haplotypes with the individuals from the Breeding center ECFA. Sequences from the same locality tend to be clustered together in minimum spanning networks (Figure 2). The nine groups according to the haplotype’s geographic distribution and the taxonomic criteria were arrange in six populations (Groups: 1 Salto and EFCA, 2 Rocha, 3: Paraná. 4: Emas, 5: Pantanal da Nhecolandia and North Argentina, and 6: Samborombón Bay and San Luıs (Table 4). This is because the haplotype diversity present in each population, where most of the haplotypes found for each, are unique and not shared with other populations.

The AMOVA results showed high levels of differentiation among populations within groups (ΦST = 0.439), among groups (ΦSC = 0.240), and as well as within populations (ΦCT = 0.262; Table 5). The Argentinean populations (Samborombón Bay and San Luis) clearly are more closely related to each other than are those from Brazil or Uruguay.

Genetic differences among and within populations. The mtDNA control region showed a rapidly evolving pattern compared with more conserved genes such as the COI and Cytb. The COI gene showed phylogenetic relationships among Pampas deer populations in the recent past, where many migrants per generation were exchanged between populations, and population genetic structure was almost non-existent. On the other hand, the characteristic D-loop hypervariable region reflects the phylogenetic relationships among contemporary populations of Pampas deer, where populations are genetically isolated without gene flow.

The results of the different mitochondrial markers showed that the D-loop has unique haplotypes within populations. The haplotype diversity value for COI gene was 0.35, the average nucleotide diversity by site was π = 0.00391 (s. d. = 0.0000003). The populations that exhibited the highest values were Paraná, Rocha, San Luis, and Salto.

We found 18 haplotypes of 63 samples in the Cytb gene sequence fragment of 417 bp in Pampas deer populations, and a diversity index of 0.28 (Table 6). These three mitochondrial markers provide resolution at different scales and allow us to elucidate different scenarios and the spatial connections of haplotypes to make inferences about demographic and evolutionary processes of this species (Figure 2, Table 7).

Figure 2 Minimum Spanning haplotype networks of mitochondrial regions showing the genetic relationships between Pampas deer individuals: A) Cytb b) COI and C) D loop. Mutational steps between haplotypes are shown as marks across connection lines. Circle sizes are proportional to haplotype frequencies and colours represent the respective sample localities. 

On the other hand, we found lower evolutionary rates for the COI gene than Cytb. From the nucleotide differences matrix between pairs of COI sequences we determined 18 haplotypes. The minimum spanning network constructed has the H1 and H5 haplotypes with an ancestral position and has a wide distribution (Figure 2). This is consistent in terms of geographic location of the different populations. The number of migrants per generation that is exchanged between these populations is less than one migrant per generation, in the case of Paraná and Santa Fe was the lowest value of migrants, Nm = 0.286.

The AMOVA results for the COI region analysis showed a lack of structure amongst populations. To explain this, we put forward the hypothesis that this lack of structure detected among populations is that millions of years ago these populations were connected. The mitochondrial control region (D-loop) analysis shows a later stage to this connection between haplotypes of the populations occurring after isolation by habitat fragmentation. The existence of high genetic variation featuring Pampas deer today indicates that the decrease in population size was recent.

Discussion

Past Genetic Variation. The direct ancestor of the Pampas deer first appeared in the Pampean Formation during the Pleistocene and may be associated with a glacial event approximately 2.5 million years ago at the boundary between the Gauss and Matuyama chrons (^{Bonadonna and Alberdi 1987}; ^{Marshall et al. 1983}). For analyzing the past genetic variation, we used two mitochondrial markers the Cytb and COI genes.

The divergence time estimated between the two species, Pampas deer and gray brocket deer (Subulo gouazoubira) was approximately 4.77 million years. The evolutionary rate found for the Cytb gene is 1.07 x10 ¯⁸ than it found for the COI gen 9.6 x10 ¯⁹ nucleotide substitution events per site per year would be happening (See Figure 2). The sequences for the D loop regions analyzed with 423 base pairs fragments showed an evolutionary rate of 1.4 x10 ¯⁸. The obtained divergence values for Pampas deer using the two mitochondrial markers, show that, as expected, the D-loop evolves faster than the COI gene.

The inference through the calculated evolutionary rate is that the COI gene could reflect another stage of the Pampas deer population history being more appropriate to analyze phylogenetics relationships. Cytb and COI genes are not able to resolve recent demographic events linked to habitat fragmentation (^{Tobe et al. 2010}). Estimation of migrants per generation showed a high gene flow among the populations studied, and because of this, we found less genetic population structure according to the COI gene.

On the other hand, nucleotide diversity per site for Pampas deer through the COI gene is 0.391 %. The Paraná population had the highest nucleotide diversity value (0.87 %) while the lowest one was evidenced in the Emas population (0.16 %). The nucleotide diversity among populations of Pampas deer, through analysis of the control region obtained by ^{González et al. (1998}), ranged between 1.1 % to 2.5 %, and Argentina's population had the lowest value. Their slower evolutionary rate than D-loop showed that existing populations in the past maintained a closer relationship, showing that the structure is a recent phenomenon and could be a consequence of genetic populations' isolation. A similar finding was observed in another endangered species like the Asian elephant (Elephas maximus), whose populations are rapidly declining and D-loop marker was suitable for analyzing genetic variation and to infer other processes such as introgression/hybridization (^{Srikulnath et al. 2023}).

Our results confirm the previous findings by ^{González et al. (1998}) that the control region of the Pampas deer is one the most polymorphic of any mammal. The large number of haplotypes and high level of nucleotide diversity in the Pampas deer suggest that this species was more abundant and widespread in the recent past. Since variability was lost rapidly from these populations, sizes have remained small for long periods of time (^{Ballou 1994}). However, currently Pampas deer are endangered in Argentina, south of Brazil and Uruguay, with fewer than 2,500 mature individuals. The levels of genetic diversity in populations from these locations suggest that historic population sizes were several orders of magnitude larger, and that recently populations have decreased dramatically, thus providing a strong mandate for restoration and augmentation. This population decline was due to habitat loss and unregulated hunting beginning in the last century and, most recently, to control efforts by ranchers who believe that deer compete with livestock. Pampas deer numbers might increase if protected from poaching in areas where natural habitats remain and if some grazing land, as a buffer, could be designated for dual use by deer and livestock (^{Castro et al. 2021}).

^{González et al (1998}) had estimated the historic population size based on the relation Ѳ = 2Nµ where N is the effective number of females and µ is the mutation rate per site. Using coalescent likelihood methods incorporated in the COALESCE program by ^{Kuhner et al. (1995}), the parameter Ѳ can be calculated from a population sample of DNA sequences. Our estimate of 2Nµ is 0.173 and assuming a mutation rate of 2.5 x 10^-8 per nucleotide site per year for the control region (based on sequence divergence between the Pampas and brocket deer, as above), the effective number of breeding females would be approximately 3,460,000. The total census size of females is likely to be at least double this value (e. g.^{Nunney and Elam 1994}). Therefore, both comparative and theoretical estimates indicate a substantial reduction in population size has occurred since the total present-day number of deer is between 64,000 and 80,000 individuals (^{González et al. 2010}).

Figure 3 Evolutionary relationships of populations. The optimal tree with the sum of branch length = 1.43701786 is shown. The tree is drawn to scale, with branch lengths in the same units as those of the evolutionary distances used to infer the phylogenetic tree. The evolutionary distances were estimated with Kimura 2 parameter and the analyses were conducted in MEGA11. 

Conservation Implications. Genetic units for conservation have been based on criteria such as reciprocal monophyly (evolutionary significant units) or differences in genotype frequency (management units; ^{Moritz 1995}). The numerous reticulations in the minimum spanning network (Figure 2) show that none of the neighboring Pampas deer populations are reciprocally monophyletic and indicate the occurrence of past episodes of migration.

The highest level of genetic differentiation was observed in Uruguayan populations, especially the Rocha population. This differentiation is explained, with the UPGMA algorithm, establishing a divergence time of 2 million years and considering mean genetic distance and D-loop mutation rate (Figure 3). Probably this separation is linked to events occurring during the Pleistocene.

However, except for the Argentinean populations from Samborombón Bay and San Luis, belonging to O. b. celer, and North Argentinean and Pantanal belonging to O. b. leucogaster, all the other populations are significantly or marginally differentiated, thus they might be classified as management units experiencing low to modest rates of gene flow. A pronounced sequence divergence exists between Brazilian populations from Emas National Park and Pantanal da Nhecolandia and North Argentinian populations (Figure. 3), corresponding to the different subspecific designation recognized by ^{Cabrera (1943}). These populations may have been historically isolated in different habitats. In fact, the population in Emas National Park is located in the cerrado of central Brazil, 650 to 1,000 m elevation, which has a distinct dry season, whereas the population in El Pantanal is found in wetlands below 100 m. Although based on limited evidence, these physiological differences may indicate differences in the timing of the reproductive cycle and hence, if of genetic origin, may be an important reason why the populations should not be interbred or used as a source for cross-translocation. Supporting differentiation between these populations are discrepancies in their physiology. In the Emas National Park population, antlers are shed in April, whereas in the Pantanal this occurs in June and July (^{González et al. 2010}). The individuals belonging to Paraná showed differentiated genetic distance and low number of migrants less than 1. This population is critically endangered, making it urgent to reduce activities such as poaching that may be severely affecting its conservation, and to promote conservation management actions and design mitigation measures to assure long-term survival.

Our results suggest that Pampas deer have the potential to exist over a much greater area and historical data demonstrate a much wider distribution for the species. Therefore, if the goal of conservation is to maintain long term population stability and preserve genetic variation, conservation efforts should focus on the restoration of deer habitats over a wide geographic area. Finally, we conclude that the genetic dynamic shown by Pampas deer allows us to identify the D-loop as the marker of choice for defining management units for conservation of this species.