Phylogeography of the Central american red brocket deer, Mazama temama (Artiodactyla, Cervidae) in southeastern Mexico

Anthropogenic threats have increasingly isolated the populations of Mazama temama (Erxleben, 1777) and limited the gene flow in this species. Knowledge of the phylogeographic structure of this species is therefore essential for its conservation. Thus, in this study, we describe the phylogeographic structure of two M. temama populations of Veracruz and Oaxaca, Mexico. We sequenced the D-Loop region of the mitochondrial DNA of 16 individuals, in order to estimate the diversity and genetic differentiation (FST), Tajima’s D index, "Mismatch distribution" test; a phylogram and a haplotype network was constructed and we performed multidimensional scaling analysis to test the hypothesis of association between geographic distance and genetic diversity. The haplotypic and nucleotide diversity was high, indicating divergent populations (FST = 0.223), while the Tajima’s D index (-1,03300; P > 0.10) determined disequilibrium in the D-Loop region, derived from a population expansion that was evidenced in the "Mismatch distribution" test and confirmed with the haplotype network in the form of a star. Four lineages were identified in the phylogram (Veracruz n = 3, Oaxaca n = 1), evidencing geographic Neotropical Biology and Conservation 16(2): 369–382 (2021) doi: 10.3897/neotropical.16.e58110 Copyright Ricardo Serna-Lagunes et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Ricardo Serna-Lagunes et al. 370 and reproductive isolation between the two populations. This was confirmed by the multidimensional scaling analysis, which evidenced recent evolutionary divergence between the populations analyzed, which are considered evolutionary units of conservation.


Introduction
The Central American Red Brocket, Mazama temama (Erxleben, 1777) is a deer species that has been targeted by illegal hunters for local commercial purposes or as a source of protein (Leopold 1959;Gallina and Mandujano 2009), as well as for its habitat (Gallina and Mandujano 2009;Bello-Gutiérrez et al. 2010). Even so, this deer species is considered Data Deficient by IUCN (Bello et al. 2016). This category is due to uncertainties about the threats that the species has suffered (Bello et al. 2016). Although threats such as those reported above are known, the degree of danger posed to the species is not clear (Bello et al. 2016). Thus, basic information about the species is essential to support conservation strategies.
Since certain populations of M. temama are found isolated in remnants of vegetation, it is important to study the phylogeography of this species in order to understand the micro-evolutionary processes these populations have undergone historically.
Phylogeography is a multidiscipline that has been used as a management tool for the conservation of species, and involves the study of micro and macroevolutionary processes that govern the geographic structure of genealogical lineages of species (Avise et al. 1987;Arbogast and Kenagy 2001;Vázquez-Domínguez 2007) and features two important components: demographic structure and genetic structure (Avise 1992). An approach in the phylogeographic context can thus help to clarify changes in genealogical lineages based on their geographic distribution (Avise and Hamrick 1996;Avise 2000).
On the American continent, particularly in Mexico, phylogeography has been used to study populations of wild fauna, especially endemic mammals, or those sought by hunters, or those included in some category of risk, leading to important implications for their conservation (Vázquez-Domínguez 2005;Vega et al. 2007;Vázquez-Domínguez et al. 2009;Escobedo-Morales et al. 2016). There is evidence that phylogeographic studies have had a considerable impact on species management (Vázquez-Domínguez and Vega 2002), especially in the ungulates, describing cervid lineages in order to define particular strategies of management (Serna-Lagunes et al. 2021), identifying bottlenecks and low gene flow in species of interest to hunters (Castillo-Rodríguez et al. 2020) and evidence of genetic divergence among lineages of the same species (Fuentes-Hurtado et al. 2011). In this sense, in this study we describe the phylogeographic structure of two populations of M. temama.

Collection of biological material
We were able to generate molecular data from 10 specimens of M. temama, among two neighboring populations in the states of Veracruz and Oaxaca, in Mexico (Table 1, Fig. 1). Eight live specimens were confiscated by the environmental authority (Procuraduría Federal de Protección al Ambiente-PROFEPA), which authorized the collection of blood and skin samples. The specimens were seized in five localities of the Sierra de Zongolica, Zongolica, Tehuipango, Tequila and Tezonapa in Las Montañas region, in the center of Veracruz, where the main threats to populations of this species are the fragmentation of their habitat, the isolation of their populations in secondary vegetation, as well as the hunting and illegal possession of these deer as pets (Salazar-Ortiz et al. 2020). The specimens are currently being kept in the wildlife management unit (UMA, by its Spanish acronym) "El Pochote", in Ixtaczoquitlán, Veracruz. The two other skin samples come from two localities in the state of Oaxaca and were provided by the Coleccion Nacional de Mastozoología of the Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional (CIIDIR), Oaxaca unit, of the Instituto Politécnico Nacional, with registry of origin in Oaxaca (Table 1).

Extraction of DNA, PCR and sequencing
The blood samples were dehydrated at ambient temperature and the skin samples fragmented in order to extract the DNA in accordance with instructions from the DNeasy Blood and Tissue kit of QIAGEN. DNA extraction was verified by applying electrophoresis in a horizontal 1% agarose gel with ethidium bromide at 70V for 40 min in 1X buffer TAE. The bands of the extracted DNA were visualized in a transilluminator and documented with the program GelCapture.
From each extracted DNA sample, 2.5 µL were taken to perform the Polymerase Chain Reaction (PCR) at the final phase or end-point. From this DNA, the Con-  trol or D-Loop region of the DNA mt was analyzed using the universal primers: DL-H16498 (5' CCTGAACTAGGAACCAGATG 3') and Thr-L15910 (5' GAATTCC-CCGGTCTTGTA 3') (Vilà et al. 1999). The reactions were performed in a final volume of 25 µl and contained the following: 2.5 µl of DNA (< 250 ng), 1.25 µl of the Forward and 1.25 µl of the Reverse primers, both at a final concentration of 10 µM, 12.5 µl of 2X Promega PCR Master Mix (25mM Tris-HCl pH 9, 25mM NaCl, 2.5mM MgCl2, 100 µM each of dATP, dGTP, dCTP, dTTP, 0.5 U of Taq DNA polymerase, 0.05 mg/ml BSA) and 7.5 µl of nuclease-free water. The PCR was conducted in a thermocycler, using the following program: initial denaturalization at 95 °C for 1 min, followed by 35 cycles of amplification: 94 °C for 1 min, annealing at 51 °C for 2 min, extension at 72 °C for 1 min 30 sec and a final extension step at 72 °C for 7 min. This procedure was conducted in the Laboratorio Cordobés de Diagnóstico Pecuario, S.C., located in Córdoba, Veracruz, Mexico. The PCR Products were verified using Green DNA Dye. Verification of the PCR products was corroborated using the same DNA extraction method. The PCR products were sent to the company MacroGen (Korea) for purification and sequencing. Macrogen employs the Sanger technique, with capillary technology and an Applied Biosystems 3730XLs sequencer in order to perform simple sequencing of nucleotides (www.macrogen.com).

Statistical analysis
A database was generated with the obtained sequences, grouping them according to population. The sequences were aligned with the algorithm Clustal W, installed in the program MEGA v. 6 (Tamura et al. 2013). For each population, the nucleotide composition (proportion of each nitrogenous base) of the sequences was estimated. The number of polymorphic sites (S), number of identified haplotypes (H i ), number of haplotypes shared between populations (H s ), haplotypic diversity (H d ) and nucleotide diversity (π) were calculated. The index of genetic differentiation (F ST ) was also calculated as a measure with which to distinguish between the populations under study. The Tajima's D neutrality test (Tajima 1989) was applied as an indicator of historic changes in population expansion. The "Mismatch distribution" test was performed to measure the effect on population size of the historic demographic changes that molded the genetic diversity observed. This analysis is represented in a graph that compares the observed and expected nucleotide frequencies (Rogers et al. 1996). The fit of this test was verified with the values of R 2 , initial theta (θ) and tau (τ) (Ramos-Onsins and Rozas 2002). These analyses were performed with the software DnaSP v.5 (Librado and Rozas 2009).
The model of nucleotide substitution that best fitted the sequences of M. temama was identified (Nei and Kumar 2000) using the software MEGA v.6 (Tamura et al. 2013). The model of Kimura 2 (K2) was that which best described the pattern of nucleotide substitution by presenting a lower number of parameters (18) under the coefficient of the Bayesian Information Criterion (BIC = 3394.84), with a mean value of the Akaike Information Criterion (AICc = 3301.260) and the highest value of the parameter of Maximum Likelihood (lnL = -1632.379) (Nei and Kumar 2000;Tamura et al. 2013). The most parsimonious phylogram was constructed with Maximum Parsimony (MP), based on the Branch and Bound algorithm (Nei and Kumar 2000). The fit of the phylogram branches was evaluated with 1000 Bootstrap replicates (Felsenstein 1985) using the software MEGA v.6 (Tamura et al. 2013), following the protocols of Hall (2004 and2013). A haplotype network was constructed with Median-Joining Network and MP. This determines the mutational steps between haplotypes and was generated with the software Network v.5 (Fluxus Technology Ltd 1999. Finally, nucleotide distances were estimated between pairs of sequences with the software DNAsp (Librado and Rozas 2009), in order to determine the degree of current genetic similarity or divergence, where short genetic distances indicate close relationships while long distances suggest distant relationships. Likewise, geographic distances (km: between collection localities) were obtained between pairs of sequences using ArcMap v.10.3. With these distances (nucleotide and geographic, see Table 4), a multidimensional scaling analysis of two dimensions (MDS; Borg and Groenen 2003) was conducted to determine the similarity (or dissimilarity) among the sequences. This analysis is a test of interdependence that calculates a matrix of paired Euclidean distances, producing a graph with a geometric space of two dimensions (Lessa 1990), on which are represented the proximities that exist between the nucleotide and geographic distances of each sequence; the fit of the graph was evaluated with the values of "Stress" (>0.2 good or >0.7 very good) and R 2 (>0.3 good fit or >0.7 very good fit). In addition, a principal component analysis (PCA) was applied to the matrix of nucleotide and geographic distances with the Varimax rotation method (Yakubu et al. 2011). Both analyses were conducted with the software SPSS 15.0 (SPSS 2004).

Results
The analysis of the nucleotide composition of the obtained sequences indicated a higher average frequency of Cytosine bases in the M. temama population of Veracruz, while the Oaxaca population presented a greater proportion of Adenine bases (Table 2).  The two populations of M. temama showed differences in genetic diversity. In the population of Oaxaca, there was a lower number of polymorphic sites compared to the population of Veracruz, which explains why the number of haplotypes was different. That is, each sample analyzed was found to be an independent haplotype, and haplotypes were not shared between the populations. Haplotypic and nucleotide diversity was high in both populations, indicating that the sampled individuals of each population were heterogeneous in terms of the D-loop region (Table 3). This is consistent with the results of the index of genetic differentiation (F ST ) = 0.223, which indicated genetic differences between the populations of Oaxaca and Veracruz, as a result of not sharing haplotypes (Hi = 0).   The Tajima's D index was negative (-1,03300; P > 0.10), indicating the presence of deletions or mutations, such that the populations were not found in neutral equilibrium. The "Mismatch distribution" test fitted to the model of a population of constant size (Oaxaca: R 2 = 0.500, initial theta (θ)= 14.750, tau (τ)= 48; Veracruz: R 2 = 0.217, initial theta (θ)= 14.75, tau (τ)= 82.571), evidencing a possible population expansion by the variation shown in the graphs of this test (Fig. 2).
The haplotype network presented a structure in the form of a star (Fig. 4), with 351 total non-segregating sites and 233 mutations. This is commonly observed in populations that have undergone a process of demographic expansion and was consistent with the results of the Tajima's D index and the "Mismatch distribution" test.  According to the MDS (Fig. 5), the M. temama individuals most related according to nucleotide distance were those found in southern Veracruz and Oaxaca (upper right extreme of the dimensional plane), while the least related were those of central and southern Veracruz (extreme lower left of the dimensional plane), which explains the genetic relationship based on geographic distance; the stress value was good (0.14), with a value of R 2 = 0.87, indicating a significant relationship between the pairs of sequences of the matrix of nucleotide and geographic distances (Fig. 5). Finally, the PCA determined that the greatest proportion of the variance was accommodated in the first two principal components (>80%). Principal component 1 explained 63.8 of the variances (Table 5), while principal component 2 explained 16.3% and together these components explained 80.3% of the total accumulated variation (Fig. 6).

Discussion
In this study we performed our first phylogeographic analysis on the species M. temama. Our results show that the populations analyzed underwent a genetic bottleneck and rapid population expansion.
Although we have studied a small sample of the population, these populations studied here do not share haplotypes due to the geographical distance between the two populations. On the other hand, it is possible that the reproductive behavior  of this cervid means that the frequency of mtDNA transfer decreases, due to the philopatry of the females (Purdue et al. 2000;Purdue et al. 2006). It is only the males that migrate great distances to mate and thus the lineages of the females are presented in greater proportion with proximity to their area of birth (Nelson and Mech 1992). The lack of shared haplotypes between the two studied M. temama populations is evidence of the geographic and reproductive isolation. The results of the "Mismatch distribution" test show an evident multimodal form, reflecting processes of sudden population expansion associated with a stochastic process of lineage extinction through genetic drift, recorded in the genetic structure of the populations (Rogers and Harpending 1992). This may reflect the threats that the species suffers in Mexico; M. temama is of hunting interest and subject to constant pressures of exploitation, habitat fragmentation, poaching, auto-consumption by local communities and predation by carnivores (Gallina 2005). This has led to a reduction in the effective population size where only certain lineages pass to the next generation, while others are extinguished or are present in lower frequency in the population. These results are to be expected considering that M. temama is a species that suffers from poaching, since a similar phenomenon has been reported in the cervid Pantholops hodgsonii in China, which is at risk of extinction due to the loss of genetic diversity (Xiang-Dong et al. 2005).
The populations analyzed have different evolutionary histories, as shown by the string haplotypic and nucleotide relationship between the populations, evidencing differentiated lineages between the two populations (Grant and Bowen 1998). It is concerning that the populations studied present low values of gene flow, as occurs in some other species of Neotropical cervids. This is associated with the fact that M. temama inhabits fragmented ecosystems with no biological corridors to allow migration of individuals for the purposes of mating. Phylogeographic studies must be conducted in more populations of M. temama in order to evaluate their status of genetic diversity and to generate a program of conservation of their diversity, since their ecological niche is restricted to vulnerable ecosystems due to the reduction and fragmentation of their habitat (Serna-Lagunes et al. 2014).