Identifying Priority Giant Anteater (Myrmecophaga tridactyla) Populations for Conservation in São Paulo State, Brazil

Abstract Habitat loss is the main threat to biodiversity conservation worldwide. Some species may be particularly susceptible to the effects of fragmentation and the isolation of populations. The impacts of human activity on wild animal populations may be understood through relationships between individual genetic data and spatial landscape variables, particularly when considering local population dynamics influenced by fragmented habitats. Thus, the objective of this study was to analyze the population structure and genetic diversity of the giant anteater (Myrmecophaga tridactyla) using an individual sampling scheme (ISS) on a regional geographic scale. Data were collected from 41 specimens from twenty different locations in São Paulo State, Brazil, and six polymorphic microsatellite loci were genotyped. Our results indicate that barriers to gene flow exist and have segregated individuals of the farther away areas into two spatially structured clusters. The populations were also found to have high genetic diversity. The experimental sampling approach used herein enabled an analysis of the population dynamics of the giant anteater on a regional scale, as well as the identification of priority populations for genetic resource conservation for this species. The results reflect the need for adequate management plans. The efficacy of the sampling scheme may vary based on the study model used, but we argue that the use of an ISS combined with suitable molecular markers and statistical methods may serve as an important tool for initial analyses of threatened or vulnerable species, particularly in anthropized regions where populations are small or hard to characterize.


| INTRODUC TI ON
Decreases in and losses of habitat connectivity lead to increases in the existing decline in wild animal populations, thus raising the risk of their extinction (Butchart et al., 2010). The transformation of a continuous habitat into "small pieces of habitat" may affect biodiversity at the population level, with changes in species number, distribution, reproduction, and survival (Carvalho et al., 2009;Fahrig, 2003;Rocha et al., 2018;Wolff et al., 1997).
Medium-and large-sized land mammals may be particularly vulnerable to habitat fragmentation and loss, particularly for species with large home ranges or long life cycles (Keinath et al., 2017;Morris et al., 2008). A large-sized animal averaging 31 kg, the giant anteater (Myrmecophaga tridactyla; Linnaeus, 1758) (Figure 1), is considered a vulnerable species in Brazil largely because of its habitat loss in the Cerrado and other Brazilian biomes (Miranda et al., 2015).
There is a tendency for the remaining populations of M. tridactyla to become progressively more isolated, which, in turn, further increases their vulnerability and leads to various cases of local extinction. Some of the species' biological characteristics may lead to low rates of population growth: a low metabolic rate and low thermoregulatory capacity, limited reproductive potential, prolonged parental care, a long gestation period, and solitary habits (Clozato et al., 2017). The species is a frequent victim of roadkills, as well as of fires caused by sugar cane crop burning and prolonged periods of drought. These factors threaten the persistence of giant anteater populations in anthropized areas (Diniz & Brito, 2013. The extinction of small populations is expected in many fragments over time; the study of local population dynamics influenced by the existence of highly fragmented habitat networks may provide clearer information that can be used to implement strategies to help protect and manage isolated populations (Canale et al., 2012;Weiss & Leese, 2016).
Conservation genetics focuses on the effects of contemporary genetic structuring on long-term survival of a species, and genetic analyses have been used to determine the prospective status of several endangered species (Wan et al., 2004). This study sought to characterize the giant anteater's population structure in fragments of Cerrado and Cerrado/Atlantic forest ecotones, to evaluate its genetic diversity, to detect any barriers to gene flow, and to determine the current degree of gene flow in anthropized areas in São Paulo State, Brazil.

| Study area
The animals sampled herein were largely distributed in the northwestern and central-southern region of São Paulo State, which is the most industrialized state in Brazil. This region is currently experiencing intense expansion and urban development, with native vegetation being highly affected by human activity (Durigan et al., 2007).
Before colonization, this region was covered by a mosaic of Cerrado and seasonal semi-deciduous forest; however, the forests have since been almost entirely depleted (Myers et al., 2000). The vegetation cover of São Paulo State has been drastically reduced, and the remnants correspond to small areas surrounded by pasture, sugarcane, soybean, reforestation, perennial crops, and urban zones (Kronka et al., 2005).
The study area included some animals that had been obtained, most frequently after car accidents, in cities that border the Noroeste Paulista Conservation Area (EENP) located in the cities of São José do Rio Preto and Mirassol. Other animals were obtained in search-and-capture campaigns in and around the Santa Bárbara Conservation Area (EESB) located in the city of Águas de Santa Bárbara. The use of animals from these locations permitted inferences about the status of individuals and populations in habitats in and around these conservation areas.
The EENP conservation area is 500 ha in size when combined with the Noroeste Paulista State Forest and is considered to be seriously threatened by urbanization and real estate speculation.
It is surrounded by pastures, farms, and urban neighborhoods (IF -Instituto Florestal, 2014). The EESB conservation area is 2,712 ha in size, and most of the area is occupied by native vegetation. There is little evidence of human disturbance, and the area is surrounded by pastures and reforestation sites (IF -Instituto Florestal, 2011).

| Sampling
The current study included a local analysis involving an ISS-individual sampling scheme (Prunier et al., 2013). The ISS used allowed for giant anteaters in the northwestern and central-southern regions of São Paulo State to be analyzed.
The samples included in this study (n = 41) were obtained from different sources. Some were collected by veterinarians and park rangers after the collection of animals that had been run over or had been victims of crop burning (n = 35). In other cases, blood samples F I G U R E 1 The giant anteater, Myrmecophaga tridactyla (Pilosa: Myrmecophagidae). Photo credit: Jason Woolgar were obtained from animals captured alive in field studies for other purposes and were provided to our laboratory (n = 6).
Samples were obtained from a total of twenty different sites.
The geographic distances between the sampling points ranged from 10 km to 340 km ( Figure 2). The maps were generated in DIVA-GIS, version 7.5 (Hijmans et al., 2012).
The analyses of interindividual information used in this study may be applied on practically any spatial scale, including in analyses of genetic data that use one or several individuals from different sampling sites (Guillot et al., 2009;Mank & Avise, 2004;Miller, 2005;Safner et al., 2011). The number of individuals required to evaluate genetic structure depends on the extent of population differentiation, as well as on the degree of polymorphism in the markers; these factors mean that a small number of sample animals may be appropriate in certain circumstances

| Amplification conditions
Six pairs of microsatellite-specific oligonucleotides were used in the molecular analyses, as detailed in Table 1 (Garcia et al., 2005).
Genomic DNA was extracted using the QIAmp DNA kit (Qiagen) following the manufacturer's instructions. It was quantified using an Epoch™ Spectrophotometer System (BioTek) to measure DNA quality and concentration.  Amplification via PCR was performed with the addition of the M13 tail, which allowed for the use of the tailed primer method (Schuelke, 2000). The amplifications were performed using touchdown PCR in a total volume of 17.2 μl, which contained 1 μl of DNA (~50 ηg), 0.2 μl of the forward primer (5 µM), 0.4 µl of the M13 primer dyed with fluorophore (5 µ), 0.6 µl of the reverse primer (5 µM), 5 µl of the Gotaq ® Colorless Master Mix (Promega; 2X), and 10 µl of ultrapure water (Don et al., 1991;Korbie & Mattick, 2008).
The amplification conditions included denaturation at 95°C for 2 min, 12 cycles at 95°C for 1 min, between 64°C and 52°C for 40 s (−1°C per cycle), 72°C for 30 s, plus 25 additional cycles (1 min at 95°C, 40 s at 52°C, and 30 s at 72°C), and the final extension for 5 min at 72°C. After the amplifications were confirmed in 3% agarose gel, the fragments were identified in an ABI 3500 Genetic Analyzer (Applied Biosystems). The DS-33 GeneScan™ Installation Standard (Applied Biosystems) was used for the capillary electrophoresis run, and the reactions were organized into two groups based on the expected fragment sizes.
The results and peak sizes were analyzed using GeneMapper ® , version 4.1 (Applied Biosystems, 2012). The correct identification of the alleles is crucial for a reliable interpretation of the microsatellite data (Arif et al., 2010). No artifacts or stutter peaks negatively affected the differentiation between homozygote and heterozygote, and the alleles exhibited well-resolved peaks. The fragments were larger in size relative to the size of the original markers due to the addition of the use of the M13-tailed primer, but the sizes corresponded to the expected sizes based on the information available on their development (Garcia et al., 2005).

| Statistical analyses
GENEPOP Software, version 4.2 (Raymond & Rousset, 1995), was used to determine the presence of null alleles based on the maximum likelihood method, (Dempster et al., 1977). No genotyping failure was detected and the original dataset was used for the remaining analyses.
To characterize the population structure, individual-based Bayesian analyses were employed using the programs STRUCTURE v.2.3.4 (Pritchard et al., 2000) and Geneland 4.0.5 (Guillot et al., 2005) to determine the most likely number of clusters (K).
In STRUCTURE, we ran 25 replicate runs for each potential number of genetic clusters with 100,000 burn-in steps, 20,000 MCMC repetitions, and 25 iterations per run. Values of K = 1 to K = 5 were tested using the model with admixture, correlated allele frequencies, and the ancestry prior option "separate α for each population" (Falush et al., 2003). The results were then used in STRUCTURE HARVESTER (Earl & Vonholdt, 2012) to select the K value associated with the highest mean posterior probability of the data (Evanno et al., 2005) with postprocessing in the CLUMPAK program (Kopelman et al., 2015).
Unbalanced sampling has a large impact on STRUCTURE analysis. It reduces the quality of individual assignments to populations and the accuracy of the estimated number of populations. However, these adverse effects of unbalanced sampling on Structure analysis can be largely overcome by simply switching to the alternative ancestry prior, at least when the number of populations is not large and STRUCTURE can yield highly accurate inferences of individual ancestries (Wang, 2017).
Population structure was also analyzed with the software GENELAND 4.0.5 (Guillot et al., 2005) and R version 3.4.2 (R Core Team, 2019). GENELAND differs from STRUCTURE in that geographical information can be incorporated to produce more accurate inferences of population structure based on the spatial distribution of individuals (Levy et al., 2012).
The spatial structure and posterior probability of belonging to a given cluster were determined using the correlated and null alleles models, for runs of K = 1 to K = 5, and 10 multiple runs were done to check convergence. Each run consisted of 100,000 MCMC iterations with a thinning of 100 and a burn-in of 200. The correlated alleles model may be more powerful at detecting subtle differentiation. Thus, the correlated alleles model could be more informative than the default model of uncorrelated frequencies and hence get TA B L E 1 Microsatellite-specific oligonucleotides (Garcia et al., 2005) more accurate results, in particular in presence of low differentiation and when K is unknown (Guillot, 2008). The uncertainty of coordinates was set to 10, as this is a computational step that allows for the inclusion of information on the mobility of the species in question (Guillot et al., 2011), according to the characteristics of the species (Medri et al., 2006) and the landscape (Durigan et al., 2007).
Many types of statistical analyses have been applied to understand how the landscape structures genetic variation in natural populations (Guillot et al., 2005;Manel et al., 2003). Bayesian spatial cluster analysis and edge detection methods are frequently useful for describing patterns and revealing microevolutionary processes between individuals and within populations (Guillot et al., 2009;Safner et al., 2011). Thus, in addition to the Bayesian spatial cluster analysis implemented in GENELAND, we performed the edge detection method in the Alleles in Space (AIS) Software (Miller, 2005).
Monmonier's algorithm (Monmonier, 1973) was implemented in the individual-based program AIS to detect the presence of barriers to gene flow and evidence of isolation by barrier (IBB). Monmonier's algorithm is supervised, and the number of biogeographical barriers to be calculated must be specified beforehand (Miller, 2005).
We computed the two first barriers that were contiguous to each other, once the forming boundary has closed on itself by forming a loop around a population (Manni et al., 2004). The barriers were inferred using a Delaunay triangulation-based connectivity network and residual genetic distances (Manni et al., 2004;Monmonier, 1973).
This combination of spatial statistical analyses involving Bayesian clustering and barrier detection methods involving Monmonier's algorithm may be a powerful tool for maximizing accuracy and minimizing type-1 errors in barrier detection (Blair et al., 2012).
Aggregation indices were also calculated in AIS. These indices are usually used in ecological studies to quantify spatial patterns and can also be useful for testing genetic diversity patterns over a given landscape. All of the analyses in AIS were calculated with 100,000 permutations/iterations (Miller, 2005).
To determine whether there was evidence of isolation by distance (IBD), the test provided by Mantel (Mantel, 1967;Sokal, 1979) was performed in the VEGAN package of the R software, Cluster analysis was also performed in PAST (Hammer et al., 2001) using the same subdivision and Bray-Curtis index (Bowlby et al., 2016;Shirk et al., 2017), according to the Neighbor-Joining (NJ) method (Saitou and Nei, 1987). The nodes were supported by bootstrap analysis with 10,000 replicates.
Parameters for the evaluation of the microsatellite loci were studied. Estimates of genetic diversity and microsatellite variation in the population of São Paulo State (all sampled individuals) and within the subpopulations of the conservation areas, from which the average allelic richness and gene diversity (HE) were calculated.
Finally, genetic divergence between populations was estimated by calculating the unbiased estimate of F ST , and the gene flow analysis was computed using the number of migrants per generation based on the method involving private alleles (Barton & Slatkin, 1986;Weir & Cockerham, 1984) using the GENEPOP program, version 4.2 (Raymond & Rousset, 1995).

| Bayesian clustering analyses
The results of all Bayesian clustering methods yielded the same optimal number of clusters. We identified K = 2 in STRUCTURE as the most probable number of clusters after calculating ΔK (Figures 3 and   4). These clusters corresponded to the individuals of the northwest-  This result indicates that the samples presented a clumped spatial distribution (RjAVE < 1) and supports the structure found in the Bayesian spatial analysis.

| Isolation by barrier and Isolation by distance
In other words, the analyses as a whole indicate that spatial structuring does exist, but that distance is not an isolated factor that led to this population structure pattern.

| Principal component analysis
The The EESB population showed less variation, but a similar average to the NP population (ellipses centers almost together), and the EENP population stands out mainly to PC2. Each ellipse comprises (in general) 95% of the population, and this explains the elements that are external to the ellipse (outliers).

| Clustering methods
The dendrogram based on the NJ algorithm revealed two discrete

clusters.
Results presented high bootstrap support values (Figure 8), a finding that corroborates the results of the individual-based clustering analyses.

| Genetic diversity and population differentiation
The Descriptive analyses indicate that all of the loci were highly polymorphic. After the application of Bonferroni's correction, two loci exhibited deviations from Hardy-Weinberg equilibrium, and the tests for linkage disequilibrium were found not to be statistically significant (p > .05). Analyses were also performed between individuals found close to the EENP conservation area (n = 10) and those found close to the EESB conservation area (n = 6) to investigate the genetic diversity and differentiation between the population of the conservation areas. The results reflected greater genetic differentiation between these two groups (F ST = 0.1739) and relatively low/moderate F IT (0.1348) and F IS (0.0912) values.

| Gene flow analysis
Gene flow analyses were carried out for testing levels of gene flow among individuals belonging to different areas, and the results showed that the number of migrants per generation is lower between two areas that are farther away from each other. This reduction in the migration rate at regional levels may have led to population differentiation. Also, we were able to infer migration between the animals that were obtained in the portions of the landscape surrounding the conservation areas and the result indicated a lower level of migration between the populations of the EENP and EESB (N m = 0.615729).

| Genetic diversity of Myrmecophaga tridactyla
The use of microsatellite markers in this study aided in the identification of substantial intrapopulation variance, population subdivision (K = 2), and high rates of genetic diversity. The genetic structure of the giant anteater populations in this region is best explained by the presence of barriers to gene flow than by geographic distance. Few studies have been performed to evaluate population dynamics and the genetic variability of giant anteater populations (Clozato et al., 2017;Collevatti et al., 2007). Garcia et al. (2005) developed six species-specific microsatellite primers and evaluated allelic diversity of markers from fifteen individuals (eight animals that had been run over in the states of São Paulo and Mato Grosso, and seven that had been collected in the state of Goiás).
They obtained levels of genetic diversity similar to those found in our study (H o = 0.61, H e = 0.63, and 26 alleles total), considering the small sample size.
Using the same microsatellite markers, Collevatti et al. (2007) analyzed 27 giant anteaters from Emas National Park (which is also located across the states of São Paulo and Mato Grosso) and found low levels of genetic diversity (mean number of alleles = 3; tridactyla, or it may be a consequence of population bottlenecks due to the decades-long history of recurrent fires in the area (Collevatti et al., 2007;Silveira et al., 1999).
The high rates of genetic diversity found were unexpected: The populations are structured, exhibit genetic differentiation, and are located in a fragmented habitat, factors which can decrease genetic diversity (Frankham, 2005). However, this relationship is not always linear; thus, not all large populations will exhibit high genetic diversity, and not all bottlenecked or decimated populations will experience a reduction in genetic diversity (Torres-Florez et al., 2014).
Furthermore, changes to genetic diversity associated with habitat fragmentation may be subtle due to recent environmental disturbances (Dixo et al., 2009).
Though some species are highly sensitive to human activity and experience local extinction and/or population reduction when urban development begins, other species are able to persevere as small and relatively isolated populations (Ramalho et al., 2018). As a consequence, in a short period and without the occurrence of extinction, many populations are expected to exhibit fixation more rapidly, but to retain greater overall genetic diversity than a single large population. However, in the long term, when small populations become extinct, the remaining small populations will retain less genetic diversity than the original single large population (Frankham et al., 2002). Therefore, the genetic consequences of habitat fragmentation in a single large population are different from the genetic consequences for various small and isolated population fragments, even in large protected areas.
The genetic patterns observed in this study may be explained by environmental, evolutionary, and biological factors that have led to regional differences in the genetic frequency of isolated populations. Levels of genetic diversity may be maintained or lost through various types of selection, as well as by the effects of genetic drift in selectively neutral mutations in finite populations (Charlesworth et al., 1997;Crow & Kimura, 1970;Gillespie, 1991;Kimura, 1983).

| Barriers to gene flow
Barriers to gene flow were identified, and it can be inferred that the In general, the current results show that barriers have prevented gene flow at the regional level, but migrations may have occurred at shorter distances and may have increased intrapopulation variation at the local level contributing to the sustained genetic diversity observed herein (Barton, 2008). Following this logic, differences between the most geographically distant populations of São Paulo State would have increased over time and led to cluster differentiation, but with substantial genetic variation at short distances due to probable short-distance migrations.
Population subdivision may also lead to a geographic structure that affects allele frequencies in space and proportions of different genotypes in local populations (Garnier-Géré & Chikhi, 2013). Thus, gene flow does not depend on distance alone, but also on the nature of the landscape surrounding and between the populations; the separation of populations may have favored the retention of genetic variants in small and isolated populations (Templeton et al., 2001).

| Biological characteristics of the species
In addition to evolutionary and environmental factors, certain species' biological characteristics may influence the degree of genetic diversity loss in populations; for example, a long life span seems to contribute to sustained genetic variability in wild populations (Torres-Florez et al., 2014). The processes that lead to the loss of genetic diversity may be buffered by intrinsic biological characteristics, and long generation times may result in small populations that seem genetically diverse despite periods of substantial decline (Hailer et al., 2006).
The estimated lifespan of the giant anteater is more than 15 years in the wild and 20-30 years in captivity. The species' generation time is estimated to be approximately seven years (Gaudin et al., 2018;Knott et al., 2013;Medri et al., 2006;Miranda et al., 2014;Nowak, 2018). Studies on long-living mammals have identified high rates of genetic variability in many small and isolated populations; in some cases, even species that are almost extinct still exhibit high genetic diversity (Dinerstein & Mccracken, 1990;Swart et al., 1994;Taylor et al., 1994;Zavala-Páramo et al., 2017).
Species with shorter generation times and smaller body sizes typically tend to experience faster genetic responses to man-made barriers (Epps et al., 2005). Therefore, a loss in genetic diversity in a long-living species cannot be easily detected in a short period of time, and this factor should be considered when evaluating the species' risk of extinction. Thus, it is plausible that the long life span of the giant anteater combined with its demographic and evolutionary history delayed the loss of the genetic diversity in these populations.
Other hypotheses to explain the processes that have maintained genetic diversity in these populations should not be discarded, but regardless of the hypotheses and theories that explain population patterns, the analyses performed herein revealed that giant anteaters in the state of São Paulo should be a priority in conservation measures due to their high genetic diversity and variability.

| Conservation issues
Genetic criteria and allelic richness in particular are crucial tools for selecting candidate populations of wild species to be prioritized for conservation (Petit et al., 2008). Thus, despite this persistently high genetic variability, factors such as habitat fragmentation, habitat destruction, and demographic instability may lead wild populations to extinction before the genetic response to environmental impact becomes clear, particularly in long-living species (Lacy, 1997).
Conservation plans are therefore essential to increase the chances of giant anteater survival and to enable the conservation of genetic resources for subsequent species recovery (Théry, 2011).
Due to the use of convenience sampling in this study, the degree of genetic diversity observed may be the result of the use of various study sites with one or more individuals; but in any case, these individuals as a whole represent substantial genetic richness for a species living in an anthropized region. The findings reported herein indicate that management plans must be implemented at local and regional levels because, in the context of intense human activity, the giant anteater population is likely to be limited to a few individuals isolated in the remaining suitable habitats around the state (Bertassoni et al., 2019).
The frequency at which giant anteaters are run over on local highways shows that road ecology studies are necessary to determine local mitigating measures to consider when implementing fauna passage systems; the ultimate goal of this research and these policy measures should be to prevent the decline of animal populations, including that of the giant anteater (Teixeira et al., 2017).
Finally, conservation genetic studies on small geographic scales may provide important information for determining the status of threatened species populations. It is also crucial that interactions between scientists and policymakers improve so that the results of genetic and ecological research can be applied in the establishment of effective conservation strategies (Montgelard et al., 2014;Vernesi et al., 2008).

| FINAL CONS IDER ATIONS
This study has shown that population of M. tridactyla in the state of São Paulo, including the conservation areas, maintain high levels of genetic diversity. Though the sample size is relatively small, it seems to be the largest sample size in conservation genetics studies performed on the giant anteater at the regional level available in the literature.
The results emphasize the importance of fragmented areas for the maintenance of genetic diversity in anthropogenically modified landscapes and support the creation of adequate conservation plans that can work to prevent the local extinction of genetically important populations.
Environmental causes such as roadkills and the increasing agricultural borders can lead populations to local extinctions before a perceptible genetic response to contemporary habitat fragmentation. Thus, the populations in the state of São Paulo should take priority in the conservation of the genetic resources of the species.

ACK N OWLED G M ENTS
We would like to thank the people and institutions that contributed directly or indirectly to this study, especially the São Paulo Research

CO N FLI C T O F I NTE R E S T
None declared.