The Role of Directed Dispersal in Driving Genetic and Morphological Structure in Invasive Smallmouth Bass

Diedericks, Genevieve; Broeckhoven, Chris; von der Heyden, Sophie; Weyl, Olaf L. F.; Hui, Cang

doi:10.3389/fevo.2021.790829

ORIGINAL RESEARCH article

Front. Ecol. Evol., 17 January 2022

Sec. Behavioral and Evolutionary Ecology

Volume 9 - 2021 | https://doi.org/10.3389/fevo.2021.790829

The Role of Directed Dispersal in Driving Genetic and Morphological Structure in Invasive Smallmouth Bass

$\r\nGenevieve Diedericks,*$ Genevieve Diedericks^1,2*

Chris Broeckhoven³

Sophie von der Heyden²

Olaf L. F. Weyl^4† $Cang Hui,\r\n$ Cang Hui^5,6

¹Centre for Invasion Biology, Department of Botany and Zoology, Stellenbosch University, Stellenbosch, South Africa
²Evolutionary Genomics Group, Department of Botany and Zoology, Stellenbosch University, Stellenbosch, South Africa
³Department of Biology, Laboratory of Functional Morphology, University of Antwerp, Antwerp, Belgium
⁴Centre for Invasion Biology, South African Institute for Aquatic Biodiversity, Makhanda, South Africa
⁵Centre for Invasion Biology, Department of Mathematical Sciences, Stellenbosch University, Stellenbosch, South Africa
⁶Mathematical Biosciences Group, African Institute for Mathematical Sciences, Cape Town, South Africa

Dispersal is an essential life-history trait crucial to species persistence and diversification. This is particularly important in spatiotemporal fluctuating environments such as freshwater habitats, where species movement is confined to the dendritic network and wetted boundaries. To persist in such fluctuating environments, a species can modify, adaptively and plastically, its phenotypic variation to better match the environment or escape via directed dispersal to a more suitable habitat (i.e., matching habitat choice). We use the invasive smallmouth bass, Micropterus dolomieu, sampled at 10 km intervals, to assess the effect of directed dispersal on the fine scale genetic and phenotypic variation in populations of M. dolomieu along a river course. Gene flow was used as a proxy for dispersal. By comparing population genetic structure, morphological variation (of linear traits and geometric landmarks), and environmental heterogeneity, we discovered a clear correlation between environmental variation and morphological traits. Although isolation by distance seemed to have shaped the overall genetic pattern detected among the populations, the strong genetic structuring observed within the Ratel tributary appeared to be non-random. These results provide novel insights into the potential mechanisms promoting the spread and establishment of invasive species and the possible influence multiple introductions may have on fine scale genetic structuring.

Introduction

Dispersal is a life-history trait essential to the persistence, diversification and evolution of species (Edelaar and Bolnick, 2012; Saastamoinen et al., 2018). Although exchange of individuals between populations has typically been considered to occur at random (Patterson et al., 2008; Clobert et al., 2009), a paradigm shift is changing the way dispersal is perceived. No longer thought to occur at random, active dispersal can reflect the response of individual phenotypes to the surrounding environment (Bowler and Benton, 2005; Edelaar et al., 2008; Clobert et al., 2009; Edelaar and Bolnick, 2012; Hui et al., 2012; Jacob et al., 2015). As many environments are known to be both spatially heterogeneous and temporally variable (Bonte and Dahirel, 2017), individuals need to sample and compare different environments while making decisions on their next move (Edelaar et al., 2008; Zhang and Hui, 2014). Dispersing away from core populations not only benefits individuals with reduced competition intensity (Dytham, 2009; Rubenson and Olden, 2017) but may also provide opportunistic niches and resources for survival (Ronce, 2007; Clobert et al., 2009).

As a phenotype-environment mismatch is bound to occur in any heterogeneous environment, species could mitigate the mismatch by fine-tuning its phenotypic variation or simply preserving a general-purpose phenotype via phenotypic plasticity or local adaptation (Richards et al., 2006). Individuals can also move into more suitable environments via matching habitat choice (Edelaar et al., 2008; Edelaar and Bolnick, 2012; Camacho and Hendry, 2020). Matching habitat choice, as proposed by Edelaar et al. (2008) states that individuals will select and settle in environments that complement their phenotype, and in doing so drive directed dispersal. All three mechanisms can jointly affect local performance, population structure and genetic diversity (Arendt, 2015; Nicolaus and Edelaar, 2018). For instance, increasing the rate of random dispersal can genetically homogenize populations, by reducing genetic differentiation between populations (Haldane, 1948; Endler, 1973; Slatkin, 1985) and potentially limiting local adaptation (Edelaar and Bolnick, 2012; Jacob et al., 2015) and even speciation (Lucek et al., 2013). In contrast, directed dispersal (associated with matching habitat choice) is thought to promote genetic differentiation and potentially local adaptation when coupled with natural selection (Benkman, 2016), ultimately increasing local population genetic structure (Edelaar et al., 2008; Edelaar and Bolnick, 2012). Indeed, in a recent simulation study by Nicolaus and Edelaar (2018), incorporating both neutral and functional traits in addition to mild environmental fluctuations, the authors observed random dispersal to erode population genetic structuring, while non-random dispersal maintained strong population differentiation with little within-population variance. However, despite directed dispersal receiving increased attention, few empirical studies have sought to test matching habitat choice (Benkman, 2016; Camacho et al., 2020), especially for introduced species in their invaded ranges.

Lotic systems in hydrologically diverse regions are characterized by extensive spatiotemporal fluctuations in both biotic and abiotic conditions (Winemiller et al., 2010). Thus, fish that inhabit such environments are expected to match their environment and adapt swiftly to these pressures by maximizing their locomotive capabilities essential for feeding, spawning and predator avoidance (Via et al., 1995), which may be of particular importance in the face of climate change and anthropogenic activities. Nevertheless, most freshwater fishes often face movement constraints, as dispersal is limited to the wetted boundaries, and in riverine environments, the dendritic network (Garant et al., 2007; Olden et al., 2010). In this regard, numerous studies on freshwater fishes have revealed the correlation between functional traits and environmental variations (e.g., swimming and feeding performance; Keast and Webb, 1966; Vogel, 1994; Langerhans et al., 2003; Langerhans, 2008; Langerhans and Reznick, 2010). Indeed, Drakou et al. (2009) demonstrated that dispersal limitations play a bigger role in structuring fish communities than habitat heterogeneity.

Within freshwater ecosystems, invasive fish have been particularly favored as model organisms due to their superior dispersal abilities, essential for invasion success, when compared to native fishes (Lodge et al., 1998; Sala et al., 2000; Beisel, 2001; Rubenson and Olden, 2017). Moreover, it has been shown that dispersal opportunities, rather than species composition or favorable abiotic conditions of the recipient community, are an important component mediating aquatic invasions (Ricciardi and MacIsaac, 2000). This is likely because introduced fish are often from mixed genetic stock (Diedericks et al., 2018a,b), and unlikely to exhibit fine genetic structuring or matching habitat choice, as one would expect for native fish.

One such freshwater invader is the smallmouth bass, Micropterus dolomieu. Native to the east-central parts of the United States of America and two Canadian provinces, M. dolomieu is currently recognized as an invasive species in at least 12 countries worldwide, including South Africa (Loppnow et al., 2013). Although M. dolomieu has been described as a sedentary species when inhabiting riverine environments, establishing and guarding home ranges (Todd and Rabeni, 1989; Ridgway and Shuter, 1996), studies have observed remarkable fluctuations in dispersal behaviors when populations are subjected to environments with fluctuating habitat suitability, often leading to extensive movement in an attempt to colonize new habitats (Lyons and Kanehl, 2002; Hui and Richardson, 2017; Rubenson and Olden, 2017). Moreover, in a recent study by Diedericks et al. (2018a), significant population genetic structure was detected for M. dolomieu across the invaded South African range. However, little is known about the species’ genetic diversity and phenotypic variability at fine spatial scales, both essential components to the successful establishment and spread of an invasive species (Prentis et al., 2008).

Considering the aforementioned, we investigate the fine scale genetic- and phenotypic variation in populations of M. dolomieu within an invaded river system in South Africa. Following the matching habitat choice hypothesis, we predict: (1) M. dolomieu to display strong genetic population structuring with little within-population variance; and (2) a strong correlation between phenotype and environment along the river course. In contrast, should dispersal occur at random, we would expect homogenous genetic population structure across sampled populations and no phenotype-environment correlation. Understanding the interplay between morphology, genetics and the environment and how species movement may affect it, is of fundamental importance when considering future range expansions of invasive species.

Materials and Methods

Study System and Sampling

The Olifants River is one of the largest river systems within the Western Cape (WC) Province of South Africa, with an approximate length of 285 km and a catchment area of 46,220 km². The WC has a Mediterranean climate, with the Olifants River catchment receiving ∼ 90% of its rainfall in the winter months (Soderberg, 2003). In addition, the rain shadow effect sees the Cederberg Mountains to the East of the Olifants River receive double the amount of rain when compared to the main river valley (Soderberg, 2003; Figure 1). Extensive agricultural development and farming activities has occurred within the valley, leading to an increase in alien invasive vegetation and sedimentation, modified river bed, and reduced water quality (DWA, 2006; Supplementary Figure 1). Moreover, due to large amounts of water being used for irrigation in the middle reaches of the Olifants, large stretches of the main stem often dry up entirely during the summer months (DWA, 2006; Supplementary Figure 1). Despite the dire state of the Olifants River main stem, the tributaries and upper reaches still harbor near-pristine ecosystems (DWA, 2006). Two mainstream impoundments, namely Bulshoek and Clanwilliam Dam, constructed in 1919 and 1935, respectively, are situated within the Olifants River (Figure 1). Clanwilliam Dam possesses sluice gates and as such, prohibits upstream fish movement. Micropterus dolomieu, initially introduced for recreational fishing purposes, was released into the upper reaches of the Olifants River system in 1943 and 1945 (from the same breeding stock) and has since successfully invaded and established itself throughout the Olifants River system and its tributaries, including the two impoundments (Van der Walt et al., 2016). Although only one introduction event was historically recorded for M. dolomieu, Diedericks et al. (2018a) propose that multiple introductions offer a more likely scenario.

FIGURE 1

Figure 1. The sampling localities within the Olifants River system of the Western Cape, South Africa. All letters (A – J) represent sampled localities used for the linear and geometric morphometric measures, while localities sampled for the genetic analyses are indicated with an asterisk (*). Localities A – C represent those sampled in the Jan Dissels River tributary, while locality D represents the Clanwilliam Dam. Localities E – I represent sampled localities situated within the mainstem of the Olifants River and locality J represents the Ratel River tributary. The Clanwilliam Dam sluice gates are indicated with an arrow.

Approximately 20 M. dolomieu specimens were collected from ten localities spanning ∼100 km (total sample size, n = 203) in the austral summer (January – March) of 2015. Sampling occurred at approximately 10 km intervals and encompassed two tributaries, the Jan Dissels River (localities A-C) and the Ratel River (locality J) as well as the main stem of the Olifants River, including the Clanwilliam Dam (localities D-I; Figure 1). All specimens were collected by angling (CapeNature permit number: 0056-AAA043-00004) using standardized lure sizes, before being euthanized with clove oil (Ethical clearance reference number SU-ACUM14-00011, Stellenbosch University). A piece of muscle tissue was excised from the abdomen and stored in 70% ethanol for DNA extraction.

Environmental Data

Data characterizing the river topology, soil classification and hydroclimate for the Olifants River system were obtained from a near-global data set consisting of freshwater environmental variables (in a standardized 1 km grid) derived by Domisch et al. (2015). Given the potentially large spatial and temporal variation (Soderberg, 2003), and consequently large measurement error, localized environmental variables (e.g., pH, dissolved oxygen, salinity, etc.) were omitted. For each sampling locality, 16 variables (Supplementary Data Sheet 1) were extracted using the R packages RASTER 2.5-8 (Hijmans, 2016) and SP 1.2-4 (Pebesma, 2016). Firstly, four variables related to river topography were calculated. These include (1) average elevation across sub-catchment – “ELE,” (2) average river slope across sub-catchment – “RS,” (3) flow length (i.e., the number of upstream stream grid cells) – “FL” and (4) flow accumulation (i.e., the total number of upstream grid cells) – “FA.” Second, four variables that classify the soil type were obtained, including (5) sand content mass fraction across sub-catchment – “SNDPPT,” (6) silt content mass fraction across sub-catchment – “SLTPPT,” (7) clay content mass fraction across sub-catchment – “CLYPPT” and (8) coarse fragments (>2 mm fraction) volumetric across sub-catchment – “CFRVOL.” Lastly, we extracted eight bioclimatic variables, four characterizing air temperature and four describing precipitation averaged across sub-catchment, which were derived from the WorldClim database (Hijmans et al., 2005). Bioclimatic variables comprised (9) annual mean temperature – “BIO1,” (10) maximum temperature of warmest month – “BIO5, (11) temperature annual range – “BIO7,” (12) mean temperature of warmest quarter – “BIO10,” (13) annual precipitation – “BIO12,” (14) precipitation of wettest month – “BIO13,” (15) precipitation seasonality (i.e., coefficient of variation) – “BIO15” and (16) precipitation of wettest quarter – “BIO16” (Domisch et al., 2015). As Domisch et al. (2015) showed a strong correlation between water- and air temperature, the latter was used as a proxy for water temperature. Due to the highly correlated nature of the environmental variables, a principal component analysis (PCA) was conducted to avoid potential problems associated with multicollinearity. The PCA reduced the dataset to a smaller number of independent PC scores (hereafter referred to as “PC_ENVIRONMENT”), which were used in subsequent analyses.

Morphometric Measurements

Both linear and geometric measurements of morphometrics were employed to test for a correlation between phenotype and environment. All measurements were taken on fresh, unpreserved specimens. Firstly, a standardized photograph of the left lateral side of each fish was taken with a digital Sony DSC-HX100V camera. To ensure standardized images, the unpreserved specimens were pinned to laminated graph paper and a tripod, set to a standard height, was used. To quantify and visualize the body shape change between localities, 12 unambiguous and homologous landmarks were digitized using tpsDig2 version 2.17 software (Rohlf, 2013; Figure 2 and Supplementary Data Sheet 1). Twenty-five linear morphological characters, representing both swimming and feeding traits (Winemiller, 1991; Table 1, Figure 2, and Supplementary Data Sheet 1), were measured on site using digital calipers (model IP67, Mitutoyo Ltd., Japan) with a precision of 0.01 mm. In addition, an otolith from each fish was removed, cleaned and stored for aging analysis. Following the procedure described by Taylor and Weyl (2013), each otolith was sectioned and the growth rings counted to determine the age of each specimen. To ensure accurate age estimation, the specimens were aged by two of the authors (GD and OW) independently.

TABLE 1

Table 1. Mean and standard deviation for each of the linear morphometric measurements.

FIGURE 2

Figure 2. The 12 geometric morphometric landmark positions ( www.frontiersin.org ) and 21 of the 25 linear morphometric measurements taken for each specimen. The four linear measurements associated with the mouth (mouth height – MoH, mouth width – MoW, Length of snout with mouth open – LSO and Length of snout with mouth closed – LSC, are not illustrated. Linear morphometric measurement numbers correspond to those in Table 1.

A PCA conducted on all the linear morphometrics revealed one axis, coined “body size,” which explained 92.8% of the total variation (Eigenvalue = 23.2). To correct for size variation among populations, each linear morphometric variable was regressed against body size. The resulting residuals were kept for further analyses. Subsequently, a PCA was conducted on these residuals, retaining the most significant PC axes (hereafter referred to as “PC_LINEAR”). The PC_LINEAR axes were also compared between the localities using univariate F-tests, with age as covariate and the interaction effect of age*locality included. Bonferroni post hoc analyses were conducted to identify differences between localities.

To analyze the geometric landmark coordinates representing body shape variation, MorphoJ (Klingenberg, 2011) was employed. Firstly, a Procrustes superimposition was performed on the “raw” coordinates to remove orientation, size and position biases from the data (Klingenberg, 2011). The resulting Procrustes coordinates were regressed against centroid size (CS) to obtain size corrected residuals for each of the landmark coordinates. These residuals were then used to construct a covariance matrix, before a PCA was performed to reduce the dimensionality of the data (hereafter referred to as “PC_SHAPE”). Furthermore, univariate F-tests were employed to compare the PC_SHAPE axes between localities, with age and the interaction effect of age*locality included as covariates. Successive Bonferroni post hoc tests were conducted to determine differences in body shape between localities.

A generalized linear model was used to test for the effects of environment on morphology, with PC_LINEAR and PC_SHAPE scores as dependent variables and PC_ENVIRONMENT as independent variables. Statistical analyses were conducted in SPSS STATISTICS v. 20.0.0 (SPSS Inc., Chicago, IL, United States).

Genetic Analyses

For the genetic analyses, seven of the ten localities were included (n = 139) and encompassed the four populations in the two tributaries (localities A-C: Jan Dissels River; locality J: Ratel River), Clanwilliam Dam (locality D), and two main stem populations (localities H – I; Figure 1). Total genomic DNA was extracted from M. dolomieu muscle tissue using the NucleoSpin Tissue extraction kit (MACHEREY-NAGEL, Separations, Cape Town, South Africa) following the manufacturers protocol. All DNA extractions were kept at −20°C. To confirm that the collected individuals do not include hybrids of smallmouth x largemouth bass (see Diedericks et al., 2018b), two partial mitochondrial (mtDNA) gene regions [cytochrome oxidase subunit I (COI) and cytochrome b (cytb)] were sequenced. PCR reactions and cycling conditions followed those in Diedericks et al. (2018b).

To characterize genetic variation, population structure and gene flow among the sampled populations, individuals were genotyped using nine nuclear microsatellite markers (Mdo3, Mdo4, Mdo5, Mdo7, Mdo8, Mdo9, Mdo10, Mdo11, Lma21; Colbourne et al., 1996; Malloy et al., 2000). Three multiplex reactions were used to amplify the loci and followed the protocol described in Diedericks et al. (2018a,b); Supplementary Data Sheet 1. Microsatellites were genotyped on an automated sequencer (ABI 3730 XL DNA Analyzer, Applied Biosystems, CAF, Stellenbosch, South Africa) before being visually inspected, aligned and scored in Geneious^® 10.0.2 (Biomatters, Auckland, New Zealand).

To evaluate the dataset for amplification errors associated with stuttering and large allele drop-out, Microchecker 2.2.3 (Van Oosterhout et al., 2006) was used. Genepop 4.2.1 (Rousset, 2008), employing the Bonferroni method to correct for multiple comparisons (Rice, 1989), was used to assess all loci for deviations from Hardy-Weinberg equilibrium (HWE) and linkage disequilibrium. As most of the populations deviated from HWE, the EM algorithm (Dempster et al., 1977) implemented in FreeNA 1.2 (Chapuis and Estoup, 2007), was utilized to check for the presence of null alleles. Nuclear genomic variability within sampled populations was assessed as allelic richness (AR), number of alleles (Na), observed heterozygosity (H_O), expected heterozygosity (H_E), and Wright’s inbreeding coefficient (F_IS), calculated in FSTAT 2.9.3.2 (Goudet, 1995), Genepop 4.2.1 and Arlequin 3.5.2.2 (Excoffier and Lischer, 2010). The statistical significance of F_IS was determined after 1,000 permutations in FSTAT 2.9.3.2.

Genetic population structure was investigated using a five-pronged approach. Firstly, FSTAT 2.9.3.2 was used to calculate F_ST to assess genetic differentiation between sampled localities. Statistical significance was calculated using a Bonferroni correction. Secondly, an analysis of molecular variance (AMOVA), implemented in Arlequin 3.5.2.2 and statistical significance ascertained using 10,000 permutations, was conducted to reveal the extent of population subdivision among sampled populations. Thirdly, using allele frequencies we explored the genetic associations within and among the sampled localities with a PCA, conducted in the R package Adegenet 1.3.1 (Jombart and Ahmed, 2011). Fourthly, the Bayesian clustering approach, as implemented in STRUCTURE 2.3.4 (Pritchard et al., 2000), was used to identify and visualize the spatial population structuring along the river. Simulations employing the admixture model and assuming correlated allele frequencies, were run five times for each K (1 < K < 10). The initial burn-in was set to 70,000 Markov Chain Monte Carlo (MCMC) generations, followed by 300,000 MCMC iterations. The results were collated and the most probable K determined using the Evanno method (Evanno et al., 2005), as implemented in Structure Harvester 0.6.94 (Earl, 2012). Finally, Clumpp 1.1.2 (Jakobsson and Rosenberg, 2007) and Distruct 1.1 (Rosenberg, 2004) were used to generate and visualize the final results. Lastly, as the underlying STRUCTURE model assumes HWE and, according to the STRUCTURE manual, may not be suited to spatially distributed data with localized dispersal (Pritchard et al., 2010), we tested for spatial autocorrelation within our dataset, as suggested by Schwartz and McKelvey (2009). Spatial autocorrelation (r) investigates the spatial genetic structure among sampled individuals, using pairwise genetic distances over a given geographic distance (i.e., watercourse distance, km) (Smouse and Peakall, 1999). GenAlEx 6.5 (Peakall and Smouse, 2006) was used to assess the spatial autocorrelation among sampled individuals, analyzing the data at two different distance classes, i.e., (a) five classes of 20 km and (b) 10 classes of 10 km, with 9,999 permutations. This was done (a) to ensure large enough sample sizes for statistical power and (b) to prevent failure to detect fine scale genetic structure. Bootstrap re-sampling (10,000) was run to calculate the 95% confidence interval.

Connectivity and Dispersal

To test for isolation by distance (IBD) among populations, a Mantel test (Mantel, 1967) was conducted in GenAlEx 6.5. Significance was assessed using 10,000 permutations. Assignment tests have great power when wanting to detect migrants using multilocus genotypes (Manel et al., 2005). As such, two Bayesian approaches, both implemented in Geneclass 2.0 (Piry et al., 2004), were used to estimate connectivity and migration between sampled populations. Geneclass uses both genetic distance and allele frequencies, in conjunction with a MCMC resampling algorithm, to identify migrants (Piry et al., 2004; Putman and Carbone, 2014). Firstly, an assignment test (Cornuet et al., 1999), using the Bayesian criterion of Rannala and Mountain (1997), was conducted to determine whether all specimens originated from the localities in which they were sampled. This was done by comparing each individuals’ genotype to a simulated dataset. The simulation algorithm by Paetkau et al. (2004) was selected, as this method has been shown to accurately detect dispersing individuals even when all potential source populations were not sampled (Cornuet et al., 1999; Paetkau et al., 2004). The Type I error was set to 0.01 and 10,000 individuals were simulated. Secondly, to assess the dispersal of M. dolomieu within the Olifants River system, we used the “detect migrant” function in Geneclass 2.0 with the Bayesian criterion defined by Rannala and Mountain (1997). This method explicitly seeks and identifies first generation (“direct”) migrants (Paetkau et al., 2004; Piry et al., 2004), i.e., individuals originating from a population other than the one they were sampled from. As not all potential source populations may have been sampled, we selected the L_home/L_max likelihood computation method (Paetkau et al., 2004). The simulation algorithm (Paetkau et al., 2004), number of simulations (10,000) and the type I error (0.01) mimicked those used in the assignment test.

Lastly, to ascertain the effect of genetic variation and dispersal on the morphology-environment relationship, four canonical correspondence analyses (CCA) were run. All analyses were conducted at the individual level and comprised the residual linear morphometric values. Due to the correlated nature of the environmental variables, five variables were chosen to represent the river topology, soil classification and hydroclimate; annual mean temperature – BIO1, annual precipitation – BIO12, average river slope – RS, flow accumulation – FA, coarse fragments volumetric across sub-catchment – CFRVOL. The first CCA was used to explore the correlation between the morphological traits and the environmental variables (Supplementary Data Sheet 1). The second, third and fourth analyses, all partial CCA’s, investigated the same morphology-environment relationship, but considered the additional influence of genetic variation, dispersal and a combination of genetic variation and dispersal, respectively (Supplementary Data Sheet 1). Individual allele frequencies were used to represent the genetic variation, while the probabilities associated with the assignment test (i.e., probability of an individual originating from the sampled locality versus other localities) were used as a proxy for dispersal. Analyses were conducted using the “vegan” R package (Dixon, 2003), while the significance of these associations was calculated using the anova.cca function (Dixon, 2003) with 999 permutations.

Results

Fish Morphology

Significant variation in size (ANOVA; F₉,₁₉₁ = 38.37, P < 0.001) and age (ANOVA; F₉,₁₉₁ = 14.03, P < 0.001) was observed among localities. Specimens collected at locality D were approximately twice the size of specimens from other localities (Bonferroni post hoc tests, all P < 0.001). Size variation was rather limited after excluding specimens from locality D (15.1 – 21.8 cm SBL). Fish belonging to locality A and D were significantly older (Bonferroni post hoc tests, all P ≤ 0.03) than those of the other localities (Supplementary Figure 2), although no significant difference in age was observed among the other populations.

The PCA conducted on the Procrustes coordinates provided four PC_SHAPE axes (eigenvalues PC_SHAPE ₁ = 4.90; PC_SHAPE ₂ = 1.64, PC_SHAPE ₃ = 1.41, PC_SHAPE ₄ = 1.21), jointly explaining 67% of the variation, which were retained for the geometric morphometric data analyses. Shape PC1 (PC_SHAPE ₁) corresponds to a dorsal to ventral redistribution in body shape, while PC_SHAPE ₂ corresponds to rounded bodies in addition to shortened caudal peduncles (Figure 3). Similarly, PC_SHAPE ₃ represented fish with wider bodies, while PC_SHAPE ₄ denoted slender fish with large heads and shortened caudal peduncles (Supplementary Figure 3). The univariate F-tests revealed significant differences in PC_SHAPE ₁, ₂, ₄ among localities, although a statistically significant interaction effect age*locality was present in PC_SHAPE _2–3 (Supplementary Table 1). Bonferroni post hoc tests indicated that fish from locality A had higher bodies (PC_SHAPE ₁) than the other populations, with only marginal variation being observed with regards to locality B and D (Figure 3).

FIGURE 3

Figure 3. Graphs illustrating the mean and standard deviation for each locality with regards to (A) PC_{SHAPE 1} (dorsal-ventral redistribution), (B) PC_{SHAPE 2} (rounded bodies, shortened caudal peduncle), (C) PC_{LINEAR 1} (body depth), and (D) PC_{LINEAR 2} (gape size). Different alphabetical letters above the bars indicate values with statistically significant differences (Bonferroni post hoc test; P < 0.05), while identical letters indicate no significant differences among means. The fish outline drawings on the left of graph (A,B) depict the variation in fish body shape for PC_{SHAPE 1} and PC_{SHAPE 2}, respectively. The light blue line representing the average shape for all fish, while the dark blue line represents the upper and lower body shape extremity (scale factor set to 0.08 and – 0.08, respectively).

The PCA of the 25 linear morphometric measurements retained nine PC axes (PC_LINEAR _1–9; Supplementary Table 2). Principal component 1 (PC_LINEAR ₁) represented fish with higher bodies (MBH), higher and narrower caudal peduncles (CPdH and CPdW, respectively), an increased head height (HH), and decreased eye height (EH), anal fin height (AFH) and pelvic fin length (PFL). Hence, PC_LINEAR ₁ represents fish with “deeper bodies” (Figure 3). Principal component 2 (PC_LINEAR ₂) corresponded to “gape size,” as the highest loading scores were related to mouth height and width (MoH, MoW; Supplementary Table 2). The univariate F-tests revealed differences in all PC axes among localities except PC_LINEAR ₉ (Supplementary Table 1). Bonferroni post hoc tests further revealed that locality A had less deep bodies when compared to the other localities (Bonferroni post hoc tests, all P < 0.05, Figure 3). However, for PC_LINEAR _1–2 body depth differed between the age categories and a statistically significant interaction effect (Supplementary Table 1) was present between locality and age, indicating that the size-age relationship was locality specific.

The PCA of the 16 environmental variables retained two PC axes which explained 94.8% of the total variation. The first PC axis, PC_ENVIRONMENT ₁, was negatively correlated with elevation, slope, percentage of coarse fragments and temperature, but positively correlated with flow length, flow accumulation, precipitation and clay content (Supplementary Table 3). Therefore, PC_ENVIRONMENT ₁ represents a gradient from fast-flowing streams at high altitude to wide, slow-flowing rivers at lower altitude. PC_ENVIRONMENT ₂ was positively correlated with percentage of sand, as well as air temperature across the sub-catchment, but negatively correlated with elevation and silt content (Supplementary Table 3). Environment (PC_ENVIRONMENT ₁) had a significant effect on body depth (PC_LINEAR ₁; standardized β = 0.77; P = 0.01) and gape size (PC_LINEAR ₂; standardized β = −0.75; P = 0.01) (Figure 4) but not on the remaining size (PC_LINEAR _3–9; standardized β = −0.28 – 0.35; all P ≥ 0.32) or shape measurements (PC_SHAPE _1–4; standardized β = −0.26 – 0.46, all P ≥ 0.10). In contrast, PC_ENVIRONMENT ₂ had no effect on size (PC_LINEAR _1–9; standardized β = −0.28 – 0.54; P ≥ 0.11) and shape (PC_SHAPE _1–4; standardized β = −0.45 – 0.07; P ≥ 0.19) measurements.

FIGURE 4

Figure 4. Graphs illustrating the relationship between (A) body depth (PC_{LINEAR 1}) and (B) gape size (PC_{LINEAR 2}) and the environmental variables (PC_{ENVIRONMENT 1}) for each sampled locality (A – J).

Genetic Diversity and Population Structure

A total of 139 specimens, collected from seven localities, were successfully sequenced (for cytb and COI) and genotyped at nine microsatellite loci. When comparing these sequences to the Diedericks et al. (2018b) alignment and NCBI BLAST results, no hybrid or introgressed individuals were obtained; hence all 139 specimens were used in subsequent analyses. No amplification errors (i.e., stuttering, large allele drop-out), linkage disequilibrium or null alleles were detected in the microsatellite dataset. The majority of loci and populations did, however, not conform to Hardy-Weinberg expectations, likely due to the significant heterozygote deficit associated with the high levels of inbreeding (F_IS; Table 2). Observed heterozygosity (H_O) ranged from 0.150 (locality D) to 0.950 (locality H) across all loci, while the expected heterozygosity (H_E) ranged from 0.152 (locality A) to 0.756 (locality D). Marginally higher allelic richness values, based on a minimum of 15 specimens, was observed for the two mainstem localities (locality H and I) and the Ratel tributary (locality J; Table 2).

TABLE 2

Table 2. Genetic diversity measures (n = number of genotyped individuals; Na = number of alleles; AR = allelic richness for a minimum of 15 individuals; H_E = expected heterozygosity, H_O = observed heterozygosity; F_IS = inbreeding coefficient) for nine microsatellite loci amplified for each of the seven localities.

Pairwise F_ST values revealed varying levels of differentiation (F_ST = 0.0269 – 0.4320) between localities, with all but two comparisons statistically significant after Bonferroni correction (Supplementary Table 4). The AMOVA results assigned the largest amount of genetic variation to within each locality (78.98%), with limited variation being observed among tributaries (14.64%) or among localities within tributaries (6.38%). All three variance components were statistically significant (P < 0.001). The genetic PCA analysis, incorporating allelic frequencies, revealed two distinct clusters along the first axis (PC1 = 13.6%); the first cluster encompassing all specimens collected in the Ratel River tributary (locality J) and a small proportion of outlier individuals from locality C, H, and I, with the second cluster representing all other sampled localities (Supplementary Figure 4). Less variation was observed for the second axis (PC2 = 6.6%; Supplementary Figure 4). The Bayesian clustering method implemented in STRUCTURE identified two genetic clusters [K = 2, LnP(D) = −2461.26; Figure 5 and Supplementary Figure 5], one encompassing nearly all of the Ratel tributary specimens, and the other cluster representing all other localities, corroborating the genetic PCA results. Both spatial autocorrelation analyses revealed that M. dolomieu specimens displayed significantly more genetic structuring in the shortest distance classes (up to ∼18 km) than what is to be expected under the null model (Figure 6; all P < 0.007). In contrast, pairwise multilocus genotypes were significantly less genetically similar for all distance classes greater than 60 km (all P < 0.197; Figure 6). These results remained the same when removing locality J from the analyses (results not shown).

FIGURE 5

Figure 5. A graphical representation of the seven sampled localities used in the genetic analyses (A – D, H – J). The percentage correctly assigned individuals, obtained from GENECLASS 2.0, for each sampled locality is depicted with a pie chart, with colors corresponding to those in the legend on the right. The STRUCTURE plot represents all seven sampled localities, with each line in the plot corresponding to an individual within that sampled locality. Colors represent the proportion of an individual’s genotype assigned to a given genetic cluster (K = 2).

FIGURE 6

Figure 6. Multiallelic correlograms (blue line) showing the spatial genetic structure of M. dolomieu across two distance classes (10 and 20 km). The 95% null hypothesis confidence intervals are indicated by the red dotted lines. Positive values represent individuals more genetically similar than expected in a given distance class, while negative values imply that individuals in that distance class are less genetically similar than expected. The value in brackets denote the number of pairwise comparisons. The asterisk denotes the statistically significant pairwise relatedness in a specific distance class.

Connectivity and Dispersal Among Populations

A significant association (rxy = 0.462, P = 0.001, R² = 0.214) between geographical distance and genetic divergence was observed across all sampled localities (Supplementary Figure 6). Nine first generation migrants were detected with GENECLASS 2.0; two each for localities A (A3, A9), I (I4, I5) and J (J9, J18) and one each for localities C (C14), D (D19) and H (H2). No migrant individuals were detected in locality B. The majority of these first-generation migrants were males (6/9), with 8/9 migrants belonging to the 2–3-year age class (Supplementary Data Sheet 1). All migrants originated from either locality H, I or J. The assignment test revealed that 54% of all individuals were correctly assigned when considering the locality of highest probability (Figure 5). The majority of misassigned individuals were collected from localities B, C, D, and H (Figure 5). All nine first generation migrants formed part of the misassigned individuals.

The first CCA revealed that the overall morphology corresponded significantly to environmental variation (F = 3.47, P = 0.001; Figure 7). River slope (RS: F = 10.39, P = 0.001), and to a lesser extent flow accumulation (FA: F = 2.45, P = 0.02) and annual mean temperature (BIO1: F = 2.50, P = 0.02) were statistically significant correlates of variation in morphology. The partial CCA, accounting for the additional effects of genetic variation, was once again statistically significant (F = 2.51, P = 0.001), but revealed less influence of climatic variation on morphology (RS: F = 6.67, P = 0.001; FA: F = 2.31, P = 0.03; BIO1: F = 1.95, P = 0.06; Figure 7), confirming the lack of genetic influence on the morphology. When accounting for the effect of dispersal and genetics*dispersal, however, the morphology-environment correspondence was no longer statistically significant (dispersal: F = 1.08, P = 0.30; genetics*dispersal: F = 0.93, P = 0.53; Figure 7). In these two partial CCAs the morphology-environmental correspondence became weak (dispersal – RS: F = 2.38, P = 0.02; FA: F = 0.82, P = 0.05; BIO1: F = 0.99, P = 0.35; genetics*dispersal – RS: F = 1.98, P = 0.04; FA: F = 0.81, P = 0.57; BIO1: F = 0.70, P = 0.68).

FIGURE 7

Figure 7. The canonical correspondence analysis biplots showing the relationship between morphological traits and environmental variables (A) without considering genetics, (B) taking genetic variation into account, (C) taking dispersal (via gene flow) into account, (D) considering both genetic variation and dispersal. Five environmental variables representing river topology, soil classification and hydroclimate were used, namely: annual mean temperature – BIO1, annual precipitation – BIO12, average river slope – RS, flow accumulation – FA, and coarse fragments volumetric across sub-catchment – CFRVOL. Morphological trait abbreviations correspond to those used in Table 1.

Discussion

Lotic systems in hydrologically diverse environments, such as in Mediterranean climate regions, are known to be spatially and temporally variable (Winemiller et al., 2010; Cooper et al., 2013). Thus, in order to persist in such changing environments, species either need to adapt or move to more suitable habitats when changes occur. However, despite directed dispersal gaining attention in the literature, few studies to date, particularly those concerning invasive species, have attempted to unravel the effect directed dispersal may have on fine scale genetic- and phenotypic variation among populations. In the present study, we used the invasive smallmouth bass, M. dolomieu, as study organism to investigate the influence of individual dispersal on the fine scale genetic- and phenotypic variation along a river course. Our results revealed a strong correlation between phenotype and environment, independent of genetic variation. Although significant genetic population differentiation was detected among populations, limited dispersal (as inferred from microsatellite data) suggests that matching habitat choice is unlikely to have contributed to the successfully establish of the smallmouth bass in the Olifants River.

Correlates in Morphological and Environmental Variation

Combining morphometric measures with environmental variables obtained from a near-global database, our results support the prediction that comparable environments will have similar phenotypes, particularly with reference to the linear morphometric measures, with both body depth (PC_LINEAR ₁) and gape size (PC_LINEAR ₂) being correlated to environmental variables (PC_ENVIRONMENT ₁; Figure 4). Despite the significant interaction effect between locality and age, examination of the within population variation revealed that all individuals belonging to locality A had high scores for PC_LINEAR1 regardless of age (Supplementary Figure 2). No significant correlation was, however, detected between any of the environmental variables and body shape (PC_SHAPE _1–4), likely due to the large amount of within-locality variation. The arching of the body, as described by PC_SHAPE ₁, may have been due to slight differences in specimen posture during landmark capture, as shown by Valentin et al. (2008).

The Jan Dissels tributary (localities A – C) is characterized by high-flow environments, evident from the steeper slope and the higher percentage of pebbles and rocks (CFRVOL; Johnson, 2004). It is thus not surprising that M. dolomieu from locality A (and subsequently localities B and C) were significantly more streamlined when compared to the other populations (Figure 3), supporting the hypothesis of Langerhans (2008) stating that streamlined fish will be present in high-flow environments, as this body shape promotes steady swimming in fast flowing waters. Despite the distinct morphology of the fish collected in the Ratel tributary (locality J – PC_SHAPE ₂, and to a lesser extent PC_LINEAR ₁; Figure 3), no association could be made with habitat (Figure 4). This is likely because of the broad scale environmental variables (1 km² grid) used in the present study and the close proximity of locality J to the main stem. Similarly, substantial variation between localities, though not statistically significant, was observed in gape size (PC_LINEAR ₂; Figure 3), with localities A – C, characterized by a relatively smaller gape size, clustering together. Relative gape size was found to decrease with a decrease in clay content mass (CLYPPT; Figure 4), and consequently an increase in coarse fragments. Invertebrates, especially crabs, are known to prefer coarse substrate habitats situated in pristine, high flow environments (Zimmerman and Covich, 2003; DWA, 2006), as present in the Jan Dissels tributary. A smaller gape and consequently stronger suction (Carroll et al., 2004; Wainwright et al., 2007; Day et al., 2015), might be advantageous in this environment as a stronger suction is required to feed on benthic prey items, such as crabs (Day et al., 2015). In contrast, a relatively large gape is better suited for ram feeding and predating on fish in the water column (Carroll et al., 2004; Wainwright et al., 2007; Day et al., 2015). This may well be the case in the mainstem of the Olifants River system, as numerous co-occurring fish species are present (Van der Walt et al., 2016). Moreover, the fact that M. dolomieu from the Jan Dissels tributary have more fusiform bodies compared to the other localities further supports this notion, as steady swimming is needed to get close enough to benthic prey (Sass and Motta, 2002; Higham, 2011). In contrast, thrust and maneuverability, associated with relatively deep bodies (as seen within localities D – J; Figure 4), provides M. dolomieu with a fast approach and attack, essential to capture prey in the water column (Sass and Motta, 2002; Tran et al., 2010).

Population Structuring and Connectivity Among Micropterus dolomieu Populations

In the present study, we predicted that, following directed dispersal, one would observe marked population structuring and little within population variance. Our results corroborate this to some extent, as significant population differentiation was detected. However, this was most likely a result of IBD (Supplementary Figure 6), with more genetic structuring than expected over short distances, and less genetic structuring than expected over large distances, as revealed by the spatial autocorrelation (Figure 6). Furthermore, the low number of first-generation migrants detected further support this notion. These findings are in line with previous work, describing riverine M. dolomieu inhabitants as sedentary, establishing and guarding home ranges and displaying nest site fidelity (Todd and Rabeni, 1989; Ridgway and Shuter, 1996). In contrast to our prediction, high levels of within locality variation were observed. This can potentially be explained by two alternative hypotheses, as suggested by Barluenga and Meyer (2005). Firstly, a large heterogenous lineage of smallmouth bass could have colonized the Olifants River, differentiating into several groups along the diverse river course. Secondly, multiple genetically diverse lineages of smallmouth bass colonized the Olifants River, homogenizing the gene pool post introduction. Considering M. dolomieu’s invasion history and the supposedly multiple introductions of smallmouth bass into South Africa, and likely the Olifants River system (Diedericks et al., 2018a), the latter hypothesis seems most probable. Moreover, the existence of a Wahlund effect (Wahlund, 1928), i.e., the presence of (cryptic) cohorts in what is believed to be one genetically homogeneous sample, further supports the hypothesis of multiple introductions. The underlying causes for a Wahlund effect are thought to be geographic barriers and genetic drift, although relatedness and inbreeding could also contribute (Castric et al., 2002). Several aspects of the current study system are in line with the characteristics of a Wahlund effect. Firstly, the presence of a dispersal barrier, the Clanwilliam Dam sluice gates, prohibits movement for fish from localities A – C to the main stem. Secondly, genetic drift could have been promoted through sudden immigration, as would be expected during an introduction event (Lee, 2002). Thirdly, as fish are known to shoal with kin (Brown and Brown, 1996; Gerlach et al., 2001; Castric et al., 2002; Dehais et al., 2010) it is highly probable that genetically related individuals were sampled in the present study, as suggested by the similarity in age classes within each locality (Supplementary Data Sheet 1). Lastly, the high levels of inbreeding detected could be ascribed to non-random mating, driven by nest site fidelity (Dehais et al., 2010), a trait displayed by smallmouth bass (Todd and Rabeni, 1989; Ridgway and Shuter, 1996). Despite the aforementioned, the Ratel tributary (locality J) was shown to be genetically distinct when compared to the rest of the individuals, as suggested by the STRUCTURE and genetic PCA results (Figure 5 and Supplementary Figure 4). Moreover, the assignment test showed limited dispersal from the Ratel to downstream populations despite the absence of any geographical barriers. Numerous studies have revealed similar patterns. For example, higher genetic differentiation was found in introduced smallmouth bass inhabiting the tributaries of the Susquehanna River basin compared to the main stem (Schall et al., 2017). In addition, Hendry et al. (2002) found high genetic divergence and little gene flow among stickleback populations inhabiting an inlet stream and a lake. Likewise, Carlsson et al. (1999) demonstrate fewer brown trout exchanges between a tributary and main stem than within the main stem itself. Tributaries are often ecologically differentiated thereby supporting locally adapted populations that rarely disperse into the main stream (Dehais et al., 2010). In contrast to the Ratel and Jan Dissels tributaries which sustain permanent water flow, individuals inhabiting the mainstem of the Olifants River face the risk of becoming trapped as sections of the river become disconnected and dry up during summer. This may force M. dolomieu to disperse either up- (locality J) or downstream (locality D and subsequently A-C if the sluice gates are opened) to other freshwater habitats with better environmental conditions (Whitledge et al., 2002; Westhoff et al., 2016), and could ultimately result in the limited population structuring within the main stem.

Conclusion

The results of our study reveal distinct genetic and morphological differences between populations of M. dolomieu. On the one hand, a strong correlation was found between morphological traits and environmental conditions, irrespective of the observed genetic variation among the populations. On the other hand, the populations were genetically differentiated due to isolation-by-distance (IBD). This IBD pattern suggests that gene flow, and thus dispersal, is more likely to occur between closely situated populations than those at greater spatial distances. Limited dispersal between habitats could result from matching habitat choice (Edelaar et al., 2008), provided individuals have the opportunity to select between different habitat types. In our study system, such instances appeared to be rare as habitats of neighboring localities were often very similar, or separated by a barrier (e.g., Jan Dissels). A notable exception is the Ratel tributary, which represents a morphological and genetically distinct population within the Olifants River. The strong genetic structuring observed within the tributary appears to be non-random (i.e., directed gene flow from the neighboring populations toward the tributary as revealed by the migrant analysis; Figure 5), reflecting the proposed genetic pattern suggested by the matching habitat choice hypothesis (Edelaar et al., 2008; Nicolaus and Edelaar, 2018).

In conclusion, our study highlights the importance of gene flow and dispersal in shaping the fine scale population genetic variation, particularly with reference to invasive species. The use of whole genome markers such as SNPs (Tiffin and Ross-Ibarra, 2014), coupled with more comprehensive spatial sampling schemes and detailed spatio-temporal environmental data, will be pivotal to understand the relative importance of directed dispersal, local adaptation, phenotypic plasticity and random genetic differentiation for trait variation.

Data Availability Statement

The datasets presented in this study can be found in the Supplementary Material.

Ethics Statement

The animal study was reviewed and approved by the Research Ethics Committee: Animal Care and Use (REC: ACU), Stellenbosch University.

Author Contributions

GD, SH, OW, and CH contributed to conception and design of the study. GD and CB conducted the sampling and took all the linear measurements. GD conducted all the genetic lab work and analyses, landmarked all specimens, and wrote the first draft of the manuscript. GD and OW aged the fish. CB performed all the statistical analysis. All authors contributed to manuscript revision, read, and approved the submitted version.

Funding

This work was supported by the Department of Science and Technology (DST) and National Research Foundation (NRF) of South Africa (CH: 89967 and 109244; OW: 110507 and 109015) and the DST-NRF Centre of Excellence for Invasion Biology (GD).

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.

Acknowledgments

In loving memory of OW who passed away on 14 November 2020 while chasing fish in mountain streams. We thank Liz Alter and Erika Rubenson for their useful comments and suggestions.

Supplementary Material

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

References

Arendt, J. D. (2015). Effects of dispersal plasticity on population divergence and speciation. Heredity 115, 306–311. doi: 10.1038/hdy.2015.21

PubMed Abstract | CrossRef Full Text | Google Scholar

Barluenga, M., and Meyer, A. (2005). Old fish in a young lake: stone loach (Pisces: Barbatula barbatula) populations in Lake Constance are genetically isolated by distance. Mol. Ecol. 14, 1229–1239. doi: 10.1111/j.1365-294X.2005.02468.x

PubMed Abstract | CrossRef Full Text | Google Scholar

Beisel, J. N. (2001). The elusive model of a biological invasion process: time to take differences among aquatic and terrestrial ecosystems into account? Ethol. Ecol. Evol. 13, 193–195. doi: 10.1080/08927014.2001.9522785