Riverscape Genomics Clarifies Neutral and Adaptive Evolution in an Amazonian Characin Fish (Triportheus albus)

Hay, Abbie C.; Sandoval-Castillo, Jonathan; Cooke, Georgina M.; Chao, Ning L.; Beheregaray, Luciano B.

doi:10.3389/fevo.2022.825406

ORIGINAL RESEARCH article

Front. Ecol. Evol., 04 March 2022

Sec. Biogeography and Macroecology

Volume 10 - 2022 | https://doi.org/10.3389/fevo.2022.825406

This article is part of the Research Topic The Role of Rivers in the Origins, Evolution, Adaptation, and Distribution of Biodiversity View all 13 articles

Riverscape Genomics Clarifies Neutral and Adaptive Evolution in an Amazonian Characin Fish (Triportheus albus)

$\r\nAbbie C. Hay$ Abbie C. Hay¹

Jonathan Sandoval-Castillo¹

Georgina M. Cooke²

Ning L. Chao^3,4 $Luciano B. Beheregaray*\r\n$ Luciano B. Beheregaray^1*

¹Molecular Ecology Laboratory, College of Science and Engineering, Flinders University, Adelaide, SA, Australia
²Molecular Ecology Lab, Department of Biological Sciences, Macquarie University, Sydney, NSW, Australia
³Bio-Amazonia Conservation International, Brookline, MA, United States
⁴National Museum of Marine Biology and Aquarium, Pintung, Taiwan

Understanding the role of natural selection in the evolution of wild populations is challenging due to the spatial complexity of natural systems. The richest diversity of freshwater fishes in the world is found in the Amazon Basin, a system where marked hydrochemical differences exist at the interface of major rivers with distinct “water colors” (i.e., black, white, and clear water). We hypothesize that divergent natural selection associated with these “aquatic ecotones” influences population-level adaptive divergence in the non-migratory Amazonian fish fauna. This hypothesis was tested using a landscape genomics framework to compare the relative contribution of environmental and spatial factors to the evolutionary divergence of the Amazonian characin fish Triportheus albus. The framework was based on spatial data, in situ hydrochemical measurements, and 15,251 filtered SNPs (single nucleotide polymorphisms) for T. albus sampled from three major Amazonian rivers. Gradient Forest, redundancy analysis (RDA) and BayPass analyses were used to test for signals of natural selection, and model-based and model-free approaches were used to evaluate neutral population differentiation. After controlling for a signal of neutral hierarchical structure which was consistent with the expectations for a dendritic system, variation in turbidity and pH were key factors contributing to adaptive divergence. Variation in genes involved in acid-sensitive ion transport pathways and light-sensitive photoreceptor pathways was strongly associated with pH and turbidity variability. This study improves our understanding of how natural selection and neutral evolution impact on the distribution of aquatic biodiversity from the understudied and ecologically complex Amazonia.

Introduction

Natural selection is a ground-breaking conceptual framework to our understanding of evolution (Darwin, 1859). However, challenges with empirically studying adaptive divergence in nature has ultimately led to a bias toward using spatial factors to explain evolutionary patterns (Wang and Bradburd, 2014). Selectively neutral models of population divergence include Isolation by Distance (IBD), where gene flow aligns with gradients of distance between populations (Wright, 1943; Sexton et al., 2014), and neutral hierarchical structure, where punctuated barriers restrict gene flow to particular geographic zones (Meirmans, 2012). When gene flow and differentiation better align with similar environments, an adaptive Isolation by Environment (IBE) model, where divergent natural selection acts on evolutionary advantageous phenotypes, might provide a more superior explanation for inferred patterns of population structure (Sexton et al., 2014; Wang and Bradburd, 2014).

The empirical challenge of clarifying patterns of adaptive divergence in wild populations can be addressed using “landscape genetics,” a multi-disciplinary study that relies on explicit statistical quantification of environmental and geographic variables to test for the relative influence of each on genetic variation across the landscape (Storfer et al., 2007). More recently, the arrival of “landscape genomics,” which utilizes genome-wide information to identify environmental features influencing genetic variation, has considerably increased the power to study local adaptation and adaptive divergence in wild populations—the initial outcomes of the process of divergent natural selection (Wang and Bradburd, 2014; Rellstab et al., 2015; Grummer et al., 2019). This includes assessment of a population’s ability to evolve and maintain fitness in response to environmental changes, or their “adaptive capacity,” an important issue in conservation management and applied evolution (Grummer et al., 2019). This is particularly timely as habitats across the globe are rapidly shifting in environmental properties, with many niches under notable threat (Thurman et al., 2020).

Limited research has explained how natural selection acts in aquatic environments (Kelley et al., 2016; Grummer et al., 2019). In rivers and streams, a constrained, unidirectional spatial flow of genes along well-defined migratory paths provides an ideal environment for the study of natural selection, as the impacts of geographic separation on genetic divergence can be controlled for (Hughes et al., 2009; Brauer et al., 2018; Attard et al., 2022). The Amazon, also known as Amazonia, holds the highest diversity of freshwater fishes in the world, yet little is understood about their evolution (Corlett and Primack, 2010; Beheregaray et al., 2014). Patterns of diversification in the Amazon have been mostly studied in the context of geographic isolation based on neutral genetic markers (Albert et al., 2011, 2020). Very few studies have assessed the role of natural selection in population divergence and evolution in the Amazon (reviewed in Beheregaray et al., 2014), an endeavor that requires information about adaptive genetic variation (Grummer et al., 2019).

The Amazon River network is composed of three distinctive river types with characteristic attributes resulting from the geological origin of water they contain, groundwater level, soil types, and vegetation (Val, 1995). Each water color shows a unique composition of pH, turbidity, water temperature and dissolved oxygen; thus “water color” can be used as an umbrella term to summarize differences in hydrochemical properties of water bodies. Variances in these hydrochemical parameters can drive constraints on the ecological communities inhabiting water bodies, and as such water color can be considered a major “aquatic ecotone” (Beheregaray et al., 2014). Water colors and their boundaries in the Amazon have been demonstrated to influence the distribution of forest floodplain avifauna (Laranjeiras et al., 2021) and the migratory paths of important fish species (McClain and Naiman, 2008).

From a landscape genetics perspective, several studies have investigated how genetic divergence in Amazonian fish populations might be influenced by water color (Cooke et al., 2012a,b,c, 2014; Beheregaray et al., 2014). This included the proposal of a framework capable of distinguishing the relative contributions of natural selection associated with aquatic ecotones vs. biogeographic history in the origin of fish population diversity (Beheregaray et al., 2014). These studies generally found support for water color as a selecting agent in Amazonian fish divergence. However, the limited number and the anonymous nature of the DNA markers used (amplified fragment length polymorphisms, AFLPs) offer little power to disentangle adaptive from neutral genetic variation (Luikart et al., 2003) and to gain insights about the functional relevance of putative adaptations (Grummer et al., 2019). Advancements in genome-wide screening technologies using single nucleotide polymorphisms (SNPs) has enabled a higher genotyping efficiency, data quality, analytical simplicity, and coverage across the genome (Luikart et al., 2003; Morin et al., 2004, 2009). This has allowed for greater power to differentiate adaptive from neutral genetic variation, enabling the identification of candidate loci influenced by natural selection (Vitalis et al., 2001; Morin et al., 2004, 2009; Seehausen et al., 2014).

This study focuses on the Amazonian benthopelagic fish Triportheus albus (Characiformes: Triporthidae) (Cope, 1871). Triportheus albus is a non-migratory species that has been recorded to grow up to 26 cm with a weight of up to 200 g (Giarrizzo et al., 2015). Little is known about the population dynamics and ecology of T. albus, a deficiency that applies to most other fish species in the Amazon (Prestes et al., 2010). This work aims to assess the relative contribution of spatial factors and environmental variation driving evolutionary divergence of T. albus, within the context of riverine “water color.” We hypothesize that, after controlling for spatial covariables, evolutionary divergence due to divergent selection should be stronger between populations of T. albus found in different ecotones than between populations from the same ecotone. This possibility derives from the strong physiological constraints that major hydrochemical and environmental gradients impose upon aquatic communities in Amazonia (Junk et al., 1983; Saint-Paul et al., 2000). To test this hypothesis, we generated a genome-wide dataset for T. albus populations that was integrated analytically with environmental and spatial datasets within a landscape genomics framework. We employed model-based and model-free approaches to assess neutral population differentiation and several Genotype-Environment Association (GEA) analyses to test how genetic variation might be associated with environmental predictors embodying different water colors. We also identify gene regions involved with adaptation to different water habitats and predict their role in driving natural selection. To the best of our knowledge, this is the first landscape genomics study for an Amazonian fish, providing valuable insights into how spatial and environmental factors simultaneously influence adaptation and population-level divergence within the ecologically rich Amazon Basin.

Materials and Methods

Sampling and Study Site

Our fieldwork explored over 1,100 km of total riverine distance in the Amazon Basin in February 2005 and in February 2008. Sampling was done at 10 sites representing five major river systems. These systems cover all three major hydrochemical environments of the Amazon Basin; the Amazon and the Madeira (white), the Negro (black), the Branco (seasonally black), and the Tapajós (clear) rivers (Figure 1). The white water Amazon River (known as Solimões upstream of Manaus) is central to the sampling transect and drains west to east, with the black water Negro, white water Madeira and clear water Tapajós rivers flowing into it. The black water portion of the Branco River that flows into the Negro was also included in the transect. Our sampling design consists of two ecological gradients, where black (Negro, Figure 1i) and clear (Tapajós, Figure 1ii) water meets white (Amazon) water, and two controls, where rivers of the same water color meet (Branco to Negro and Madeira to Amazon). Sampling was carried out in sandy and shallow beaches along the river banks and floodplains. This is the preferred habitat for T. albus (Figure 1iii), a species that has not been recorded in deeper river channels. Fish were caught with seine nets, euthanized, and muscle tissue preserved in 95% ethanol. Geographic coordinates of each site were obtained using a global positioning system (GPS). Measurements of pH, temperature (°C), turbidity (cm), dissolved oxygen (mg/L), and oxygen saturation (%) were collected in situ at each locality to assess environmental heterogeneity (Table 1).

FIGURE 1

Figure 1. Sampling localities of Triportheus albus in the Amazon Basin. Shape of site points correspond to water color of the locality (square = black water, triangle = clear water, circle = white water). Black water sites are labeled in blue points (B1 = Branco, N1 = Negro), clear water in purple (T1 = Tapajós), the white water Amazon in a gradient of yellow to orange (A1–A6), and the white water Madeira in pink (M1). Inset (i) shows the confluence of the Amazon and Negro rivers, and (ii) shows the confluence of the Amazon and Tapajós rivers. Inset (iii) depicts the study species T. albus. Inset (iv) displays the study area in northern South America.

TABLE 1

Table 1. Sampling locations of Triportheus albus in different water “colors.”

Laboratory Procedures

A modified salting-out protocol (Sunnucks and Hales, 1996) was used to extract DNA from muscle tissue samples. Double digest restriction-site associated DNA (ddRAD) libraries were developed using a modified protocol by Peterson et al. (2012), as detailed in Brauer et al. (2016). Libraries were sequenced at Novogene across two lanes of Illumina HiSeq4000 at 150 base-pair paired-end reads. Details of the DNA extraction method and the ddRAD library protocol are in Supplementary Appendix 1.

Bioinformatics and Single Nucleotide Polymorphism Filtering

Illumina HiSeq4000 sequencing produced a total of 902,555,657 sequence reads, each of around 150 bp. Quality control tests were initially performed using FASTQC 0.11.8 (Andrews, 2010). Sequences were demultiplexed and barcodes and adapters were removed using process_ragtags 2.4 from STACKS (Catchen et al., 2013). Low quality bases (Q < 20) and reads (more than 5% Ns, less than 40 bp, or < 30 average Q) were removed using TRIMMOMATIC 0.39 (Bolger et al., 2014), SNP variants were characterized using the DDOCENT 2.7.8 pipeline (Puritz et al., 2014), as detailed in Sandoval-Castillo et al. (2018). The SNPs were filtered in a pipeline of steps from VCFTOOLS (Danecek et al., 2011). A detailed summary of the filtering process is found in Supplementary Table 1.

Genetic Diversity and Population Structure

Percentage polymorphic loci, expected heterozygosity (He) and observed heterozygosity (Ho) were calculated across the filtered SNP dataset for each sampling locality using ARLEQUIN 3.5.2.2 (Excoffier and Lischer, 2010).

Population structure was assessed using model-based and model-free methods. Pairwise FST among sampling localities was estimated with the Jukes and Cantor (1969) method using 1,000 permutations. P-values were adjusted for a false discovery rate (FDR) of 5%, using the false discovery script.¹ A principal component analysis (PCA) of the total SNP dataset was generated using the program ade4 (Dray and Dufour, 2007) in R. Population structure was also inferred using Admixture (Alexander et al., 2009). The number of ancestral populations (K) was assessed by comparing log-likelihood ratios for multiple independent runs of each K (K = 1–7) before using a cross-validation (CV) procedure with 10 replicates to identify the K-value that best explains the number of random mating populations (Alexander et al., 2009). Admixture coefficients were plotted using the ggplot2 package (Wickham, 2009) in R. Relationships among individual fish were then assessed based on the SNPs with a Neighbor-joining (NJ) tree in PAUP4 (Wilgenbusch and Swofford, 2003) and the GRT + G model (Tavaré, 1986), which best fitted our data according to the Bayesian information criterion implemented in ModelTest-NG (Darriba et al., 2020). The aim of this analysis was not to infer phylogenetic history in T. albus, but to instead obtain a visual representation of the relationships among individuals across localities that can be contrasted with results of other approaches used to infer population structure.

Controlling for Spatial Genetic Structure

Patterns of genetic structure are known to be influenced by the physical structure of river catchments, including the spatial separation between populations via IBD, and the hierarchical structure of the dendritic system (Excoffier et al., 2009; Fourcade et al., 2013; Brauer et al., 2018). Riverine distance between sampling localities was estimated using Google Earth Pro 7.3.3 (2020) to test if any inferred genetic differentiation is consistent with the IBD model. A global Mantel test based on 999 permutations was performed to assess correlation between riverine distance and linearized genetic distance between all localities (FST/1—FST) using the ade4 package in R (Dray and Dufour, 2007). A second Mantel test (999 permutations) exclusively assessed the relationship between riverine distance and linearized genetic distance of white water Amazon and Madeira River localities. Localities A5 and A6 were omitted from this analysis due to insufficient sample size. These relationships were visualized using regression plots. A pairwise matrix of riverine distances between localities was calculated and translated to a set of synthetic IBD coordinates for each locality in R using the isoMDS function in the MASS library (Venables and Ripley, 2002) for use as a control for IBD in subsequent GEA analyses.

Hierarchical structure can be identified through distinct breaks in allele frequencies between clusters of populations, which may coincide with geographical features such as river catchments (Meirmans, 2012). A matrix of population covariance across locality allele frequencies (Omega matrix) was calculated through the core model in the BayPass program [explained below under “Genotype-Environment Association (GEA) Analyses; Gautier, 2015]. This matrix infers gene flow resulting from elements of demographic history, including neutral models of hierarchical structure. This matrix was translated to a set of covariance variables using the MASS library (Venables and Ripley, 2002) in R. The Omega coordinates were used as a control for hierarchical structure in successive GEA analyses.

Genotype-Environment Association Analyses

Five hydrochemical variables sampled in situ were considered in the GEA analyses; pH, temperature (°C), turbidity (cm), dissolved oxygen (mg/L), and oxygen saturation (%) (Table 1). These variables have been suggested as possible sources of ecological selection in fishes between Amazonian water colors (Beheregaray et al., 2014; Borghezan et al., 2021). The association between genetic and environmental heterogeneity was assessed using two multivariate approaches, redundancy analysis (RDA) (Legendre and Legendre, 1998) and Gradient Forest (Ellis et al., 2012), as well as by the univariate approach in BayPass (Gautier, 2015).

BayPass is a program for Bayesian population association analyses (Gautier, 2015). This Bayesian hierarchical model was proposed by Coop et al. (2010) as an FST-based approach to evaluate associations between variation in ecological variables and genetic markers, and to identify candidate loci for adaptive divergence while controlling for neutral covariance of alleles across populations. BayPass was used to investigate the signal of neutral and adaptive variation through three models; the core model, the standard covariate model, and the auxiliary model (Günther and Coop, 2013). Explanations of these three models are in Supplementary Appendix 2. The auxiliary model was used to identify candidate adaptive loci significantly correlated with each environmental variable over a log10(BF) threshold of 30.

The RDA is an asymmetric ordination method used to explore the relationship between multivariate response data and a set of explanatory variables using multiple linear regression (Legendre and Legendre, 1998). The RDA has been shown as a robust method to identify polygenic adaptive divergence through a strong ability to detect minor changes in allele frequency over many covarying loci (Grummer et al., 2019), while providing an optimal balance of low false-positive and high true-positive rates (Forester et al., 2018). An RDA was carried out using the vegan package in R studio (Dixon, 2003) comparing the response allele frequencies for each SNP with the five environmental variables (Table 1). To avoid collinearity, highly correlated variables were excluded using a variance inflation factor ≥10 (Dyer et al., 2010). First, we performed a standard RDA comparing variance in allele frequencies with variance in environmental variables. Subsequently, two partial RDAs were performed accounting for different aspects of spatial neutral population structure. These included synthetic IBD coordinates as a conditional variable and the synthetic Omega coordinates (see “Controlling for spatial genetic structure” section above) to control for other elements of demographic history, including hierarchical structure. A hierarchical analysis of molecular variance (AMOVA) and marginal AMOVA were performed to determine the significance of the partial and non-partial RDA models under a 0.05 threshold at 1,000 permutations in the R package vegan (Dixon, 2003). The SNPs exhibiting high contribution to environmental associations in the RDA (p = 0.05) were identified as candidate loci for adaptive divergence.

Random Forest (Breiman, 2001) is a machine-learning regression tree technique that can effectively handle thousands of DNA markers simultaneously to identify regions of the genome accounting for complex polygenic traits (Brieuc et al., 2018; Grummer et al., 2019). The original randomForest package was developed by Liaw and Wiener (2002). This was later modified by Ellis et al. (2012) into the gradientForest package, which extends the random forest method to the community level. Gradient Forest can pinpoint where compositional turnover occurs along an environmental gradient, which allows identification of important environmental thresholds that correlate to distinct changes in allele frequencies (Ellis et al., 2012). Gradient Forest was run to determine the correlation of allele frequencies with the predictor variables of environmental variation, synthetic IBD coordinates, and synthetic Omega coordinates (ntree = 750, mtry = number of variables/3, corr. threshold = 0.5). The cumulative importance of each predictor variable in shaping genetic populations was also assessed, with threshold values of each predictor in determining allele frequency breaks identified. Candidate adaptive loci were identified by assessing the R² weighted importance distribution of SNPs (Supplementary Figure 7). Loci above the upper elbow of this distribution curve (>0.86) were selected as candidates. As there is not a specific method implemented in random forest analyses to extract outlier SNPs (Goldstein et al., 2011; Laporte et al., 2016), the upper elbow approach used above has proven valuable in genome-wide association studies (Batley et al., 2019).

Functional Annotation

Candidate adaptive SNPs identified through RDA, Gradient Forest and BayPass methods were compared to select loci for functional analysis. The flanking sequence of each candidate SNP was aligned to the UniProt² Actinopterygii protein database using BLASTX 2.11.0 + (Altschul et al., 1997) with e-value threshold set to 1 × 10^––6. Gene ontology (GO) terms were assigned with a modified version of the python script GAWN.³ The functional roles of annotated candidate loci were investigated further by exploring the relevant literature.

Results

Sequencing Quality and Genetic Diversity

The SNP calling resulted in a total dataset of 743,123 raw SNP variants. After the filtering process (detailed in Supplementary Table 1), 15,251 high-quality SNPs were retained for analysis in 98 T. albus individuals (including four replicates). Genotyping error rates for the dataset were < 1%, with ∼0.05% missing data.

Adjacent Amazon River (Amazonas) localities with small sample sizes (A4, A5, and A6) were merged into a single sample to estimate genetic diversity. Genetic diversity was high and similar among localities, with average expected heterozygosity of 0.263 (He; 0.252–0.300), average observed heterozygosity of 0.240 (Ho; 0.221–0.275), and average polymorphic loci of 64.9% (54.4–76.1%) (Supplementary Table 2).

Population Structure

Population structure was strong across the region sampled in the Amazon Basin. In general, genetic differentiation among localities was highest between white (A1–A6, M1) and black (B1, N1) localities, and between white and clear (T1) localities. Genetic differentiation was nil to low between localities within the white and black selective environment (Supplementary Table 4). Pairwise FST estimates were significant (p < 0.05) after 5% FDR correction between black water localities, between black and clear localities, and between black or clear and white water localities. Pairwise FST estimates between white water localities were not significant (p < 0.05) after 5% FDR correction (Supplementary Table 4).

For the PCA and subsequent GEA analyses, Amazon River localities A5 and A6 were merged to a single sample (n = 4), with A4 (n = 6) kept as a separate sample. The PCA plot (Figure 2i) supported the pattern of population differentiation estimated via FST, showing distinct clusters of individuals among white, black, and clear water. A stronger level of separation was evident between white water T. albus and individuals from black and clear water localities, as explained by PC1 (17.412%). Clustering was also observable among black and clear water localities, explained by PC2 (2.871%)(Figure 2i). Two putative hybrids between white and clear water clusters were observed, originating from sites T1 and A5/A6. Possible migrants were also detected between black and white water, and between the clear water Tapajós and white water cluster (Figure 2i).

FIGURE 2

Figure 2. Population genomic structure of Triportheus albus based on 15,251 SNPs. Site labels are described as B1, Rio Branco (black water); N1, Rio Negro (black water); A1–A6, Rio Amazonas (white water); M1, Rio Madeira (white water); T1, Rio Tapajós (clear water). (i) PCA plot with each point indicating an individual and color a sampled site (shades of blue = black water sites, shades of orange = white water sites, purple = clear water site); (ii) admixture plot with each vertical bar depicting an individual and color representing ancestry from each water color catchment (blue = black water, orange = white water, purple = clear water). (iii) Neighbor-joining tree generated with PAUP4 based on the GRT + G substitution model. Catchment water color is labeled.

The Admixture analysis provided highest support for two genetic populations of T. albus, closely followed by three populations, as indicated by CV errors displayed in Supplementary Figure 1. The K = 2 plot combined the geographically distant black and the clear water samples into a single genetic population, with the second population represented by white water samples (Supplementary Figure 2). The K = 3 plot divided the black and clear catchments into two separate populations (Figure 2ii). An individual caught at A4 had ancestry to the black water population, supporting its migrant status suggested by the PCA. Two individuals from sites A5 and A6 shared ancestry between white and clear water populations, suggesting the possibility of an admixture zone (Figure 2ii).

The topology of the NJ tree was in agreement with results of population structure analyses (Figure 2iii). Individuals from the clear water Tapajós clustered as a single clade, whereas individuals from black and white water sites also comprised separate clades. Three individuals deviated from this pattern; one from white water A4, and one from each of the black water sites B1 and N1 (Figure 2iii). The outlier in site A4 corresponded to the potential migrant identified in the PCA and Admixture, supporting its suggested ancestry from the black water tributary. Potential hybrids from sites A5 and A6 identified in PCA and Admixture population structure analyses grouped closely to the clear water cluster in the NJ tree (Figure 2iii).

Genotype-Environment Associations and Controls for Models of Neutral Structure

The composition of environmental variables varied substantially between rivers of different water “color” (Table 1), and as such these are considered here as predictors of putative adaptive divergence within T. albus. The pH was considerably lower in the Negro River (5.2) than all other localities (6.7–7.2), including the seasonally black Branco River (6.9) (Table 1). Turbidity was substantially higher in the white water Amazon (14.5 cm) and Madeira (5.5 cm) sites compared to both black (79.7 cm) and clear (118 cm) water sites. Temperature remained relatively consistent across all systems. Dissolved oxygen was lowest at Amazon 3 (4.8 mg/L) and highest in the clear water Tapajós (7.1 mg/L). Oxygen saturation was greatest in the Tapajós (97.5%) and lowest at Amazon 5 (82.0%).

The Mantel test performed globally across white, black and clear water localities was not significant (p = 0.14), suggesting a relatively low effect of IBD at the broader regional scale. This was reflected in the global matrix regression plot (Supplementary Figure 3i) which shows no distinct pattern between riverine distance and FST among T. albus populations. The Mantel test among white water localities was also not significant (p = 0.09), although a trend of increasing genetic distance with geographic distance was observable in the regression plot (Supplementary Figure 3ii).

After controlling for highly correlated environmental predictors, pH and turbidity were included in the RDA model. The standard RDA separated individuals from white water localities from those in black and clear water sites (Figure 3i). This RDA model indicated that pH and turbidity accounted for 48.2% of total genetic variation and was globally significant (p = 0.010) (Supplementary Table 6). After controlling for IBD, which explained 37.7% of variation, pH and turbidity accounted for 32.0%, with the RDA model being globally significant (p = 0.036) (Supplementary Table 6 and Supplementary Figure 4). A clear pattern was less apparent in the RDA plot controlled by Omega (Figure 3ii). In that model, turbidity and pH explained 15.9% of the total genetic variation after accounting for demographic history (including hierarchical structure), which explained 53.9% of variation, with the RDA model not being significant (p = 0.384). From the RDA model conditional to Omega, we identified 1,749 candidate adaptive loci which associated highly with environmental predictors.

FIGURE 3

Figure 3. Genotype-environment associations of T. albus in the Amazon Basin based on RDA. The RDA shows how genetic differentiation between localities is explained by non-correlated environmental variables pH and turbidity (turb) collected in situ at each site. (i) Unconstrained, without controlling for spatial structure, and (ii) constrained to population covariance calculated with BayPass (Gautier, 2015) to control for neutral hierarchical genetic structure. Each colored point represents a sampling locality. Black water Branco River (B1) and Negro River (N1) localities = shades of blue, white water Amazon River localities = gradient of yellow to orange (A1–A5), white water Madeira River (M1) = pink, clear water Tapajós River (T1) = purple. Based on 15,251 filtered SNP loci.

The Gradient Forest identified that Omega 2, pH, Omega 1, and turbidity all correlated significantly with allele frequency at a similar level to one another (Figure 4i). The synthetic IBD coordinates, dissolved oxygen, temperature, and oxygen saturation did not correlate highly with allele frequency (Figure 4i). Demographic history predictors Omega 1 and 2 were the strongest in influencing breaks in allele frequencies of T. albus across Amazonia (Supplementary Figure 5). This was closely followed by turbidity and pH, which both influenced compositional turnover of allele frequency to a high level (Supplementary Figure 5). Turbidity classified one distinct compositional turnover in allele frequency, defined by ∼42 cm (Figure 4ii), which separated white water localities from black and clear water localities (Table 1). The pH defined two allele frequency breaks occurring at threshold pH values of 6.2 and 7.0 (Figure 4iii), which generally corresponded to pH differences between each of the three water colors (Table 1). The Branco River (pH 6.9; Table 1) was an outlier to this pattern, grouping with the clear water population. Dissolved oxygen, oxygen saturation and temperature did not substantially define any turnover in allele frequency (Supplementary Figure 6). The Gradient Forest identified 53 candidate adaptive loci associated with environmental variables over threshold R²_C value 0.86.

FIGURE 4

Figure 4. Gradient Forest (Ellis et al., 2012) results for associations between environmental variables and genetic diversity of Triportheus albus across the Amazon Basin. Based on 15,251 genome-wide SNPs. R² weighted importance (i) clarifies the level of correlation between each predictor variable and allele frequency variation. Omega 1 and Omega 2 denote control variables of neutral hierarchical structure, measured through a matrix of population covariance translated to two synthetic coordinates. The IBD1 and IBD2 are controls of isolation by distance, measured using a matrix of geographic distance between localities. Turbidity (turb), pH, dissolved oxygen (OD), oxygen saturation (Osat), and temperature (temp) are hydrochemical variables measured in situ at each sampling locality of T. albus. Cumulative importance displays the allele frequency compositional turnover functions for important environmental predictor variables of Triportheus albus divergence across the Amazon Basin, turbidity (ii) and pH (iii).

The auxiliary model in BayPass (Gautier, 2015) identified 268 candidate adaptive loci, mostly associated with pH (67) (Supplementary Figure 8). Of the 1,749 candidate loci identified by RDA and 268 by BayPass, 53 loci were consistent between the two methods (Supplementary Table 1). The 53 candidate SNPs identified by Gradient Forest were not discovered by RDA or BayPass as outliers (Supplementary Table 1).

Functional Annotation

For the 53 unique candidate loci identified by both BayPass and RDA, and the 53 loci identified by Gradient Forest, 29 were annotated and assigned to 58 GO terms. Many genes annotated were components of cellular membranes, with the most common GO terms “integral component of the membrane” and “plasma membrane” (10 and four of 29 annotated genes, respectively). Twelve of these 14 genes had molecular functions associated with the binding of ions, nucleic acids, or ATP. The most common biological processes identified were related to signaling pathways (four of 29 annotated genes) (Supplementary Table 7). Most candidates were highly associated with pH (14 of 23) or turbidity (7 of 23) (Supplementary Table 7). Further information about GO Terms assigned to genes, environmental variables that genes were associated to, and regression coefficients of the correlation between genetic variation and variation in the associated environmental variables are reported in Supplementary Table 7.

Discussion

Disentangling adaptive from neutral genetic diversity in wild populations is a difficult task, but the ease and power to do so has improved substantially with the development of genome-wide screening techniques and approaches in landscape genomics (Yeaman et al., 2016; Grummer et al., 2019). Here, we were able to identify the contributions of spatial and environmental factors in shaping genomic divergence in populations of the characin tetra T. albus from central Amazonia. Using a robust approach that integrated multivariate RDA and Gradient Forest, and univariate BayPass approaches, we have demonstrated that in addition to strong neutral hierarchical structure, adaptive divergence is likely to be occurring in T. albus populations across black, white and clear water selective environments.

Evolution and Population Genomic Structure

Our study is among the first aimed at assessing population genomic structure of an Amazonian aquatic species using genome-wide approaches. To the best of our knowledge, a study by Torati et al. (2019) examining stock structure in the teleost Arapaima gigas between South American basins using 393 SNPs represents the only prior attempt to this. Landscape genetic methods have been used to identify population structure among Amazonian catchments for the characin T. albus (Cooke et al., 2012a), the silver croaker Plagioscion squamosissimus (Cooke et al., 2012b), the riverine puffer Colomesus asellus (Cooke et al., 2012c), and the barred knife fish Steatogenys elegans (Cooke et al., 2014). Genetic structure between Amazon and Madeira River populations has also been identified for the catfish Brachyplatystoma platynemum (Ochoa et al., 2015). None of these studies have employed datasets and statistical approaches capable of inferring genome-wide neutral and adaptive signal (Luikart et al., 2003; Grummer et al., 2019). As such, comparisons of genome-wide diversity with our T. albus dataset are hampered by the remarkably small number of aquatic studies (including for characiform tetras from other basins). For instance, for the least-concern (IUCN) characiform Astyanax mexicanus (blind cavefish) in Northeast Mexico (Bradic et al., 2013), population-level observed heterozygosity ranged from 0.05 to 0.19. For the abundant atheriniform Melanotaenia fluviatilis (Murray River rainbowfish) in Australia, heterozygosity values ranged from 0.111 to 0.317 (Brauer et al., 2018). For our study, genetic diversity was not only moderate to high but was markedly similar across all sampled sites (observed heterozygosity ranged from 0.221 to 0.275), suggesting high adaptive potential for the inferred populations in our study region.

When assessing evolutionary patterns of freshwater biodiversity it is important to consider the history and timing of geological formation of river basins. The Characiformes originated over 100 million years ago (Ma) in Gondwana (Ortí and Meyer, 1997), and as such deep evolutionary divergence in the Neotropics should be assessed in relation to the uplift of the Andes and the initial development of the Amazon drainage basin (Rull, 2008; Hoorn et al., 2010a). Prior to the influence of the Andes, the Amazon was an undivided sedimentary basin dominated by the fluvial system of the cratonic shields of which the Negro and Tapajós rivers are derived from Figueiredo et al. (2009) and Hoorn et al. (2010b). Increased sedimentation, rising sea levels, and overfilling from the late Miocene (∼6.4 Ma) led to the establishment of west to east transcontinental flow and formation of the modern Andean-dominated Amazon River network by the late Pliocene (∼2.5 Ma) (Val, 1995; Campbell et al., 2006; Figueiredo et al., 2009).

Geomorphological evolution is thought to have impacted the phylogeography of T. albus. The work by Cooke et al. (2012a) based on 360 AFLP loci supported two historically distinct populations of T. albus in central Amazon; a white water Amazonas/Madeira population, and a black water Negro and clear water Tapajós population. Phylogenetic and phylogeographic analyses indicated a recent divergence between white water and black/clear water ecotypes of T. albus followed by demographic expansions dated to the late Pleistocene (Cooke et al., 2012a). These findings, paired with geomorphological and paleoecological evidence for the basin, indicate that the black/clear water ecotype is older, with the west to east transcontinental formation of the Amazon River promoting the colonization and demographic expansion of the white water T. albus ecotype (Lundberg et al., 1998; Campbell et al., 2006; Cooke et al., 2012a). This scenario is supported by comparative phylogeographic studies of co-distributed fishes; Colomesus asellus, Plagioscion squamosissimus, and Steatogenys elegans (Cooke et al., 2012b,c, 2014). Although assessing phylogeographic history of T. albus is outside the scope of this study, our results support the strong divergence between white and black/clear ecotypes, as well as the similarity between black and clear water populations detected by Cooke et al. (2012a). The latter was found despite the large contemporary riverine distance between the black water Negro and the clear water Tapajós (up to 1,100 km) and their different geological settings (Latrubesse and Franzinelli, 2005; Hoorn et al., 2010b), pointing to the strong influence of phylogeographic history (Cooke et al., 2012a) on contemporary patterns of population structure in T. albus.

In addition, the greater power offered by our genome-wide dataset revealed weaker, but biologically meaningful, genomic divergence between black and clear water populations. Only a low level of admixture was noted between black and white, and between clear and white populations, suggesting downstream migration following the direction of river flow; from tributaries into the main Amazon channel (Figure 2ii). This is consistent with the isolation with migration model results of Cooke et al. (2012a), which inferred a higher probability of unidirectional gene flow in T. albus from black and clear water into white water than of any other alternative scenario. Overall, population structure analyses indicated three contemporary genetic clusters for T. albus in central Amazonia that correspond to rivers of distinct water color. As genetic structure was not observed between the white water Amazon and Madeira rivers, divergence due to the confluence of a tributary alone remains an unsupported explanation. Similarly, a simple model of IBD seems insufficient to explain divergence patterns across the vast study region. Instead, our study results suggest that population divergence is also associated with the boundaries of catchments that harbor divergent hydrochemical conditions.

Neutral and Adaptive Contributions to Genomic Divergence

Species that inhabit contrasting selective environments provide a valuable opportunity to study divergent natural selection (Schluter, 2000). The dendritic arrangement of the Amazon Basin encompasses distinct ecotones of varying environmental conditions characterizing different water “colors,” which may theoretically promote divergent selection and IBE across populations of aquatic organisms. For instance, Pires et al. (2018) described divergence and reproductive isolation between Amazon and Negro River lineages of the sailfin tetra Crenuchus spilurus. A landscape genetics study by Cooke et al. (2012a) suggested a putative pattern of IBE for T. albus across Amazonia, as well as several co-inhabiting fishes (Cooke et al., 2012b,c, 2014; Beheregaray et al., 2014). In those studies, evidence for putative divergent selection lay in patterns of diversification that aligned best with catchments of similar water properties, irrespective of riverine distance. The genomic dataset used in this study not only provides comparatively higher resolution than the above genetic datasets to infer fine- and broad-scale population divergences (Luikart et al., 2003), but also contains sequence information around candidate adaptive regions that might prove useful to identify mechanisms influenced by divergent selection.

The IBD model is not uncommon in freshwater systems and has been observed for species with varying dispersal potential (e.g., Hrbek et al., 2005; Hubert et al., 2007; Zieritz et al., 2010; Crookes and Shaw, 2016). The GEA analyses strongly indicated that genetic divergence among T. albus populations cannot be explained by a simple IBD model. Broadly speaking, we found no clear relationship between genetic and riverine distance across the entire study area in central Amazon, which was also reflected in both the RDA and Gradient Forest results. A strong pattern of genetic structure was found between white water and black/clear water localities, but not between white water Amazon and Madeira Rivers. Similarly, very low differentiation and a low but non-significant signal of IBD was present along the Amazon River, despite the ∼750 km separating sites A1 and A6. Critically, the low levels of neutral differentiation and corresponding reduced genetic drift along white waters are not expected to impact on our ability to characterize putative adaptive divergence within that selective environment (Grummer et al., 2019).

The stream hierarchy model (Meffe and Vrijenhoek, 1988) should be routinely considered when assessing gene flow among populations in complex dendritic spatial environments (Brauer et al., 2018). Consistent with the broad inferences of population structure, the RDA analysis demonstrated hierarchical structure as the prevailing factor driving genetic divergence among populations (Figure 3ii). Although pH and turbidity explained a substantial proportion of the variance in all RDAs, their contributions were not statistically significant as model complexity increased following the addition of demographic and spatial factors, such as hierarchical structure (Supplementary Table 6). This is expected due to the complex links between spatial-demographic factors and environmental variables and the resulting difficulties in separating adaptive from neutral variation (Grummer et al., 2019). As such, analyses did not suitably disentangle whether divergence followed a neutral hierarchical pattern resulting exclusively from the arrangement of tributaries, or a combination of neutral and non-neutral influences due to the varying environmental conditions within each of the catchments.

The Gradient Forest analysis also revealed neutral hierarchical structure as the greatest driver of genetic differentiation, however, turbidity and pH were revealed to correlate with allele frequency at comparably high levels (Figure 4i). Heterogeneity in turbidity created two distinct genetic populations of T. albus, separating the white water localities from the black and clear sites (Figure 4ii). The white Amazonian headwaters are subject to copious erosion in the Andes, and are consequently rich in suspended sediments (Sioli, 1984; Val, 1995; McClain and Naiman, 2008). Associated with this heavy sediment load is an abundance of nutrients and organic matter which facilitate a higher productivity and diversity of species in white waters than adjoining black and clear tributaries (McClain and Naiman, 2008; Hoorn et al., 2010a). However, characiformes are generally visually-orientated, surface-dwelling diurnal fishes (Tejerina-Garro et al., 1998), and consequently may be offered an ecological advantage in the transparent, mineral-deficient black and clear water tributaries (Wallace, 1854; Furch, 1984; McClain and Naiman, 2008). A detailed proposal of how Amazonian water types may promote differentiation in the sensory mechanisms of fish has been outlined by Borghezan et al. (2021).

Variation in pH across the basin is associated with the split T. albus into three distinct genetic populations corresponding to each water color (Figure 4iii and Table 1). The dark color of the Negro River is attributed to staining from tannins and humic acids leached by decaying vegetation (Wallace, 1854; McClain and Naiman, 2008), and as a result exhibits lower pH than their pH-neutral white and clear counterparts. The acidic, nutrient-poor properties of black water systems forms a challenging environment for aquatic species (Val, 1995), however, research has suggested that phenotypic plasticity of mechanisms conferring acid-base regulation has enabled T. albus to survive in harsh Negro waters (Araújo et al., 2017). The Branco River has white headwaters of mountainous origin, however, sampling took place near its junction with the Negro, thus this site was heavily influenced by cratonic black water. Even so, geochemical analyses have revealed the Branco River to be chemically and sedimentologically intermediate between black and white water (Küchler et al., 2000; Evangelista and Tosi, 2015). The clear water Tapajós represents a second intermediate physiochemical condition to the Amazon and Negro rivers (Duncan and Fernandes, 2010), hence why the pH of the Branco and Tapajós tributaries are similar. However, the influence of environmental discontinuities on population divergence presented by the Tapajós should be cautiously interpreted given the lack of replication in this river system.

A key limitation of the study was that sample size was small in terms of individuals for some white water sites, and in terms of locality replication for the Tapajós tributary. However, simulations by Gaughran et al. (2018) demonstrated that accurate estimates of genetic differentiation in highly structured populations can be obtained using thousands of SNPs and only 2–5 individuals per locality. Thus, merging sites A5 and A6 for analyses, should be sufficient to provide a reasonable representation of their genetic patterns, and is justified by their genetic similarity identified in population structure analyses. A greater number of sampled sites in the Tapajós would provide an improved understanding of the structure of the clear water tributary as a whole. The study was also limited as a whole genome reference is not available for T. albus; the reference of a closely related species would help to better characterize regions under selection. In addition, the hypothesis would be better tested by performing a comparative analysis on multiple co-distributed species, as done by Beheregaray et al. (2014) using a landscape genetics framework.

The discovery of 106 candidate loci highly associated with variation in the environment provides further support for adaptive divergent selection of T. albus across the Amazon Basin. Although no Gradient Forest candidates were identified as outliers in the RDA or BayPass analyses, it is important to consider that methods of outlier detection all vary in assumptions and balance between low false positive and high true positive rates (Narum and Hess, 2011). The strong associations of pH and turbidity with allele frequency that correspond to environmental differences between each water color, and the discovery of candidate adaptive genes for environmental associations, are supportive that characteristics of water color can perhaps induce selection and promote adaptive divergence between populations.

Functions of Candidate Adaptive Genes

Establishing the function of candidate adaptive genes and their ecological importance should be done cautiously, as a large fraction of genes remains without annotation of ecological relevance (Pavey et al., 2012). The correlation of potentially adaptive genes with mechanisms linked to mediation of the impacts of environmental fluctuations do not necessarily imply causation, but nevertheless can be useful to understand how adaptive divergence may arise (Grummer et al., 2019). Most candidates identified here were annotated to genes associated with ion channel activity (Supplementary Table 7), which is consistent with up-regulated GO terms for T. albus reported by Araújo et al. (2017). As aquatic environments are composed of a wide range of salinities, ion compositions, and pH values, fish are required to cope with challenging osmotic and ionic gradients (Hwang and Lee, 2007). Ion channels can be sensitive to changes in extracellular pH (Holzer, 2003), and fluxes of ions play an important role in maintaining acid-base homeostasis in freshwater fish (Goss et al., 1992). The black water Rio Negro drains the nutrient-poor soils of the Amazon forest, and consequently contains low concentrations of nutrients such as Na⁺, Cl^–, K⁺, Ca²⁺, and Mg²⁺ (Furch, 1984; Val, 1995). The gene coding for the K⁺ channel protein KCNK1 was identified in T. albus to be associated with environmental fluctuations across the study region, specifically in relation to pH in the RDA (Supplementary Table 7). The KCNK1 protein was found to be sensitive to low pH in zebrafish (Christensen et al., 2016), alluding that its function in T. albus may be sensitive to the strong fluctuations in pH across Amazonia.

Diurnal fishes such as the Characiformes depend highly on vision as a source of sensory information, and as a result have large, well-developed eyes (Guthrie, 1986; Tejerina-Garro et al., 1998). There is ample research demonstrating the effects of turbidity and light intensity on the feeding ability (Gardner, 1981; Rowe and Dean, 1998; Leahy et al., 2011), and antipredator behavior (Higham et al., 2015; Kimbell and Morrell, 2015) of fishes. Calcium ions and stores play an important role in the detection and transfer of light stimuli in photoreceptors (Krizaj and Copenhagen, 2002; Križaj, 2012). The ryanodine receptors (RYRs) are a family of Ca2+ release channels (Sutko and Airey, 1996) expressed in the photoreceptors of vertebrates (Križaj et al., 2004; Križaj, 2012). Photoreceptor pathways are sensitive to light (Križaj et al., 2011), and the RYRs have been linked to visual stimuli response in fish (Frank et al., 2019). The genes RYR2 and RYR3 were highly associated with environmental heterogeneity in T. albus, with RYR3 most highly correlated with variation in turbidity in the Gradient Forest analysis. Considering this, differences in light intensity produced by varying turbidity across Amazonia may be impacting the expression of the RYRs in T. albus, supporting sensory systems research directions suggested by Borghezan et al. (2021).

Fluctuations in dissolved oxygen is an issue in many aquatic systems including the Amazon (Junk et al., 1983) and can influence species distributions (Mandic et al., 2009) and cause negative impacts on freshwater fishes (Landman et al., 2005). Hypoxia is a condition in which cells suffer from oxygen deficiency (Hughes, 1973). Hypoxia has been found to affect the RYR2 gene in the large yellow croaker Larimichthys crocea (Mu et al., 2020). This gene was statistically associated to variation in oxygen saturation in T. albus (Supplementary Table 7). The connections of candidate genes to functions sensitive to environmental heterogeneity support the role of adaptive divergence in the evolution of T. albus across catchments of different water conditions.

Conclusion

By disentangling signal of adaptive and neutral divergence in T. albus across Amazonia using a novel landscape genomics framework, we have identified genetically distinct populations in black, white and clear water catchments of the Amazon Basin. The dominant driver of population divergence across the basin was hierarchical structure, which is expected in dendritic river systems (Brauer et al., 2018). However, pH and turbidity were identified as having a significant secondary influence on genetic patterns across the system, which was supported by the discovery of candidate adaptive genes with functions linked to mediating the biological impacts of hydrochemical fluctuations. Thus, pH and turbidity, which are of distinctive compositions in contrasting water colors, may be important agents of divergent natural selection. In spite of the known limitations of genome scans of selection (reviewed in Grummer et al., 2019), this study has improved our understanding of the evolutionary processes operating in a complex tropical system, which is vastly understudied in proportion to its immense biodiversity. The common IBD model appears as not being sufficient to explain broad patterns of genetic variation in this system. Importantly, the notable impact of IBE in freshwater fishes suggests that changing hydrological conditions of the Amazon Basin may have implications to their distribution and persistence. Understanding genomic vulnerability of Neotropical fishes amidst accelerated environmental and climatic change may be a valuable avenue for future research and would further help to inform on conservation management strategies for this diverse group of vertebrates.

Data Availability Statement

The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found at: Figshare (https://figshare.com/articles/dataset/_/17109146).

Ethics Statement

The animal study was reviewed and approved by the Macquarie University Animal Welfare Committee 2007/033.

Author Contributions

AH carried out the molecular laboratory work, performed the data analyses, and wrote the first manuscript draft. JS-C provided technical and analytical advice, supervised the laboratory work and data analyses, and made suggestions for the manuscript. GC carried out fieldwork and laboratory work. NC carried out fieldwork and provided guidance about project design. LB conceptualized and designed the project, obtained funding, carried out fieldwork, provided technical and analytical advice, critically reviewed, and edited the manuscript. All authors contributed to the article and approved the submitted version.

Funding

This project was funded by the Discovery Program of the Australian Research Council (ARC DP0556496 to LB). Logistics and local arrangements were also supported in part through the Brazilian National Council of Research and Technology CNPq- SEAP 408782/2006-4 via NC (formerly at the Universidade Federal do Amazonas), and by Macquarie University through a postgraduate travel grant and research award to GC. Collection permit is under IBAMA No. 1920550, and ethical approval was received from Macquarie University (2007/033).

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.2022.825406/full#supplementary-material

Footnotes

References

Albert, J. S., Carvalho, T. P., Petry, P., Holder, M. A., Maxime, E. L., Espino, J., et al. (2011). Aquatic biodiversity in the amazon: habitat specialization and geographic isolation promote species richness. Animals 1, 205–241. doi: 10.3390/ani1020205

PubMed Abstract | CrossRef Full Text | Google Scholar

Albert, J. S., Tagliacollo, V. A., and Dagosta, F. (2020). Diversification of neotropical freshwater fishes. Annu. Rev. Ecol. Evol. Syst. 51, 27–53. doi: 10.1146/annurev-ecolsys-011620-031032