Genetic assignment of illegally trafficked neotropical primates and implications for reintroduction programs

The black and gold howler monkey (Alouatta caraya) is a neotropical primate threatened by habitat loss and capture for illegal trade in Argentina. Using multilocus microsatellite genotypes from 178 A. caraya individuals sampled from 15 localities in Argentina, we built a genotype reference database (GRDB). Bayesian assignment methods applied to the GRDB allowed us to correctly re-assign 73% of individuals to their true location of origin and 93.3% to their cluster of origin. We used the GRDB to assign 22 confiscated individuals (17 of which were reintroduced), and 3 corpses to both localities and clusters of origin. We assigned with a probability >70% the locality of origin of 14 individuals and the cluster of origin of 21. We found that most of the confiscated individuals were assigned to one cluster (F-Ch-C) and two localities included in the GRDB, suggesting that trafficked A. caraya primarily originated in this area. Our results reveal that only 4 of 17 reintroduced individuals were released in sites corresponding to their cluster of origin. Our findings illustrate the applicability of genotype databases for inferring hotspots of illegal capture and for guiding future reintroduction efforts, both of which are essential elements of species protection and recovery programs.

Similar to other countries, wildlife illegal hunting and trade are threats to Argentinian wildlife. Confiscated and surrendered animals from trafficking are transported to rehabilitation centres, and the return of these confiscated animals to the wild receives strong support from the public. Ideally, trafficked animals would be reintroduced into the population that they were extracted from or translocated to another suitable site within the species' original range. Although translocations are considered a good option for conservation and a solution to trafficked animals 1 , they might be detrimental for the animal and/or the environment if choices over where to return individuals to their natural habitats are not properly based on scientific evidence [2][3][4] . For example, when significant genetic structure exists within the species in question, translocations may inadvertently lead to admixture of distinct evolutionary lineages and act to homogenize existing diversity and biogeographic patterns instead of protect them 5 . However, the genetic consequences of translocations have seldom been studied 6 , and in Latin America, the development of biodiversity management and conservation plans as part of public policy have not yet taken advantage of newly developed genetic techniques to inform translocation policy and decisions 7 . Conservation genetics can help strengthen the links between scientists and decision makers and improve reintroduction and translocation policy.
Molecular genetic studies using individuals of known origin allow researchers to calculate levels of differentiation among populations and to assess population structure within a species in the wild [8][9][10] . Based on this, the populations studied may be grouped into clusters based on genetic similarity. Once such an assessment of differentiation between populations/localities or clusters is obtained, genetic assignment analysis of individuals of unknown provenience can be performed to identify the likely locality or cluster to which they belong -i.e., to

Results
Genotypes of 143 individuals from 10 localities (Locs) were available from previous studies 9, 28 . To these, we added newly derived and complete genotypes at 10 microsatellite loci obtained from fecal samples of 39 individuals from 5 new localities. Two fecal samples were collected from each of these individuals, and the genotype assignment for each individual at each locus was replicated twice (in the case of heterozygotes) or four times (in case of homozygotes) in order to minimize possible genotyping errors due to allelic dropout. We likewise constructed complete 10-locus genotypes for all 25 additional individual samples of unknown provenience for genetic assignment. The complete Genotype database obtained from all localities and clusters of A. caraya studied in this investigation is available at: https://doi.org/10.5281/zenodo.3660723, and the genotypes of all the assigned individuals are presented as Supplementary Table S1.
We did not observe evidence of scoring errors due to stuttering, large allele dropout, or null alleles for any locus in any population in our screening with the program Micro-Checker v2.2.3 32 .
Using Arlequin v 3.5 33 , we did not observe evidence of linkage between any pair of loci (P > 0.05). Significant deviation from Hardy-Weinberg equilibrium was only detected by GenAlEx v 6.5 34 and Arlequin v 3.5 33 softwares for the marker D8S165 in Loc 9 (Piñalito Province Park). Significant evidence of inbreeding (inbreeding www.nature.com/scientificreports www.nature.com/scientificreports/ coefficient: F IS = 0.29) was already found for this population in a previous study 9 . The numbers of different alleles, effective and private alleles, observed heterozygosity (H o ), expected heterozygosity (H e ) and unbiased expected heterozygosity (uH e ) are presented in Table 1.
Analysis using the software Structure v.2.3.4 35 and applying the method described by Evanno 36 to a dataset comprising genotypes already collected for A. caraya in Oklander et al. 9 plus the new localities sampled here showed that the most likely number of different genetic clusters was three (K = 3, Fig. 2). This new analysis resulted in the disappearance of one previously published cluster (EBCO cluster 7 , K = 4, Fig. 2), which in our expanded dataset now clusters with the localities of Chaco National Park, Chaco (Loc 4), Guaycolec, Formosa (Loc 5), and San Alonso, Corrientes (Loc 6, Fig. 2). Accordingly, the complete set of 15 locations sampled were grouped into 3 distinct clusters or regions (Fig. 1). Cluster 1 includes two localities (Locs 1 and 2) from Paraguay-Isla Rio Paraná (P-RP); cluster 2 includes four localities (Locs 3 to 6) from Formosa-Chaco-Corrientes (F-Ch-C); and finally, cluster 3 includes nine localities (Locs 7 to 15) from Misiones-Rio Uruguay (M-RU).
Aiming to evaluate the efficiency of the software GeneClass2 37 for assigning individuals to their localities or clusters of origin, we first tested the individuals included in the complete database whose origins was known. The software correctly assigned 73% of individuals in the database (Quality index 68.73%) when separated according to the 15 localities and 93.3% (Quality index 89.23%) when separated according to the three clusters. We then used GeneClass2 to assign the genotypes of each of the 22 confiscated individuals and 3 corpses of unknown origin to both localities and clusters in the GRDB (Table 2). We established a threshold value to consider the assignment to be reliable if it was higher than 70% 38,39 (Table 2).
For the 25 individuals analysed, 9 individuals' assignment values were below the threshold for assignment to a specific locality. The lack of assignment was expected because not all known populations of A. caraya are represented in the database.
Of the 15 confiscated individuals with assignment values above the threshold, almost half of them (seven) most likely came from Loc 3, two from Loc 2, two from Loc 5, two from Loc 6, one from Loc 1 and one from Loc 13. Additionally, of the 3 corpses found in cities, we could only assign one them, found in Loc 13, to Loc 5, therefore showing a different likely site of origin from the place where was found ( Table 2).
As a means to identify the approximate origins of the individuals who were not assigned to locations, our next step was to evaluate if these individuals could be assigned to one of the three broader clusters described above (Fig. 2). Assignment to clusters allowed the detection of the possible region of origin of 6 additional individuals, one showed a value very close to the threshold level (69.5%), and 1 remained unassignable ( Table 2).
We assigned more than half of the confiscated individuals (14) to 1 of the 3 clusters (F-Ch-C). Four were assigned to the P-PR cluster and 3 to the M-RU cluster.  www.nature.com/scientificreports www.nature.com/scientificreports/ There was a great deal of variation in the extent to which individuals housed at particular rehabilitation centre actually came from nearby clusters. Specifically, the centre Güirá-Oga, in Misiones province, housed only 2 individuals from the nearby cluster (M-RU where they were reintroduced), while 11 belonged to the cluster F-Ch-C and other 4 to the P-RP clusters ( Table 2). Of the 5 confiscated animals housed at Esmeralda, only 4 could be assigned to clusters. Three of them likely came from the nearby F-Ch-C cluster and one from the M-RU cluster. The confiscation sites of these 5 individuals were registered, and the animals that clustered in F-Ch-C (individuals Esmeralda 2, and 4, Supplementary Table S1) were confiscated in Santa Fe province, while the other two (Esmeralda 1 and 5) were confiscated in provinces that do not belong within the natural distribution of A. caraya (Supplementary Table S1). Regarding the 3 corpses found in cities, we could assign one that belonged to a cluster in the same region of Argentina where the city lies, and another one that belonged to a different cluster    www.nature.com/scientificreports www.nature.com/scientificreports/ (F-Ch-C, Table 2). In summary, of all 25 analysed individuals, 15 were inferred to have come from sites within the F-Ch-C cluster.
Finally, of the 17 individuals reintroduced (12 in Misiones and 5 in Santa Fe) only four were reintroduced into a site in their cluster of origin (1 into Misiones, 3 into Santa Fe, Table 2, Supplementary Table S1).

Discussion
Molecular genetic studies have allowed the identification of species confiscated from the illegal wildlife trade, as well as traded animal, timber, and wood products. These applications are even used for species identification from limited samples of body parts (e.g., teeth, feathers, processed tusks), allowing a reliable assessment of the effects of exploitation and the conservation needs of species that would be impossible otherwise. This is the case of the genetic assignment of 28 ivory samples from different elephant populations in Africa 14 between 1996 and 2014 that resulted in the identification of two major poaching hotspots, or determination of the species of dried shark fins being sold in Asian and Mediterranean commercial markets 40 , allowing the monitoring of trade for conservation assessment. Moreover, a study conducted in South Africa revealed that products labelled as "game meat" belonged to domestic species in 76.5% of cases 41 . Thus, molecular analyses are helpful for poaching detection, traffic route identification, and other crimes involving wildlife.
In a study on tortoises 12 , where researchers also developed a GRDB and correctly assigned 90% of the individuals in that database to their population of origin, the lack of assignment of confiscated individuals was attributed to the fact that they came from different locations than the sampling sites included in the GRDB. Nevertheless, the researchers were able to determine that all the confiscated individuals came from the same population. This other approach of genetic assignment of living animals (mostly confiscated or in captivity) shows how genetic tools can be used by wildlife managers to identify the most probable regions of origin of individuals as well as to determine the genetic appropriateness of potential recipient populations when designing reintroduction projects. Translocations have been used to mitigate population decline and restore locally extinct populations 6 . In these cases, genetic data are necessary to guide the selection of populations of origin to which translocated individuals should be released and subsequently evaluate the success of the restoration 2,42,43 .
In this study, the first application of the GRDB of howler monkeys, our results indicate that the most likely origins of most of the confiscated and surrendered individuals were from the areas around Locs 2 and 3, close to  www.nature.com/scientificreports www.nature.com/scientificreports/ the Argentina-Paraguay border (Fig. 1). Therefore, the largest number of illegally trafficked A. caraya originated in this area.
This area is also the location of northeastern Argentina's largest cities, Chaco and Corrientes, and National Highway 12, the main highway connecting these cities with Buenos Aires. The illegal sale of A. caraya has been reported at several locations along this highway 21,22 . This information supports a possible animal trafficking route that begins in northeastern Argentina and ends in Buenos Aires, where the majority of confiscations occur (10 of 22, Supplementary Table S2 21,22 ). Importantly, most of the confiscations and surrenders occurred in cities outside the normal distribution of the species (17 of 22, Supplementary Table S2), indicating that these animals are not only opportunistically captured by locals, but that these animals are intentionally transferred to urban centres. This example illustrates how genetic analysis helps trace wildlife trafficking routes and hotspots and thus aids in the planning and implementation of more effective control measures.
On the other hand, 15 of the 17 animals that arrived at the rescue centre Güirá-Oga, in Misiones, were assigned to either to the F-Ch-C or P-RP clusters; twelve of these individuals were subsequently reintroduced near this rescue centre on Isla Palacio, where the endemic genetic variation belongs to the M-RU cluster; thus, translocation and reintroduction resulted in the injection of genetic variation from animals belonging to different genetic clusters. The 5 individuals that arrived at Esmeralda were also reintroduced in a protected area of General Obligado, Santa Fe (Fig. 1). Although nearby localities are not sampled in the database, we would expect that of our sample areas, genetic variation in the liberation area would be most similar to F-Ch-C, similar to the southernmost area of the distribution of A. caraya. Of these five reintroduced individuals, three belonged to the same cluster and only one belonged to the M-RU cluster; therefore, translocation and reintroduction of these animals also introduced non-local genetic variability, albeit in a lower proportion.
Our findings highlight the importance of conducting genetic studies prior to the liberation of rescued animals. These results also raise the concern of establishing rehabilitation centres servicing each of the three described clusters that could be considered as management units for A. caraya if the goal is to reintroduce animals into their native populations. Based on this, future work could consist of the genetic assignment of all the individuals that are going to be part of translocation programs, as well as the extension of the GRDB in areas where releases are scheduled. In this way, the genetic variability that would be entering a reintroduction site could be evaluated and possible restoration analysis would be possible afterwards.
Conservation genetics is generally not yet well integrated with other efforts in conservation policies. In Latin America, the practical application of genetic principles for the management of threatened species and in the development and implementation of conservation plans should be emphasized 7 . One possible explanation for this disconnect may be that knowledge obtained from scientific research is often not communicated effectively to the field practitioners and/or the authorities who formulate and enact policies.
As shown in the present study, concrete and measurable genetic data represent a very effective tool to help establish and enforce adequate legislation to curb the loss of biodiversity, generate conservation guidelines, and develop population management strategies that include translocation and reintroduction projects.

Methods
Sampling for GRDB. The GRDB for A. caraya used in this study was built using a previously complied database for many locations in Argentina that contained 143 individuals 9 . The number of individuals in the initial database was increased by adding individuals sampled during a monitoring program for the epidemiological surveillance of YFV and other arboviroses in NHP from several sites in Misiones province in 2017 and 2018 (Table 1 and Fig. 1). We sampled 34 individuals from five newly sampled localities (Locs 11 to 15 in Fig. 1) as well as three individuals from Loc 8 and five individuals from Loc 9, which were resampled localities. In total, these 42 individuals sampled corresponded to 39 new individuals (3 wound up being duplicates of previously sampled in the resampled localities, 1 individual in Loc 8 and 2 in Loc 9; Table 1). Summarizing, we used 178 individuals sampled from 15 localities in the overall GRDB. This included 139 from our previous study 9 (we excluded 4 individuals from 2 localities that shared an allele at each locus and were thus considered first-order relatives) and 39 new individuals. Table 1 summarizes the number of individuals, geographical coordinates of sampling locations, and type of samples analysed.
Sampling for genetic assignments. Twenty-five samples were collected for individual genetic assignment. We received hair samples from 17 howler monkeys arriving at Güirá-Oga in 2017. Twelve of these individuals were later reintroduced into a protected area in Isla Palacio at 25°53′32″S, 54°24′38″W (Fig. 1). We also received five tissue samples from monkeys arriving at Esmeralda. All these individuals were later reintroduced into a protected area in General Obligado, Santa Fe at 28°00′12.7″S, 59°32′42.09″W (Fig. 1). A detailed description of these individuals is provided in Supplementary Table S2. Finally, we analysed three tissue samples from monkeys found dead by local authorities in Apóstoles, Posadas, and San Antonio, all in Misiones province in northern Argentina (Fig. 1).
DNA extraction. Two separate faecal samples per individual were stored at room temperature in 50 ml screw-top tubes containing solid NaCl 44 until DNA extraction (three months to one year later). DNA was extracted from faeces using the QIAamp DNA Stool Mini Kit (QIAGEN, Valencia, USA), according to the manufacturer's protocols with slight modifications. DNA was extracted from tissue and hair samples using standard SDS/Proteinase K digestion followed by phenol: chloroform (1 to 1 volume ratio) organic extraction and Microcon P-100 counter-dialysis filters 45 . (2020) 10:3676 | https://doi.org/10.1038/s41598-020-60569-3 www.nature.com/scientificreports www.nature.com/scientificreports/ Microsatellite amplification. Ten microsatellites, developed for A. caraya or other primates and previously used in studies of A. caraya population genetics, were amplified from each sample: AC14, AC17, AC45, TGMS1, TGMS2, D8S165, D17S804, LL1118, LL157 and AB07 [46][47][48][49] . Genotyping PCRs were performed in a final volume of 25 µl using 5-10 ng of DNA template for tissue samples or 5 µl of the extraction pool from stool samples and included 20 mM Tris-HCl, 50 mM KCl, 1.5 mM MgCl 2 , 0.2 mM each dNTP, 1 U of GoTaq DNA polymerase (Promega, Madison, USA), 1 pmol of each forward primer bearing an M13 tail, 4 pmol of each reverse primer, and 4 pmol of M13 labelled with a fluorescent dye (6-FAM) on its 5′ end following recommendations from previous studies 48,50 . All amplifications were performed in a Gen Amp ABI 9700 machine (Thermo Fisher, Palo Alto, USA). PCR products labelled with different fluorochromes were combined and the amplicons separated by electrophoresis on an ABI PRISM 310 Genetic Analyzer (Thermo Fisher, Palo Alto, USA). Alleles were manually scored by performing a visual inspection of electropherograms after developing the bin panel for each locus in GeneMapper ID-X v. 1.2 (Thermo Fisher, Palo Alto, USA) using HD400-ROX as internal size standard. For DNA extracted from stool samples, PCR and sizing was repeated twice (in the case of a heterozygous genotype call) or four times (in case of a homozygous genotype call) to minimize possible genotyping errors due to allelic dropout 51,52 . We recorded an allele only if it was observed at least twice in different amplifications from the same DNA extract. All amplification assays included negative controls.
Ethics statement. This study was carried out in strict accordance with Argentinean laws for research on NHP and following the recommendations of 'Principles for the Ethical Treatment of Primates' of the American Society of Primatologists (available at: https://www.asp.org/society/resolutions/EthicalTreatmentOfNonHumanPrimates. cfm). We received specific approval to conduct this study by the Consejo Nacional de Investigaciones Científicas y Técnicas (CONICET) of Argentina (no. 11420110100322CO). Sampling permits for all the locations of the original database previously complied were already published 7 . Additional specific sampling permits for the new samples presented in this study were obtained from the Ministry of Ecology, Misiones province, Argentina (Permit number: 9910-00086/17) and from the Ministry of Production Santa Fe province, Argentina (Permit Number: GT 13605). Faecal collection was conducted without capturing animals and therefore does not cause any harm to the studied species.
Statistical analysis. Genotypes were screened for null-alleles and to discriminate between errors in allele frequency estimates caused by null-alleles, allele dropout or stutter bands using Micro-Checker v2.2.3 32 . Numbers of different alleles, effective and private alleles, observed heterozygosity (H o ), expected heterozygosity (H e ), unbiased expected heterozygosity (uH e ) and inbreeding coefficient were computed with the software GenAlEx v6.5 34 for each locus and population. Deviations from Hardy-Weinberg equilibrium (HWE) were assessed by employing an exact test and F IS inbreeding coefficient using Arlequin v 3.5 software 33 . Allelic richness was calculated for each locus in a population using the equation: equation: where Ni represents the number of alleles of type i among the 2N genes, and n is sample size, using the software, Fstat v2.9.4 53 .
The new complete set of samples collected was analysed using non-spatial Bayesian clustering with the Structure v.2.3.4 35 program. A series of 20 independent runs per K (ranging from 2 to 6) was conducted using the admixture model with correlated allele frequencies, sampling locations as a prior (LOCPRIOR), and 500,000 Monte Carlo-Markov iterations after a burn-in of 50,000 replicates. The data analysis procedure was further refined using Clump software 54 and a bar plot was constructed with the Disrupt software 55 . The most likely number of K was identified using the method described by Evanno 36 .
Assignment tests give the probability of an individual's multilocus genotype of belonging to reference populations/locations or clusters 56 .
We used Bayesian methods computed by the software GeneClass2 37 to assign the origin of confiscated individuals of unknown origin into the 15 potential locations sampled and into the different three clusters described here using a leave one-out procedure, excluding self-assignment and being 0.05 the assignment threshold score (or p-value, set by the program for default). After testing all the combinations of approaches presented by the program, we chose the Bayesian criteria described by Rannala & Mountain 57 , resulting in a higher quality index and the highest number of correctly assigned individuals when tested against the database. This program is suitable for our study because it does not assume that all populations of origin have been sampled.
For each individual, the best matching location and cluster were sorted and a score in percentage was obtained. The score of an individual, i, in a population, T, is computed as follows: Once this probability of assignment of each individual to a certain location or cluster was obtained, we established a threshold value above which assignment was considered reliable. The criteria were that the scoring for a given population was superior to 70%. This threshold was based on our correctly obtained re-assignments of 73% of the database individuals to their locality and in previous works on assignments for fishes and birds 33,39 .