Broad-scale sampling of primary freshwater fish populations reveals the role of intrinsic traits, inter-basin connectivity, drainage area and latitude on shaping contemporary patterns of genetic diversity

Target species

A total of 17 Iberian endemic cyprinids were selected as target species: Anaecypris hispanica (Steindachner, 1866), Achondrostoma oligolepis (Robalo, Doadrio, Almada & Kottelat, 2005), Achondrostoma occidentale (Robalo, Almada, Sousa-Santos, Moreira & Doadrio, 2005), Iberochondrostoma lemmingii (Steindachner, 1866), Iberochondrostoma lusitanicum (Collares-Pereira, 1980), Iberochondrostoma almacai (Coelho, Mesquita & Collares-Pereira, 2005), Pseudochondrostoma polylepis (Steindachner, 1865), Pseudochondrostoma duriense (Coelho, 1985), Pseudochondrostoma willkommii (Steindachner, 1866), Luciobarbus microcephalus (Almaça, 1967), Luciobarbus bocagei (Steindachner, 1865), Luciobarbus sclateri (Günther, 1868), Luciobarbus comizo (Steindachner, 1865), Squalius carolitertii (Doadrio, 1987), Squalius pyrenaicus (Günther, 1868), Squalius torgalensis (Bogutskaya, Rodrigues & Collares-Pereira, 1998) and Squalius aradensis (Bogutskaya, Rodrigues & Collares-Pereira, 1998). Of all the cyprinids native to Portugal, only four species [Squalius alburnoides (Steindachner, 1866), Luciobarbus steindachneri (Almaça, 1967), Achondrostoma arcasii (Steindachner, 1866) and Iberochondrostoma olisiponensis (Gante, Santos & Alves, 2007)] were not included a priori due to their hybridogenetic origin (Sousa-Santos, Collares-Pereira & Almada, 2007), uncertain taxonomic classification (Robalo et al., 2006; Gante et al., 2015) or extreme scarcity in wild populations (Sousa-Santos et al., 2014a). L. microcephalus was posteriorly excluded from the analyses due to the low number of individuals sampled in each population (see below).

Sampling

Thirty four river basins were sampled in the Portuguese mainland hydrographical network. The six largest river basins (Douro, Vouga, Mondego, Tagus, Sado and Guadiana) were further sub-divided into 47 sub-basins resulting in a total number of 81 geographical units sampled (Table S1 and Fig. 1A). Populations for which less than 15 individuals were collected were excluded from analyses.

Fish were collected with standard wadable electrofishing procedures (CEN, 2003) and returned to the water immediately after non-destructive sampling. Collected fin clips were preserved in 96% ethanol and vouchers were kept at the tissue collection of MARE/ISPA for subsequent DNA extraction, amplification and sequencing. Permits for field work were given by ICNF (permit number 176/2010/CAPT and 53/2012/CAPT).

DNA extraction, amplification and sequencing

Total genomic DNA was extracted from fin clips using REDExtract-N-Amp Tissue PCR kits (Sigma-Aldrich) following the manufacturer’s instructions. The mitochondrial cytochrome b (cytb) gene was amplified using the primers LCB1-new ACTTGAAGAACCACCGTTG (adapted from the LCB1 primer described by Brito et al., 1997) and HA-CAACGATCTCCGGTTTACAAGAC (Schmidt & Gold, 1993). PCR conditions were the following: 35 × (94°C1′ + 50°C1′ + 72°C2′). PCR products were purified and sequenced in the forward direction using the LCB1-new primer, at GATC Biotech (Konstanz, Germany). Obtained sequences were trimmed at the 3′ and 5′ ends so they had the same length for all the individuals sampled (720bp) and deposited in GenBank ( KU366823– KU370500).

DNA analyses and gene diversity mapping

All sequences were edited and aligned using CodonCode Aligner v4.0.4 (CodonCode Corp., USA). The search for shared haplotypes and the listing of representative sequences was conducted by DNAcollapser (FaBox v.1.41; http://www.birc.au.dk/software/fabox).

ARLEQUIN software package V.3.5 (Excoffier & Lischer, 2010) was used to quantify the number of private haplotypes for each population (i.e., haplotypes that were exclusive for the considered population and that were not found elsewhere). These values were then used to calculate the percentage of private haplotypes per population (of all the haplotypes found in the considered population) and the average percentage of private haplotypes per population for each target species (%N_PH). Mean, standard deviation, minimum and maximum values of %N_PH were obtained with EXCEL 2013 (Microsoft^®).

ARLEQUIN was also used to estimate gene diversity (h index, defined as the probability that two randomly chosen haplotypes are different in a sample; used as a measure of haplotype diversity for haploid data), nucleotide diversity (π index, defined as the probability that two randomly chosen homologous nucleotides are different) and mean number of pairwise differences (MNPD, defined as the mean number of differences between all pairs of haplotypes in a sample) for each population. Analyses of molecular variance (AMOVA) were also performed with ARLEQUIN.

Values obtained for the h index were mapped in the sampling areas by Yris Graphics (www.yrisgraphics.com). This index, which varies between 0 and 1, reflects the probability of two randomly chosen haplotypes being different in a sample, and was selected to illustrate haplotype diversity of sampled populations. Using this index, diversity classes (common for all species) were established, allowing for an automatic visual inspection of the diversity level of distinct populations. The maps and database are available for download at the project’s webpage (www.fishatlas.net).

Data analyses and hypothesis testing

Populations, regardless of species, were considered the unit of comparison. In order to test whether the obtained genetic diversity pattern was influenced by species’ intrinsic traits and extrinsic factors, data concerning 8 variables were collected and organized in a matrix (Table S2). Six of these variables are categorical: (1) “species” (with 16 categories, corresponding to the 16 species studied: L. bocagei, L. sclateri, L. comizo, P. duriense, P. polylepis, P. willkommii, I. almacai, I. lusitanicum, I. lemmingii, S. carolitertii, S. pyrenaicus, S. aradensis, S. torgalensis, A. hispanica, A. oligolepis and A. occidentale), and the dichotomous variables (2) “hydrological regime,” (3) “latitude,” (4) “migratory behaviour,” (5) “species range,” and (6) “inter-basin connectivity.” The quantitative variables “drainage area” and “species maximum size” were also included in the matrix.

“Hydrological regime” was classified as either permanent or temporary, depending on whether or not the population occurs in river basins that maintain flowing water throughout the year; “latitude” as northern or southern (for river basins located north or south of the Central Massif of Estrela, which divides the Portuguese territory in two halves of distinct reliefs and climate); “migratory behaviour” as non-migratory or potamodromous (using ecological data compiled by Ribeiro et al., 2007); “species range” as wide or restricted if the species occurs, respectively, in more or less than 12 river basins/sub-basins (see Sousa-Santos et al., 2013 for the distribution areas of each species); and “inter-basin connectivity” as connected or unconnected to other water bodies. Drainage areas were calculated from the shape file of Portuguese hydrographical network (downloaded from http://sniamb.apambiente.pt/Home/Default.htm) by using the polygon area as implemented in QGIS 2.10.1 software. Whenever a population inhabits a sub-basin of a larger basin, the area of the sub-drainage was considered, as for the studied small fish species the main course of a large river may represent a natural barrier to gene flow.

Individual linear regressions were performed to test the effect of each of these independent variables on h, π and MNPD indices calculated for each population. Dummy variables were created for the categorical independent variable “species,” enabling the use of the 16 categories of this variable in a regression model.

Predictor variables with significant effects (α < 0.05) were selected as candidate variables for the modelling process. Then, a hierarchical linear regression was performed, first including all selected variables (excluding “species”) with a stepwise method and, finally including the categorical variable “species.” This procedure sought to first analyse the effect of each of the variables regardless of the effect of “species.” This variable was then included to test whether other aspects intrinsic to the species, not measured by the remaining independent variables, were significant. This method was adopted independently for each of the three dependent variables: h, π and MNPD. We search for the presence of outliers and excluded them to avoid spurious trends or masking of valid ones.

For the analyses of each species separately, given the lower number of samples, non-parametric tests were performed: Spearman’s ρ, to test the correlation between “drainage area” and each of the dependent variables; and Mann–Whitney’s U tests to compare groups (populations inhabiting connected or unconnected water bodies; populations occurring in permanent or temporary river basins; and northern and southern populations) concerning their genetic diversity (h, π and MNPD indices). All statistical analyses were conducted using IBM SPSS Statistics, version 22 (IBM Corporation, 2013).

Patterns of genetic diversity

In general, populations of the small sized non-migratory species A. hispanica, A. oligolepis, I. lemmingii, I. lusitanicum, I. almacai, S. aradensis, S. carolitertii and S. pyrenaicus show higher percentages of private haplotypes (average values per species ranging from 27.84% to 75.00%, Table 1) than populations of the larger sized potamodromous species L. bocagei, L. comizo, L. sclateri, P. willkommii, P. duriense and P. polylepis (average values per species ranging from 2.14% to 25.00%, Table 1).

The average values of haplotype diversity ranged between h = 0.061 ± 0.029 (for L. comizo) and h = 0.936 ± 0.010 for (A. hispanica). The average nucleotide diversity and the average mean number of pairwise differences showed the same pattern (Table 1): L. comizo presented the lowest values (π = 0.000 ± 0.000 and MNPD = 0.061 ± 0.133, respectively) and S. pyrenaicus the highest (π = 0.014 ± 0.007 and MNPD = 10.080 ± 4.616, respectively).

At the population level, the spatial distribution of haplotype diversity values revealed distinct levels of diversity among populations of the same species and distinct patterns among species (Fig. 2). Indeed, analyses of molecular variance (AMOVAs) conducted for each target species separately revealed two contrasting patterns: in eight of the 15 species most of the variation (ranging from 50.65% for L. bocagei to 91.91% for A. oligolepis) could be attributed to differences among populations, while for the remaining seven species 53.26% (for P. polylepis) to 98.95% (for L. comizo) of the variation was explained by genetic differentiation within populations (Table 2).

Genetic diversity determinants

Preliminary linear regressions retrieved no outliers when using h as the dependent variable. Contrastingly, when using the other genetic diversity indices six outlier populations were detected: S. carolitertii from Mondego-Ceira (for π only), S. pyrenaicus from Tagus-Erges (for MNPD only), and S. pyrenaicus from Tagus-Grande da Pipa, S. pyrenaicus from Sado-São Martinho, I. lusitanicum from Sado-São Martinho and A. hispanica from Guadiana-Ardila (for both π and MNPD). These populations show a higher number of point mutations, probably due to an admixture of individuals from distinct haplogroups, and those differences between the sequences are reflected in the higher values of π and MNPD obtained.

Linear regressions between selected independent variables and h values showed that “species maximum size,” “inter-basin connectivity,” “species range,” “latitude,” “migratory behaviour” and “species” had a significant effect on the diversity of populations (Table 3). The independent variables “drainage area” and “hydrological regime” had no significant effect on haplotype diversity (Table 3). After the removal of the outliers, linear regressions using π and MNPD as dependent variables retrieved the same pattern (Table 3).

Regarding the effect on the h index, the inclusion of all independent variables (excluding “species”) in the linear regression, using a stepwise method, indicated that the best fit model only includes “species maximum size,” “inter-basin connectivity” and “latitude.” Non-significant predictor variables “range” and “migratory behaviour” were excluded from the analysis (Table 4). This model explains 26.7% of the haplotype diversity variance (Table 4). The inclusion of the variable “species” in a hierarchical regression analysis had a significant effect, increasing the explained variance to 32.3% (Table 4).

When applying the same stepwise procedure to the remaining two genetic diversity indices, the results from the linear regressions showed a similar pattern: the same variables were included in the model (“species maximum size,” “inter-basin connectivity,” “latitude” and “species”) and the best fit model explained 33.7% and 26.7% of the variance, respectively, for π and MNPD (Table 4).

Concerning the maximum size attained by the species, the results showed a negative correlation with all three genetic diversity indices (Table 3). Regarding “inter-basin connectivity,” the results indicated that populations occurring in sub-basins connected with other water bodies have higher genetic diversity (mean h = 0.413 ± 0.306, N = 135; mean π = 0.00121 ± 0.00120, N = 131; mean MNPD = 0.831 ± 0.856, N = 131) than those occurring in unconnected river basins (mean h = 0.252 ± 0.267, N = 53; mean π = 0.000606 ± 0.000766, N = 52; mean MNPD = 0.508 ± 0.614, N = 52). Finally, regarding “latitude,” southern populations exhibited higher levels of haplotype diversity and mean number of pairwise differences (mean h = 0.417 ± 0.292, N = 108; mean π = 0.00126 ± 0.0.00119, N = 104; mean MNPD = 0.871 ± 0.830, N = 103) than northern populations (mean h = 0.302 ± 0.308, N = 80; mean π = 0.000749 ± 0.000974, N = 79; mean MNPD = 0.570 ± 0.747, N = 80), a tendency which was already evident from the diversity mapping depicted for the National Genetic Atlas (Fig. 2).

It is worth mentioning that the variable “species” was significantly correlated with the three genetic diversity indices (Table 3) and was included in the model which best explains the observed variance for all of the genetic diversity indices (Table 4). The inclusion of this variable in the best fit models highlights the influence of species intrinsic idiosyncrasies other than those tested explicitly herein (migratory behaviour, maximum size and range). Since it is not possible, with the present dataset, to disentangle which intrinsic traits influence genetic diversity the most, we further analysed which of the extrinsic factors (“drainage area,” “inter-basin connectivity,” “latitude” and “hydrological regime”) were determinant for the genetic diversity pattern observed for each species. The results presented in Table 5 show that the environmental variables tested have no influence on the genetic diversity for most of the species, except for four species (L. bocagei, P. polylepis, S. pyrenaicus and I. lusitanicum).

More specifically, all of the genetic diversity indices were significantly correlated with “inter-basin connectivity” for I. lusitanicum and with “latitude” for L. bocagei and P. polylepis (Table 5). The haplotype diversity obtained was also significantly correlated with “inter-basin connectivity” for S. pyrenaicus and with “drainage area” for I. lusitanicum (Table 5).

Historical determinants of genetic diversity

The present study, targeting multiple co-distributed species of primary freshwater fish, assessed the relative role of historical versus contemporary factors affecting genetic diversity. Indeed, the origin of Iberian lineages of Anaecypris, ex-Chondrostoma, Luciobarbus and Squalius dates back to the Miocene, around 19–7.7 Mya (Levy, Doadrio & Almada, 2008; Gante, 2011). Speciation within these genera and subsequent diversification must have occurred through the same available connections between paleobasins until the Pleistocene-Holocene, when the current hydrographical network became established (Sousa-Santos, Collares-Pereira & Almada, 2007; Gante et al., 2009; Almada & Sousa-Santos, 2010; Sousa-Santos et al., 2014a; Sousa-Santos et al., 2014b; Mesquita et al., 2005). As such, if populations responded identically to landscape rearrangements and climatic conditions through time, one would expect the patterns of genetic diversity to be similar for co-occurring species, despite their intrinsic traits.

It is known that during the last glacial maximum (LGM, 0.018 Mya) the ice sheet reached the central part of the Portuguese territory, presumably as far as the Tagus River (Fig. 1A, Dias et al., 2000). Previous phylogenetic data, obtained with a calibrated molecular marker (cytochrome b), indicates that, at the time of the LGM, all contemporary species were already differentiated (Levy, Doadrio & Almada, 2008; Gante, 2011). Additionally, the extirpation of fish populations at northern latitudes and the persistence of Iberian fish species in southern refugia throughout the Quaternary had already been reported by several authors (e.g., Mesquita et al., 2005; Gante et al., 2009; Almada & Sousa-Santos, 2010; Araguas et al., 2013; Sousa-Santos et al., 2014b; Perea & Doadrio, 2015). Thus, as a consequence of the LGM, species inhabiting northern river basins should exhibit similar low levels of genetic diversity. Data presented in this paper shows that although this is true for L. bocagei, it is far from being a generalized pattern among northerly distributed species. The relatively high levels of genetic diversity observed in A. oligolepis, P. duriense and S. carolitertii could be explained by posterior recolonizations of northern streams by migrants from southern refugia, as suggested for many aquatic species (e.g., Hewitt, 2004; Gómez & Lunt, 2007; Provan & Bennett, 2008; Roselló & Morales, 2010; Oberdoff et al., 2011). However, this hypothesis is not plausible for the target species since at the time of the LGM most of the connections between river basins had already ceased (Pais et al., 2012). Since differences between species with identical distribution areas were detected for the three genetic diversity indices used, one must postulate than contemporary determinants of genetic diversity must have played a more significant role than historical ones.

Potential contemporary determinants of genetic diversity

As referred above, species intrinsic traits and landscape features may influence current levels of genetic diversity. Habitat loss and fragmentation, for instance, may result in declining population sizes and genetic diversity depletion, ultimately leading to local extinctions (Frankham, 1995; Spielman, Brook & Frankham, 2004; Ewers & Didham, 2006; Lowe & Allendorf, 2010). Husemann et al. (2012) found that frequent local extinction and re-colonization cycles in seasonal pools may even obscure historical signatures in the genetic patterns of some North American cyprinids as a result of repeated bottlenecks and population expansions. Demographic and genetic changes will impact species differently, according to their dispersal ability, life-history characteristics and habitat requirements (Faulks, Gilligan & Beheregaray, 2011).

Ecological traits such as tolerance to stagnant disconnected summer pools (Husemann et al., 2012), preference for flowing headwaters (Faulks, Gilligan & Beheregaray, 2011; Buonerba et al., 2015) or reproductive strategies (benthic vs. pelagic spawning, Osborne et al., 2014) were pointed out as determinants of genetic diversity in a vast array of freshwater fish. Our results showed that there are intrinsic characteristics (maximum size attained) and environmental characteristics (inter-basin connectivity and latitude) that are clearly determinant for genetic differences between populations. The variable “species” includes additional intrinsic characteristics other than those considered explicitly that could be causing the observed patterns of genetic variation. The relevance of this component is underlined by its significant correlation with genetic diversity indices and by its inclusion in the best fit model (together with “species maximum size,” “inter-basin connectivity” and “latitude” explained 26.7%–33.7% of the variance). The causes explaining the remaining 66.3%–73.3% of the genetic diversity variance remained currently unexplained and should ideally be addressed in future studies.