We collected 74 individuals of A. gephyra (30 from the Cuieiras River and 44 from the Tarumã-Mirim River), and 39 individuals of A. pertensis (27 from the Cuieiras River and 12 from Tarumã-Mirim River) (Table 1), thus a total of 113 specimens. The fishing efforts were similar in both rivers (Figure 1). The taxonomic identification of both species was carefully performed by two taxonomists from INPA’s ichthyology group. The spinal cord of all animals was severed and the whole animal was placed in microtubes with 70% alcohol and transported to the Laboratory of Ecophysiology and Molecular Evolution (LEEM/INPA), where they were stored in a freezer at −80^°C for posterior DNA extraction. The collection permits were provided by the Instituto Chico Mendes de Conservação da Biodiversidade (ICMBio permit #34130-1). Sampling was carried out according to the Brazilian guidelines for animal care and authorized by the Committee on the Ethics of Animal Use and Care—CEUA of the Brazilian National Institute for Amazonian Research—INPA (approval number 048/2012).

The total DNA was extracted of muscle tissue using the phenol-chloroform-isoamyl alcohol protocol described by Sambrook et al. (1989). For amplification, we used 11 microsatellite primers for A. gephyra and A. pertensis (Leitão et al., 2017). To sequence the partial mitochondrial genes, we developed two sets of primers for the genes 16S rRNA and cytochrome oxidase subunit I (COI), based on the partial sequences of these mitochondrial genes for the Apistogramma species available at GenBank—National Center for Biotechnology Information (NCBI). For the cytochrome b gene (Cyt-b), We used the primers developed by Taberlet et al. (1992) and Lydeard and Roe (1997), as suggested by Farias et al., 2001 (Supplementary Table 1). The polymerase chain reactions (PCR) conditions were as follows: for a 25 μL total volume, we used 4.0 μL of total DNA (50 ng μL^–1), 1.5 μL of forward primer (5 μM), 1.5 μL of reverse primer (5 μM), 5.5 μL of H₂O, and 12.5 μL PCR Master Mix 2X [Taq DNA polymerase (0.05 U), MgCl2 (1.5 mM), and dNTPs (0.4 mM)] (Fermentas, Waltham, MA, United States). Amplifications were performed with the following PCR profile: initial denaturation at 94°C for 2 min, followed by 35 cycles with denaturation at 94°C for 1 min, annealing at 50–51°C for 1 min, polymerization at 72°C for 1 min, and a final extension at 72°C for 5 min. We confirmed the amplified PCR products by electrophoresis in a 1.5% agarose gel using GelRed^® (Invitrogen, Waltham, MA, United States), and we visualized it using an L-PIX Molecular Image transilluminator (Loccus Biotecnologia, Cotia, Brazil).

The genotyping of the microsatellites and the sequencing of the mitochondrial genes were performed in a genetic analyzer (ABI 3130xl, Applied Biosystems, Waltham, MA, United States) using GeneScan Liz-500 (–250) size standard to determine the fragment length of the microsatellite, and the Big Dye Terminator v3.1 kit for sequencing (Applied Biosystems, Waltham, MA, United States. The microsatellite alleles were scored based on the consistent pattern of their stutter peaks and on the peak intensity corresponding to each individual at each locus using GeneMapper v4.0 (Applied Biosystems, Waltham, MA, United States). The quality of the sequencing was verified using the 3,130 series data collection software v4 and Sequencing Analysis software v6 (Applied Biosystems, Waltham, MA, United States), respectively.

Using the following allelic diversity parameters of A. gephyra and A. pertensis populations, we calculated the total number of alleles (TNA), mean number of alleles (MNA), number of private alleles (NPA) using GDA v1.1 software (Lewis and Zaykin, 2002). Allelic richness (AR) was acquired using FSTAT v2.9.3.2 software (Goudet, 2001), and the number of effective alleles (NEA) was calculated using the mathematical model: NEA = 1/(1–H_E). We estimated the genetic diversity through the observed and expected heterozygosity (H_O and H_E, respectively) in Hardy–Weinberg equilibrium, fixation index (F_IS) using GENETIX v4.05.2 software (Belkhir et al., 2009).

We applied the mixture model (Admixture) ancestor to analyze the structures of the populations, which correlates the gene frequency among the studied populations. This was performed using STRUCTURE v2.3.1 (Pritchard et al., 2000; Falush et al., 2003) with a burn-in of 50,000 followed by 200,000 steps using the Markov chain Monte Carlo method (MCMC). Each analysis was repeated 10,000 times from a different randomly selected starting point, and convergence between independent runs was assessed via examination of α values and profile of posterior probabilities. The Q values from each of the 10 independent runs for each K scenario were extracted using the program STRUCTURE HARVESTER 0.6.92 (Earl and von Holdt, 2012) and summarized in the program CLUMPP 1.1.2 (Jakobsson and Rosenberg, 2007). Results were visualized in the program DISTRUCT 1.1 (Rosenberg, 2003), and a run for each analysis was performed for K = 1 and K = number of population samples + 2, with ten replicates for each K. We applied the value of ΔK suggested by Evanno et al. (2005) to identify the highest level of genetic division hierarchy.

The levels of genetic differentiation were analyzed using Wright’s F-statistics F_IS and F_ST (Weir and Cockerham, 1984) using ARLEQUIN v3.5 software (Excoffier and Lischer, 2010), and the genetic distance (GD) was calculated using GENETIX v4.05.2. We also calculated the number of migrants (Nm) between populations by applying the F_ST values and the molecular variance (AMOVA) using the ARLEQUIN v3.5 software at the significance level of 5%.

The sequences were edited and aligned using Geneious R7.1.8 (Biomatters Ltd., Auckland, New Zealand) and the Clustal W algorithm (Larkin et al., 2007). The genetic diversity, the total number of haplotypes (NH), the haplotype diversity (h) and the nucleotide diversity (π) (Nei, 1987) was obtained using DnaSP v5.10.1 (Librado and Rozas, 2009). Tajima’s D test (Tajima, 1989) and Fu’s FS test (Fu, 1997) were used to verify the occurrence of recent population size changes.

An analysis of molecular variance (AMOVA) was performed in ARLEQUIN v3.5. Population structure analysis was performed using the Bayesian approach implemented by BAPS 4.14 (Corander et al., 2007), in which the number of clusters was defined using an algorithm that estimates the distribution of haplotype frequencies of all populations tested. For these analyses, we applied two approaches: (1) two groups according to the rivers, to verify the segregation of the two species populations; and (2) two groups according to the species identification, to identify whether their groups are distinct and disconnected, and whether they present differences as populations between the rivers.

The analysis of allelic and genetic diversity between the populations of A. gephyra and A. pertensis from the Cuieiras River (CR) and Tarumã-Mirim River (TMR) revealed a similar mean number of alleles (MNA) (Table 1). The number of private alleles (NPA) ranged from 5 to 9 for A. pertensis and A. gephyra from the Cuieiras River, and 13–25 for A. pertensis and A. gephyra from the Tarumã-Mirim River, respectively.

The results of the structural analyses for the microsatellite data clearly show the separation of the two species within each studied river (K = 2) and that ΔK = 430 (Supplementary Figure 1) yield a clear separation between the species analyzed. The same result was observed with DNA mitochondrial data using BAPS. For K = 3 and K = 4, we observed the formation of two distinct populations for each species between the Cuieiras and Tarumã-Mirim Rivers, respectively (Figure 2). For K = 4, a possible sign of hybridization is suggested between the species A. gephyra and A. pertensis in the Cuieiras River, but in order to verify the actual existence of such hybridization process, it would be necessary to increase the number of samples, as well as increase the number of microsatellite markers.

The genetic distance (GD) between populations is compatible with the expected result for the genetic differentiation index (F_ST), and higher values are described between the populations of A. pertensis from the two rivers (0.228 GD and 0.285 F_ST) when compared to A. gephyra populations from the two rivers (0.120 GD and 0.128 F_ST), respectively (Table 2). The number of migrants (Nm) was 0.85 between A. pertensis from the Tarumã-Mirim River and Cuieiras River and, interestingly, 1.83 between the populations of A. gephyra from the two rivers.

The AMOVA analysis (Table 3) indicates that there are significant differences among individuals of the same population and between different populations for both species. These results are consistent with the F_ST values (Table 2), which show structure among different populations, as corroborated by the significant genetic structure when all sampled individuals were analyzed using the STRUCTURE program (Figure 2).

The number of sequences obtained for the three mitochondrial genes were unequal among populations: 69 sequences for 16S, 62 for COI, and 112 for Cyt-b (see Table 4). Despite this, the size of the amplicons and the quality of the DNA were acceptable for our main goal, which was to verify the structure of the two populations from the two species supposedly separated by the river capture. Due to this unequal sampling, we performed population analysis for each mitochondrial gene separately. The analysis of molecular variance (AMOVA, p < 0.001) for both species (A. gephyra and A. pertensis) indicates high variation within the species (Table 5), confirming that populations from both rivers are structured within each species for all genes.

Haplotype diversity (h) for Cyt b and COI genes was similar for A. gephyra and A. pertensis in the two rivers. For the 16S gene, lower values (0.209 and zero) were observed for the populations of A. gephyra and A. pertensis, respectively, in the Tarumã-Mirim River when compared to the Cuieiras River. The genes 16S and COI showed higher nucleotide diversity (π) in the Cuieiras River for both species when compared to the Tarumã-Mirim River. However, a low nucleotide diversity (π) was found in the Cyt-b gene in both A. gephyra and A. pertensis (Table 4).

The neutrality tests of Tajima (D) and Fu (Fs) showed negative values in the two populations of both rivers, and one positive value was observed in the Fu test for A. pertensis (Cyt-b gene). The neutrality tests were not significant for the Cyt-b gene in either species (Table 4). Therefore, the Tajima test indicates an excess of recent mutations in both populations of Apistogramma, while the Fu test did not suggest any recent population expansion in either of the A. pertensis populations based on the Cyt-b gene. The results obtained for the pairwise genetic differentiation index (F_ST) among the populations of dwarf cichlids showed significant statistical values (P < 0.001), with the only exception being detected for the COI gene in A. pertensis populations, which presented no difference in F_ST values (P > 0.05). The indexes for the genes 16S, COI, and Cyt-b were, respectively, 0.721, 0.356, and 0.352, for A. gephyra from the Cuieiras River and the Tarumã-Mirim River; and 0.377, 0.280 and 0.339 for A. pertensis from the Cuieiras and Tarumã-Mirim Rivers (Table 5). This analysis also indicates high levels of population genetic structure among them.

Analysis of population structure in BAPS based on 16S rRNA and COI genes clearly shows the separation of the two species between the Cuieiras and Tarumã-Mirim Rivers (Figures 3A,C). From the Cyt-b data, we observed a separation of all A. pertensis specimens, though we observed several clusters for A. gephyra (Figure 3E). When we analyzed the species separately, the results showed a clear differentiation between the rivers for both species. On the other hand, when we observed the Apistogramma populations for the Cuieiras River, we verified two or three clusters formed for the mitochondrial genes, and only one cluster in the Tarumã-Mirim River (A. gephrya—16S rRNA, COI, and Cyt-b; and A. pertensis—16S rRNA and COI, respectively) (Figures 3B,D,F). A similar result was observed for the microsatellite data after the analysis using the STRUCTURE program, by which we verified a higher gene flow between both species occurring in the Cuieiras River. In regards, the BAPS analysis corroborates the STRUCTURE analysis for microsatellite data (Figure 2).

Geological events may have been important in the diversification of the genus Apistogramma (Römer, 2006), as well as the biology of this genus, and the habitat characteristics influence the genetic structure of their populations (see Figures 4A–C). Thus, understanding genetic structure is critical in evolutionary and conservation biology, and discriminating fish sub-populations is essential for detecting possible reproductive isolation (Santos et al., 2016).

The results from the microsatellite data indicate low levels of gene flow between populations of the same species from the different rivers (Table 2). These results suggest diversification among populations within each species and the possibility that a local adaptive process is occurring. On the other hand, we cannot discard the possibility that genetic drift caused allele fixation. In fact, the occurrence of local adaptation and random allele fixation by genetic drift leads to a reduced number of genotypes (Blanquart et al., 2012). Moreover, local adaptation is also sensitive in small populations and to the balance between gene flow and local selection.

When gene flow is limited, specialized genotypes can be maintained in isolated populations and support local adaptation, thus providing important insight into evolutionary processes and adaptive divergence of populations (Blanquart et al., 2012, 2013). A key prerequisite for the emergence of local adaptation is the existence of a spatially heterogeneous environment, which generates a heterogeneous selective pressure (Blanquart et al., 2013). The region where the Cuieiras and Tarumã-Mirim Rivers are located is constantly changing, which involve this type of study. In addition, the drainage network of the Amazon is not static but is, in fact, constantly changing over time (Albert et al., 2018). All those factors induce allele fixation in these populations, genetic drift and local adaptation to habitat heterogeneity, thus making this a perfect ecosystem for evolutionary studies.

Interestingly, the nuclear data indicates genetic similarity between the two species of the Cuieiras River. However, this pattern is not seen in the mtDNA data and the two species form a well-supported high differentiation. We considered that these findings indicate that gene flow occurs in the Cuieiras River, which affects only the nuclear genome, though not the maternally inherited mtDNA. To have greater certainty regarding these processes, further studies must explore a higher and more variable number of markers (e.g., SNPs) and use more rapidly evolving mitochondrial genes (e.g., mtDNA control region).

Thus, the lack of support for intraspecific clades in mtDNA may be due to the slow evolutionary rate of the markers utilized in this study, while the nuclear microsatellites clearly show an intraspecific differentiation in the different drainage systems. Another explanation would be periodic connectivity or delayed disconnection between the Tarumã-Mirim and Cuieiras Rivers, as observed in ongoing drainage captures elsewhere in the Amazon (Stokes et al., 2018). This would allow eventual gene flow among sub-populations (or evolutionary significant units—ESU). Although the sample size was unequal for the populations of Tarumã-Mirim and Cuieiras Rivers, the level of polymorphism for both the microsatellites (Leitão et al., 2017) and the number of mitochondrial genes was substantial enough to reach these conclusions.

Tajima’s D test results in negative values after a recent bottleneck, which indicates population expansion (Fu. 1997). The statistical values for Fs are also very sensitive to deviations from neutrality that are promoted by demographic changes. Herein, they resulted in negative values under most conditions and were positive only for A. pertensis (Cyt b).

The distribution of Apistogramma species is higher in Cuieiras River than the drainage basin of the Tarumã-Mirim River (Figure 1). This is consistent with an increase in species range as the Cuieiras River has gained drainage area due to river capture. Thus, a larger gene flow has been detected between species of dwarf cichlids from the Cuieiras River, possibly due to mating occurring at higher levels, though with production of infertile offspring. We suggest that the mating observed between the two dwarf cichlids species in the Cuieiras River occur because these two species have not differed enough to prevent reproductive mating after the river capture event, which means they might still maintain a correlation with the original populations (Hypothesis 1, Figure 4B). In addition, both species from Tarumã-Mirim River have been genetically disconnected because they inhabit a fragmented and smaller area, and this may have caused higher environmental pressure on these populations. In addition, we must emphasize that these genetic studies are contemporary, and that the current landscapes may have a strong influence on our results. Thus, river captures may be considered a mechanism that promotes fish biodiversity, particularly in non-migratory fish species, and contributes to the genetic diversity at the population level, as in the present case (Albert et al., 2017).

An alternative explanation for the gene flow in the Cuieiras River is presented in Hypothesis 2 (Figure 4C): for the same reason as above, i.e., the populations in the Tarumã-Mirim River after river capture are the original populations prior to river capture (Figure 4C). Thus, our preferred explanation, prior to capture, is that these two species remained there, though with no gene flow.

The current pressure that affects the Tarumã-Mirim River is due to the urban expansion of Manaus, which probably affect these populations. The conservation approach over these populations must be taken to avoid disappearance in the future. In addition, the population of the Cuieiras River has evolved in a larger population with more landscape changes over time, which means that the river capture induced all the changes in this river and left the Tarumã-Mirim River behind. The latter remained a short (Figure 1) river and is the final part of the original river.

Based on the results, we are able to affirm that, based on the structured populations of both species in both rivers, these two species have given rise to four populations, and that these four populations will most probably generate new species through interruption of gene flow.

The current literature is of the consensus that a mosaic of events must have occurred at different times and spatial scales and that this created a complex scenario of spatial and temporal landscape changes that have influenced the biodiversity in the Amazon we see today (Hoorn et al., 2010, 2017; Rull, 2011, 2013; Gorini et al., 2014; van Soelen et al., 2017). Our study suggests that the populations of Apistogramma were affected by landscape changes that contributed to speciation and, eventually, species enrichment. Thus, even at smaller spatial scales, local landscape changes are important in creating biodiversity.

The cichlids are extremely interesting fish not only for adaptive radiation studies, but also for studies of behavior and reproductive characteristics (Turner, 2007). About 112 dwarf cichlids species are described as being part of the Apistogramma genus, and wide distribution throughout South America (Costa et al., 2019), particularly in the Amazon basin. This extraordinary diversity and high endemism rate found for dwarf cichlids has been explained by the geological, paleoclimatic, and water-level history that has played an important role in speciation events by isolating populations and favoring the establishment of reproductive barriers (Römer, 2006). Thus, new conservation strategies for these species and ecosystems should be rethought; mainly in relation to the vulnerability of the population structure of Apistogramma species that live in small basins and tributaries of large rivers to extreme anthropic changes.

Thus, for conservation measures to be effective, they must be based on the information between environment and organism, in order to avoid local extinction of threatened species. The integration of some types of analyses, such as species distribution, phenotypic variation, landscape change, and population genetics, can contribute to management and conservation strategies (Santos et al., 2011, 2016; Escobar et al., 2015; Gravena et al., 2015). In addition, knowledge of the biology of the species is vital for defining a conservation strategy, particularly in areas under environmental pressure. Most species of the genus Apistogramma are restricted to a single river or adjacent floodplains and nearly never found in the main river channel (Kullander, 2003; Römer, 2006).

Genetic population studies have constantly revolutionized our theories regarding evolutionary processes (De la Ossa-Guerra et al., 2020), and genetic variability in this case is considered an important factor for understanding these evolutionary processes and for planning conservation programs (International Union for Conservation of Nature [IUCN], 2020). Thus, understanding population structure is crucial for management actions for evolutionary significant units (ESU); however, developing conservation strategies for establishing priority units has historically been difficult. The low migration rate and the high population structure imply that A. gephyra and A. pertensis populations from Cuieiras and Tarumã-Mirim Rivers should be managed separately at the moment, each as an independent ESU (Avise, 2000; Baker et al., 2002) and, therefore, taxonomic units have to be conserved considering these aspects (De la Ossa-Guerra et al., 2020). This operational concept is often used to guide short-term management strategies, as it refers to current population structures and allele frequencies (Moritz, 1994).

The two dwarf cichlid species analyzed in this study (A. gephyra and A. pertensis) belong to two drainage networks that were separated by a river capture event. Based on our results, the paleovalley that now separates the Cuieiras and Tarumã-Mirim River basins is seen to have contributed to the local interruption of gene flow, and led to the development of structured populations of both species. According to De la Ossa-Guerra et al. (2020), in order to preserve evolutionary processes, species management must preserve natural networks of genetic connections among populations, rather than just the isolated populations within that network. On the other hand, our results support that geological changes induce new river courses that may cause genetic diversification in aquatic fauna, thus increasing its biodiversity. Moreover, given that events resulting in landscape changes, such as the one studied here, have been common throughout the Amazon basin during and after the Miocene (Rossetti et al., 2005, 2016; Almeida-Filho and Miranda, 2007; Hayakawa et al., 2010; Hoorn et al., 2010; Hayakawa and Rossetti, 2015), we may theorize that small geological processes are capable of producing vicariant events and should be considered within the complex framework of biodiversity drivers in the Amazon Basin. Further studies will clarify the ongoing population structure of A. gephyra in the Cuieiras River to better explain the influence of other surrounding geological events, such as the formation of the Anavilhanas Archipelago and the influence of the main channel (Negro River), on the emergence of new genetic information in this species.