Black-and-gold howler monkeys were sampled from three countries: Argentina, Paraguay, and Brazil, with a total of 256 individuals in 22 populations (Table 1 and Figure 1). Tissue or blood samples were obtained by the capture of adult individuals in the context of EIA or during flooding for dam construction in four sites (GO.BR, MT.BR, MS1.BR, and MS2.BR; Table 1 and Figures 1, 2). Consequently, we are analysing individuals that were relocated or dead.

Samples from GO.BR were collected in the EIA of Serra da Mesa Hydroelectric Power Plant, located in the Alto Tocantins Basin, Goiás state, Brazil (GO). This dam is located on the main course of the Tocantins River, and its reservoir is the largest in Brazil in terms of water volume, with 54.4 billion m³ and an area of 1,784 km². Samples from MT.BR were collected in the EIA of the Represa de Manso hydroelectric dam located on the Manso River, a tributary of the Cuiabá River, Mato Grosso state, Brazil. The reservoir has a volume of 7.3 billion m³ and an area of 427 km². During construction, many human families were displaced from the area now covered by the reservoir, and the submerged, decaying vegetation has had a negative impact on water quality and has affected edible fish stocks (Hylander et al., 2006). Samples MS1.BR and MS2.BR were collected during the flooding of the Engenheiro Sérgio Motta Hydroelectric Power Plant (also called Porto Primavera Hydroelectric Power Plant) on the Paraná River, 28 km upstream from the confluence with the Paranapanema River. This dam lies between the Municipality of Porto Primavera in the eastern part of the Mato Grosso do Sul (MS) state and southwestern-most São Paulo state. This dam has the greatest surface area in Brazil, being 10,186.20 m long with a reservoir formed by water from the Paraná River that has a flooded area of 2,250 km². The other three samples’ sites in Brazil belong to the states of Paraná (PR), São Pablo (SP), and Minas Gerais (MG), in areas where no dams have been built.

Tissue or faecal material was obtained from monkeys in sites where the habitats remain intact and accessible for sampling (Table 1; Oklander et al., 2017, 2020). For some analyses, the following study sites in Brazil GO.BR, MT.BR, MS1.BR, and MS2.BR were grouped as “Flooded.” The Brazilian sites PR.BR, SP.BR, MG.BR, and all sites in Paraguay and Argentina were grouped as “Remaining” (either Rem.BR or Rem.PY.AR; Table 1).

Sampling complied with Argentine, Paraguayan, and Brazilian laws and was conducted with permission from the National Resources Board (Brazil: SISBIO n. 19927-12 by ICMBio, Paraguay: Ministerio de Salud Pública, Mercosur/GMC/Res. n. 50/08, Argentina: CONICET n. 11420110100322CO and Ministry of Ecology, Misiones province, Argentina n. 9910-00086/17). We adhered to the American Society of Primatologists’ Principles for the Ethical Treatment of Primates.¹

DNA was extracted from faecal samples using the QIAamp DNA Stool Mini Kit (QIAGEN, Valencia, United States), according to the manufacturer’s protocols, with slight modifications (Oklander et al., 2017). DNA was extracted from tissue and blood samples using standard SDS/Proteinase K digestion followed by phenol:chloroform organic extraction and Microcon P-100 counter-dialysis filters (Green and Sambrook, 2012).

Ten microsatellites or Short Tandem Repeats (STRs) characterised for A. caraya and previously used in population genetics studies, were amplified from each sample: AC14, AC17, AC45, TGMS1, TGMS2, D8S165, D17S804, LL1118, LL157, and AB07 (Tomer et al., 2002; Di Fiore and Fleischer, 2004; Gonçalves et al., 2004; Oklander et al., 2007; Supplementary File 1). Genotyping PCRs were performed using recommendations from previous studies, with a final volume of 25 μl using 5–10 ng of DNA template for tissue samples or 5 μl of the extraction of faecal samples, 20 mM Tris–HCl, 50 mM KCl, 1.5 mM MgCl2, 0.2 mM each dNTP, 1 U of GoTaq DNA polymerase (Promega, Madison, United States), 1 pmol of each forward primer bearing an M13 tail, 4 pmol of each reverse primer, and 4 pmol of M13 labelled with fluorescent dyes (6-FAM, HEX or NED) on its 5′ end (Schuelke, 2000; Oklander et al., 2007). Singleplex PCRs were carried out for faecal samples and multiplex of 3 or 4 STRs for tissue and blood samples (Supplementary File 1).

All amplifications were performed in a Gen Amp ABI 9,700 machine (Thermo Fisher, Palo Alto, United States). PCR products labelled with different fluorochromes were combined and co-injected with HD400-ROX as an internal size standard to be separated by electrophoresis on an ABI PRISM 310 Genetic Analyzer (Thermo Fisher Scientific, Palo Alto, United States) or in MegaBACE 1000 (Pharmacia, Upsala, Seden) automated sequencer (Amersham Biosciences) according to the protocols of the manufacturer. 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, United States) or software Genetic Profiler 2.2 (GE Healthcare).

For DNA extracted from faecal samples, PCR and sizing were repeated two times (in the case of a heterozygous genotype call) or four times (in case of a homozygous genotype call) to minimise possible genotyping errors due to allelic dropout (Pritchard et al., 2000; Peakall et al., 2003). We only recorded an allele if it was observed at least two times in different amplifications from the same DNA extract. All amplification assays included negative controls. Every step of the lab work was carried out in specific laboratory spaces (DNA extraction room and PCR room), inside laminar flow hoods with negative pressure, and using aerosol-resistant filter tips to avoid between sample cross-contamination.

We amplified a 507-bp fragment of the left domain of the mitochondrial DNA Control Region (mtDNA, CR) using primers How RA-1 (5′-CTACCATCAACACCCAAAGC-3′, Ascunce et al., 2007) and RC-BugioR (5′-CCAGGTTAAGAGGGTGATAGC-3′, Oklander et al., 2017). The amplified zone corresponds to the position 15,412–15,903 bp from the mitogenome published by Finstermeier et al. (2013) and deposited at Genbank: NC_021938. Amplification reactions were performed in a final volume of 25 μl containing 25 ng of DNA extractions, 1X of GoTaq buffer (Promega, United States), 1.75 mM of MgCl2, 0.2 mM of each dNTP, 1U GoTaq DNA polymerase (Promega, United States), 4 pmol of each primer, and 4 mg/ml of bovine serum albumin (BSA). Cycling parameters were an initial denaturation (5 min at 94°C) followed by 35–40 cycles of denaturation (1 min at 94°C), annealing (30 s at 50°C), and extension (1 min at 72°C), then a final extension step (3 min at 72°C). Sequencing was carried out using the same reverse and forward primers used for initial amplification. All samples were sequenced bidirectionally using the BigDye^® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA) and visualised in an Applied Biosystems 3,500 Genetic Analyzer or the DYEnamic™ ETDye Terminator Cycle Sequencing Kit (GE Healthcare) on a MegaBACE 1,000 automated sequencer. Sequences quality was assessed with Sequencher 5.3 (LifeCodes, United States).

Sequences were edited and aligned using the CLUSTAL W algorithm (Thompson et al., 1994) in MEGA10.2.6 (Kumar et al., 2018).

Haplotype frequencies, the number of polymorphic sites (S), haplotypic (Hd), and nucleotidic (Pi) diversity were calculated in DNASP 5.0 (Librado and Rozas, 2009). Standard tests of selective neutrality and Tajima’s D (Tajima, 1989) and Fu’s F_S (Fu, 1997) with 95% confidence intervals were computed in DNASP with 1,000 simulations and a neutral infinite-sites model assuming a large constant population size (Hudson, 1990). All these parameters were performed for each study site and for the three groups (Flooded, Rem.BR, and Rem.PY.AR, Table 2).

A haplotype network was built using Median-Joining (Bandelt et al., 1999), and implemented in POPART (Leigh and Bryant, 2015).

Pairwise Fst’s and Nei’s average number of pairwise differences between populations (Nei, 1972) were determined with Arlequin v3.5 software and plotted with R-lequin (Excoffier and Lischer, 2010). This analysis involved 112 nucleotide sequences, with all ambiguous positions removed for each sequence pair (pairwise deletion option). There were a total of 491 positions in the final dataset from the original 507 pb.

As there is evidence of hybridisation between A. caraya and A. guariba in some of the sampled sites (Aguiar et al., 2007; Agostini et al., 2008), we also carried out amplifications of the same mtDNA fragment in 4 samples of A. guariba (2 from Argentina and 2 from Brazil) in order to obtain the evolutionary relationships for these two taxa. The optimal Neighbor-Joining phylogenetic tree obtained shows reciprocal monophyly between taxons (Supplementary File 2), ruling out that the differences found in our study are due to hybridisation with this species.

We conducted Mantel tests between genetic Fst (mtDNA and STRs) and geographic distances (in Km, computed from geographical coordinates) using the mantel function (method = “spearman,” permutations = 9,999) implemented in the “vegan” v2.4 ± 3 package (Dixon, 2003) in R. The proportion of shared alleles by pairs of populations was calculated with PopGenReport v2.2.2 package (Adamack and Gruber, 2014) in R. If geographic distances influence gene flow, we would expect the Mantel statistic to be significant (Chiappero et al., 2016).

R-matrix analysis and the regression model of Harpending and Ward (1982) D2 distances derived from R-matrix analysis (Harpending et al., 1973) were employed for gaining insights into the patterns of external gene flow and genetic drift due to isolation under the island model. The broken line is the expected relationship predicted by the model according to the formula hi = H (1-ri), where ri is the distance from the centroid. hi is the heterozygosity of population i, and H is heterozygosity of the total population according to Genocline v1.4 software (Peña et al., 2020). If gene flow from outside the region varies substantially from population to population, this linear relationship no longer holds. Isolated groups will be less heterozygous than the linear prediction and will therefore lie below the theoretical regression line, whereas populations receiving more gene flow from outside will be more heterozygous and lie above the line (Lakshmi et al., 2002). Consequently, the outliers above and below the theoretical regression line may give some insights into the population structure.

For the majority of the samples (n = 243 or 94.9%), all 10 STRs loci were amplified. For the remaining 13 samples, a minimum of nine loci were amplified. We found no evidence of significant deviations from the Hardy-Weinberg equilibrium, with the 256 A. caraya individuals from 22 populations exhibiting moderate levels of microsatellite diversity (Table 1).

The highest population of genetic variability statistics were found in the four flooded study sites. The highest mean number of different alleles (NA), number of effective alleles (Ne), private alleles (PA), and observed heterozygosity (Ho) were found in the two study sites, now submerged by the Porto Primavera Hydroelectric Power Plant in Brazil (MS1.BR and MS2.BR; Table 1), and the highest mean gene diversity (Hs) and allelic richness (AR) in the other two (GO.BR, and MT.BR; Table 1 and Figures 1, 2). On the other hand, the lowest values were detected in the remaining study sites in Paraguay and Argentina (Table 1). A detailed analysis of Argentine and Paraguayan study sites is given in Oklander et al. (2017, 2020).

Mean NA, Ne, Ho, Hs, and AR showed significant differences between groups (AR p = 0.004, Ho p = 0.001, Hs p = 0.001, Na p = 0.0006, Ne p = 0.00009; Table 1).

The BLAST search of the 112 A. caraya 491-bp mtDNA fragment sequences as queries retrieved only mitochondrial sequences. Comparison showed an average identity of 94.20% (Alouatta caraya), 82.44% (Aluoatta guariba, Aluoatta palliate, and Aluoatta seniculus), 77.26% (Ateles sp.), 77.66% (Brachyteles arachnoides), 77.30% (Cebus olivaceus castaneus), 77.24% (Lagothrix sp.), 76.82% (Saimiri), and 77.21% (Sapajus sp.). Therefore, the presence of nuclear mitochondrial DNA sequences (NUMTs) in our fragment sequences can be considered negligible.

The 112 mtDNA sequences were grouped in 42 haplotypes, of which 21 were new to this study (Genbank accession numbers: OK033252-OK033272) and the rest of haplotypes were taken from Oklander et al. (2017; MF095740- MF095754). For the complete dataset, 80 polymorphic (segregating) sites, an overall haplotypic diversity of 0.963 ± 0.007, and an overall nucleotide diversity of Pi = 0.0338 ± 0.001 were observed. The median-joining network showed a complex pattern of relationship among observed haplotypes (Figure 3 and Supplementary File 2). The Flooded mtDNA variation, represented by the haplotypes Hap1-Hap11 (MS1.BR and MS2.BR), Hap20 (GO.BR), and Hap21-Hap22 (MT.BR) are well distributed along the network, showing a higher number of mutation steps among these 14 Control Region variants. For the rest of the remaining populations sampled in Brazil (Rem.BR) and Paraguay-Argentina (Rem.PY.AR), 8 and 21 haplotypes were obtained, respectively, with lower values of variation indices (S, Hd, and Pi; Table 2). Although the overall pattern of the network showed strong haplotype relationships, only one haplotype (Hap8) was shared between MS1 (Flooded) and PIM (Rem.PY.AR). Similarly, the general pattern of the haplotype network supports a relative demographic stability and absence of a recent population expansion, as suggested by the values obtained for the Fs of Fu and the D of Tajima (all p-values > 0.1; Table 2).

Mantel statistic based on Spearman’s rank correlation rho showed significant results for mtDNA (r = 0.3847, p = 0.0023) and STRs (r = 0.4603, p = 0.0002) which implies an inverse proportion of effective dispersal to geographical distance.

However, for “short” distances (< 300 km) a wide range of Fsts is observed (from almost zero to 1; Figure 4), perhaps related to what is observed in the haplotype network. Particularly, close study sites with high genetic differentiation, up to ∼30 difference mutations (Hap 1 and 11 from MS1.BR and MS2.BR), and low differentiation (2 mutations) between haplotypes in distant sites (Hap 11 from MS2 and Hap 19 from MG.BR; Figure 3).

The STRUCTURE results identified a genetically different number of clusters (K) depending on the method. The parsimony method only detected one cluster grouping in the 22 study sites (Figure 5; K = 1). The Evanno method results showed K = 2 with a clear separation between the remaining study sites of Argentina and Paraguay (Rem.PY.AR, represented as green in Figure 5; K = 2) from all other studied sites—all of which were in Brazil (represented as blue in Figure 5; K = 2). Finally, the Pr[X| K] method results detected five genetically differentiated clusters (K = 5) as best explaining the nuclear genetic variation observed in A. caraya. Two of these clusters are mostly represented by the four-flooded sites for dam construction coloured as yellow and blue in Figure 5, K = 5 (Flooded: GO.BR, MT.BR, MS1.BR, and MS2.BR), although the yellow cluster is still represented in the remaining MG:BR study site. The third cluster is represented in orange in Figure 5, K = 5, and is mostly represented in the remaining study sites in Brazil (PR.BR and SP.BR). Finally, the remaining Argentine and Paraguayan study sites were grouped in two clusters (green and grey) with this method (Figure 5; K = 5). These study sites, grouped as Rem.PY.AR, were previously studied in Oklander et al. (2017, 2020) and remain currently relatively undisturbed since no major infrastructure projects or deforestation have been carried out in these sites since sampling occurred.

Pairwise Fst’s and Nei’s were used to further investigate the genetic differences between flooded (Flooded) and remaining sites (Rem.BR and Rem.PY.AR), with both markers revealing different gene flow situations in time (Figure 6). The mtDNA analysis showed that the flooded study sites showed the highest genetic diversity estimates and shared one haplotype with one Argentine site (Table 2 and Figure 3). In addition, the highest average pairwise between-population differences occur between the Flooded group and the two remaining groups (Flooded vs. Rem.BR p = 0.009; Flooded vs. Rem.AR p = 0.0000; Rem.PY.AR vs. Rem BR p = 0.0000; Figure 6A), suggesting that they historically exchanged migrants with other populations. The remaining populations of Brazil showed an intermediate variability (Table 1) and differentiated mtDNA haplotypes from both the flooded populations and the remnants of Argentina, suggesting isolation processes between all the sites (Table 2 and Figure 3).

On the other hand, STR analysis showed that the major mean within-population pairwise difference occurs in the Rem.BR group, and the major average pairwise differences between populations are from both Flooded and Rem.BR groups compared to the Rem.PY.AR group. Meanwhile, Nei’s distance shows greater values for both molecular markers between the two Brazilian groups (Flooded and Rem.BR) compared to the Argentine and Paraguayan sites (Flooded vs. Rem.BR p = 0.1273; Flooded vs. Rem.AR p = 0.0000; Rem.PY.AR vs. Rem BR p = 0.0000, Figure 6B).

Tajima’s Neutrality Test comparison of the mtDNA from three groups (Table 2) again showed significantly higher parameter values for the Flooded group compared to those in the remaining groups. According to R-matrix analysis and the regression model of Harpending and Ward (1982), the Brazilian study sites received more gene flow as they have greater heterozygosity and lie above the theoretical regression line. In contrast, Argentine and Paraguayan sites are more isolated and less heterozygous, hence they lie below the theoretical regression line (Figure 7A). When we perform the same analysis grouping in the three groups already mentioned, we see that only the Flooded group lies above the theoretical regression line (Figure 7B).

Genetic diversity is important in ensuring species survival as populations that have high levels of genetic diversity are better able to face environmental change than those that have low levels (Hoban et al., 2020). The sites of A. caraya included in the current study cover a large part of the distribution of the species (Figure 1). Our results showed: (1) genetically differentiated populations between the Brazilian and the southernmost populations of the species, (2) that a large part of the genetic variability and unique genetic lineages of the Brazilian populations were probably lost due to dam constructions, and (3) that remaining studied populations have low genetic variation, perhaps making them more vulnerable to future threats as fragmentation and disease outbreaks.

As our sampling does not cover all the occurrence sites or the entire species range of Alouatta caraya. we cannot assert that the genetic variability lost in the flooded sites is not represented in any other remaining population. However, the densities and group sizes recorded for the species on the upper and middle Parana River (Aguiar et al., 2009) were larger than those reported elsewhere in Brazil (Moura, 2007; Silva and Codenotti, 2007). Therefore, the processes of dam construction described in these areas undoubtedly affected at least the frequency of these variants as a whole for this species. On the other hand, as we have detected that genetic variability is related to geographic distance, the unsampled populations at the extreme west of the distribution are also highly likely to be differentiated from those sampled in this study.

At the mtDNA level, of the 42 haplotypes described in the 22 study sites, 14 belong to the flooded sites, of which only one is shared with one of the remaining sampled sites. In other words, these 13 haplotypes found only in the now submerged areas are absent from any other studied population and have possibly disappeared from the genetic variability of the species.

At the autosomal level, the study sites of A. caraya included in the current study had a mean microsatellite diversity of 0.531 ± 0.095. This is slightly lower than most congeneric species studied with some of the same markers employed here [Alouatta pigra: 0.588 (Van Belle et al., 2012), Alouatta belzebul: 0.640 (Gonçalves et al., 2004), Alouatta palliata: 0.584 (Milton et al., 2009)]. Flooded study sites were the most diverse, as shown by the genetic indexes and the presence of exclusive alleles and haplotypes (Tables 1, 2). If we now consider only the variability of the remaining sites (0.504 ± 0.074), it is likely that, on average, they are genetically impoverished compared to other howler species.

The results from the population structure can be hierarchical as evidenced in the different k-value estimation methods that detect one, two, or five clusters depending on the model used. The Pr[X| K] method results detect five genetically differentiated population clusters (k = 5), of which two are practically gone with the water, with only one poorly represented in the remaining MG:BR study site.

To preserve what is left of the gene pool of A. caraya, we recommend that the remaining four genetic clusters must be given high conservation priority. Our evidence supports significant genetic differentiation between populations, indicating contemporary restricted gene flow between A. caraya geographic distribution. These results have direct implications for the conservation of howlers and should be considered by policymakers in future management plans, as well as during dam construction in the habitat of these vulnerable Neotropical primates.

During the flooding to create the lake for the Engenheiro Sérgio Motta Hydroelectric Plant (sites MS1.BR and MS2.BR in the current study), 4,000 A. caraya individuals were rescued along a 180 km stretch of the Paraná River, with their subsequent release in nearby forest areas (Inglêz, 2006). However, as can be seen in Figure 2, there are almost no remnants of forest on the sides of the flooded area.

The impact of dams on primate populations is not exclusively due to the loss of the flooded forest area. Before and during construction, numerous anthropogenic modifications are also generated in the surrounding areas, such as clearings, arrival of workers, resettlement of local people, and construction of roads and infrastructure, housing, and agricultural needs for these people, all of which diminishes the quality of adjacent forests (Tan and Yao, 2006; Moore et al., 2010; Harrison-Levine et al., 2019). After construction is complete, animals that have remained in the area of the reservoir either drown, remain in treetops that rot and die shortly after, or may try to disperse by swimming. During the flooding of a dam in the Brazilian Amazonia, caiman and jaguars were observed capturing howler monkeys while swimming toward the banks, while eagles captured the animals that remained on emergent trees (Harrison-Levine et al., 2019).

Both “rescued” animals and those arriving in non-flooded areas by their own means are exposed to an unfamiliar environment, difficulties finding resources, such as food, water, and a social group, and a higher risk of stress, disease, predation, and/or human hunting (World Commission on Dams [WCD], 2000; Martin, 2004). There are reports of primates found malnourished, injured, infected with parasites, and stressed during post-flood rescue operations at dams in Thailand (Nakhasathien, 1989), French Guiana (de Thoisy et al., 2001), and Brazil (Harrison-Levine et al., 2019). De Sá (2004) observed an increase in the population density of Plecturocebus, Pithecia, and Samiri in the adjacent habitat just after the release, but an abrupt decrease in the same in subsequent years. Another example, also in Brazil, reports a telemetry follow-up that found that only 7% translocated Cebuella pygmaea survived 3 months after release (Dias et al., 2015).

Historically, our work has aimed to discern the genetic diversity of this species to provide guidelines for the future management of populations (Oklander et al., 2017, 2020, 2021). For all the reasons mentioned above, we assume that the populations analysed in the flooded areas are functionally extinct (Säterberg et al., 2013). Therefore, in the current study, we unfortunately present results of “fossil” genetics. Since the largest number of dams in the future is scheduled in South America (Fearnside, 2006; Zarfl et al., 2015), an awareness of both their long-term detrimental impacts on biodiversity and developing strategies to mitigate such effects should be considered to be of utmost importance, especially in the case of A. caraya, a species demonstrated to be of public health importance. We believe this work provides strong evidence for the need for attention when planning and implementing future hydroelectric dams.

It has already been reported that pre-construction EIA studies fail to: (1) estimate the indirect impact on surrounding areas (Tan and Yao, 2006; Moore et al., 2010), (2) include intra-specific measures of biodiversity, such as genetic variation of species (Simões et al., 2014), and (3) implement management plans to minimise biodiversity losses (Schneider, 2001; Alho, 2011). The occurrence of a severe population declines of the already “Critically Endangered” (CR) Tonkin snub-nosed monkey (Rhinopithecus avunculus) (from 130 to 40 individuals in 13 years, Hoang, 2010) was predicted by an EIA conducted prior to dam construction in Vietnam (Scott Wilson Asia Pacific Ltd [SWAPL], 2000). Despite this, the dam was still built, and few of the recommended mitigation measures were implemented. In the current study, we show an important consequence of the construction of dams for the species A. caraya as the extensive loss of genetic variation, leaving surviving populations greatly genetically impoverished.

To prevent this type of irreversible biodiversity loss from continuing, it should be mandatory for hydropower infrastructure to include all items mentioned in EIA in the final construction. Construction of future hydroelectric plants should explicitly consider appropriate release sites for animal rescues with sufficient resources, minimal anthropogenic threats, and relatively low population densities of conspecifics. In addition, affected primate populations should be part of long-term monitoring programs, including impacts on behavior, health status, and genetics (McAllister et al., 2001).