Introduction

Freshwater environments are naturally and anthropogenically fragmented because of physical features such as mountain ranges, waterfalls, cascades, dams and agriculture. The degree to which fragmentation affects species is directly influenced by their life-history traits and habitat preferences (Frankham et al., 2010; Osborne et al., 2014). For example, fragmentation can have a negligible effect in species with high dispersal abilities that are able to cross unsuitable landscapes between habitat patches (Boizard et al., 2009; Faulks et al., 2010). However, in species with low dispersal abilities, migration between patches can become impossible (Bilton et al., 2001) that can result in severe demographic reductions and local extinctions, especially when habitat fragmentation is accompanied by habitat loss and degradation (Fahrig, 2003; Fagan et al., 2005; McCraney et al., 2010). Such scenarios are expected to lead to loss of genetic variation and inbreeding depression in populations of poor dispersers, negatively affecting their adaptive potential and long-term persistence (Hanski and Gaggiotti, 2004).

Landscape genetics contributes to the management of fragmented metapopulations by providing an environmental angle to identify associations between landscape features and ecological and evolutionary processes. Landscape genetics combines information from spatial statistics and population genetic data sets to identify factors affecting population connectivity (Manel et al., 2003). As such, it can clarify relationships between a species and its contemporary environment. Riverine environments provide ideal systems for landscape genetic surveys, as dispersal of obligate aquatic organisms is necessarily restricted to the river network. Thus, dispersal corridors and the associated habitat features can be better defined compared with many terrestrial and marine study systems (Wang et al., 2009).

Landscape genetics allows assessment of connectivity, both functional (dispersal and gene flow) and structural (distance between patches and differing habitats), enabling statistical assessment of models of population structure. One such model is the stream hierarchy model that predicts that dispersal of riverine organisms is influenced by riverine distances between populations, and is therefore hierarchically constrained by the dendritic nature of the river basins (Meffe and Vrijenhoek, 1988; Fausch et al., 2002; Hughes et al., 2009). Being able to test predictions from spatial models of population structure can greatly improve our knowledge of the movements and habitat preferences of a species. This allows us to better direct conservation initiatives, predict metapopulation scenarios and plan for future threats (see, for example, Hopken et al., 2013).

To date, most landscape genetic studies on freshwater fishes (that is, riverscape genetics) have focused on medium- to large-bodied species with relatively high migration potential (Wofford et al., 2005; Whiteley et al., 2006; Faulks et al., 2010). Very few studies have been conducted on benthic, low dispersal species that display traits that make them more susceptible to habitat fragmentation (McCraney et al., 2010; Hopken et al., 2013; Roberts et al., 2013). In addition, most studies have predominantly focused on the impact that barriers or landscapes variables have on population connectivity, or gene flow, with only a few assessing the impact of the local environment on population genetic diversity (Cena et al., 2006; Faulks et al., 2010, 2011). Yet, it is these landscape genetic studies that identify which species-specific habitats are best able to support larger, and consequently, healthier populations that are likely to provide the most valuable information to guide basin-wide conservation efforts.

Here we use population and landscape genetic approaches to study a surrogate indicator species. Species-based biodiversity surrogates, such as indicator species, provide valuable alternatives for assessing the complex impacts of habitat fragmentation on ecosystems (Rodrigues and Brooks, 2007). Indicator species can be treated as proxies that monitor or reflect the presence, abundance or richness of other elements of biota, critical ecosystem processes or particular environmental impacts (Lindenmayer and Likens, 2011). We assess the influence of riverine configuration and habitat features on the population structure and genetic diversity in a poor dispersal and threatened freshwater fish from Australia’s vast Murray–Darling Basin (MDB). The MDB is Australia’s most important agricultural region, accounting for 41% of the national agricultural income and 60% of its water use (CSIRO, 2008), despite being characterized by low and variable rainfall and severe droughts. Reflecting these paradoxical attributes, the MDB is one of the world’s most degraded freshwater systems (Palmer et al., 2008). Large modification to flow regimes, extensive agriculture, land clearing, invasive species, in-stream barriers to fish passage and habitat modification throughout the MDB have together led to a highly fragmented environment (Jackson and Williams, 1980; Allison et al., 1990; Lintermans, 2007). To alleviate the anthropogenic pressures placed on threatened species in the MDB, an understanding of the key environmental features influencing population demography is clearly needed before sound conservation initiatives can be designed and implemented.

Our study species is the river blackfish Gadopsis marmoratus (Percichthyidae), and specifically the small-growing form found in the MDB (size <350 mm total length; candidate species ‘NMD’, sensu Hammer et al., 2014). All members of the G. marmoratus species complex are regarded as a habitat specialist, with a strong preference for microhabitats with low levels of sediment, slow-flowing, deep pools and high structural integrity from undercut banks, woody debris and large rocks (Khan et al., 2004; Koster and Crook, 2008). Furthermore, they have a relatively small home range, low fecundity and demersal larvae that inhabit dense cover (Jackson, 1978; Khan et al., 2004; Koster and Crook, 2008). Populations of G. marmoratus across the MDB have declined dramatically in their extent of occurrence and abundance since European settlement, especially for heavily regulated large river environments and areas with less reliable rainfall (Lake, 1971; Lintermans, 2007). This decline is formally recognized at the northern and southern edges of the MDB with the species protected under fisheries legislation in Queensland and South Australia, respectively, with other regional populations outside the MBD also listed as threatened (that is, upper Wannon River, Victoria and Snowy River, New South Wales). Several genetic data sets are available for the G. marmoratus complex (Sanger, 1984; Ovenden et al., 1988; Miller et al., 2004; Ryan et al., 2004; Hammer et al., 2014), but these were not directed at analyses of population genetic variation, connectivity and landscape genetics.

Limited dispersal capabilities, threatened status and anthropogenically altered habitat make G. marmoratus an ideal surrogate indicator species on which to use a landscape genetics approach to explore how habitat features and fragmentation influence ecosystem-wide contemporary genetic architecture. To provide the appropriate detailed genetic resolution, we employed a suite of recently developed microsatellite markers (Arias et al., 2013) and a comprehensive sampling of catchments and localities across the MDB. We expect G. marmoratus populations to exhibit low genetic diversity and high levels of population structure. We further predict that genetic diversity will be positively influenced by the frequency of slow-flowing pools and perenniality as measures of microhabitat availability (see, for example, Bond and Lake, 2005), and negatively influenced by sedimentation and riparian damage (habitat limiting). The anthropogenic pressures currently placed on G. marmoratus in the MDB are likely to persist into the future, a scenario exacerbated because of global climate warming (Balcombe et al., 2011). Therefore, we use our findings to develop targeted management recommendations that prioritize conservation resources across the basin while maintaining and enhancing population genetic diversity. Our results should shed further light on the influence that habitat quality and quantity imposes on habitat specialists with similarly low dispersal, plus provide practical guidelines for prioritizing conservation management actions, at a range of spatial scales, for organisms inhabiting highly fragmented aquatic ecosystems.

Materials and methods

Study system and samples

The MDB is an expansive river system, occupying most of south-eastern Australia (>1 million km2). The broad climate shifts from cool temperate in the south to semi-arid in the north and west, with a fringe of cooler upland habitat in two areas of topographic relief: (1) the Great Dividing Range in the east and south, the third longest land-based ridge in the world; and (2) the smaller Mount Lofty Ranges (elevation 200–480 m). The majority of stream flow is generated in these upland areas (Leblanc et al., 2012), with a series of major tributaries arising along the Great Dividing Range draining into two main rivers, the Murray and the Darling. The lower Murray River exits into the Southern Ocean via the Lower Lakes. Small tributary streams flow off the eastern Mount Lofty Ranges into the Lower Lakes and lower Murray River (Figure 1).

Figure 1
figure 1

Sampling sites for the river blackfish G. marmoratus in the MDB. The top inset shows the location of the MDB study region in Australia, the lower insets show further detail on sites in (a) the Mount Lofty Ranges and (b) central Victoria. Sites are grouped into one of three major regions in the MDB (see key) and three letter site codes follow Supplementary Table S1.

Full size image

A total of 332 G. marmoratus individuals (whole fish in liquid nitrogen or fin clips in 100% ethanol) were collected between 1999 and 2013 from 29 sites encompassing 14 catchments and 3 broad regions of the MDB, namely the Upper Darling, Upper Murray and lower Murray (Figure 1 and Supplementary Table S1). Five sites in the lower Murray, the subject of regular threatened species monitoring (Hammer et al., 2013), were represented by multiple temporal-sampling events to reach a suitable sample size (MAR, ROD, ARS, NAN, TOO: Supplementary Table S1). In three small stream systems, three samples represented a composite of spatially proximate sites (<1 km apart: TOO, ARS, KBC). All genotypes from the 29 sites were retained for individual-level analyses (Supplementary Table S1), whereas population-level analyses were carried out only on 12 site samples that had ⩾8 individuals (Table 1).

Table 1 Summary statistics for Gadopsis marmoratus in the MDB based on 12 microsatellite DNA markers
Full size table

Genetic variation

Genomic DNA was extracted using a modified salting out protocol (Sunnucks and Hales, 1996). Twelve polymorphic microsatellite loci developed specifically for the MDB G. marmoratus were amplified via the PCR: Gama01, 02, 03, 04, 05, 06, 07, 08, 09, 10, 12 and 13 (Arias et al., 2013). Reagent concentrations and the touchdown PCRs follow Beheregaray and Sunnucks (2000). The PCR products were screened on an ABI 3730xl automated DNA sequencer (Applied Biosystems, Foster City, CA, USA) and alleles called using GENEMAPPER 4.0 (Applied Biosystems). MICRO-CHECKER 2.2.3 (Van Oosterhout et al., 2004) was used to assess scoring errors and null alleles.

Fisher’s exact test for linkage disequilibrium between all pairs of loci and an exact test for deviations from Hardy–Weinberg equilibrium were carried out in GENEPOP 4 (Raymond and Rousset, 1995), adjusted for multiple comparisons using sequential Bonferroni (Rice, 1989). Rarefied allelic richness (Ar) and rarefied private allelic richness (Pa) (controlled for a population size of five) was calculated using HP-RARE (Kalinowski, 2005). The latter uses rarefaction to compensate for sampling bias due to a larger number of alleles expected in large sample sizes compared with small sample sizes (Leberg, 2002). The observed (HO) and expected heterozygosity (HE) per population and per locus was calculated in ARLEQUIN 3.5.1.3 (Excoffier and Lischer, 2010). The inbreeding coefficient (Fis) was calculated using FSTAT 2.9.3 (Goudet, 2001) for each population. A nonparametric Wilcoxon signed-rank test was used to assess differences in genetic diversity (Ar and HE) between populations.

Genetic structure and assignment tests

The allele size permutation test in SPAGeDi 1.3d (Hardy and Vekemans, 2002) indicated that Fst was more appropriate than Rst as an index for assessing population differentiation. Thus, population pairwise comparisons were carried out in Arlequin using Fst that was also used to assess temporal variation in population structure at sites with multiple samples (Supplementary Table S1). The latter disclosed no significant differentiation and thus samples from different years were pooled for further analyses. In addition, an analysis of molecular variance (AMOVA) based on Fst was conducted in ARLEQUIN to assess the fit of predictions from the stream hierarchy model by partitioning genetic variation by differences between sites within streams, as well as streams within catchments and catchments within regions. The significance levels for all Fst and AMOVA comparisons were assessed using 10 000 permutations.

STRUCTURE 2.3.4 (Pritchard et al., 2000) was used to determine the most likely number of population clusters using the correlated allele frequency model, and not using sampling locations as priors. Twenty STRUCTURE runs were carried out, each with 20 separate iterations and 1 million Markov chain Monte Carlo replicates after an initial burn-in period of 300 000. The most likely number of clusters (K) was chosen following Evanno et al. (2005) and assessed with STRUCTURE HARVESTER (Earl and vonHoldt, 2012). The results were combined using CLUMP1.1.2 (Jakobsson and Rosenberg, 2007) and visualized using DISTRUCT 1.1 (Rosenberg, 2004). All above analyses were conducted using 12 population samples.

Individual multilocus assignment tests were carried out using the Bayesian approach of Rannala and Mountain (1997). Assignment tests were used to determine the population of origin of fish from sites with small sample sizes and to assess the proportion of individuals assigned to their sampling population. The latter assignment is based on higher than average probability of being born locally and provides a measure of natal site philopatry (see, for example, Möller and Beheregaray, 2004). This was carried out in GENECLASS 2 (Piry et al., 2004) by simulating 10 000 individuals and using an assignment threshold of 0.01

Spatial population structure

Isolation by distance (IBD) predicts genetic differentiation as a function of geographic distance and is considered the baseline for evaluating more complex metapopulation scenarios that may exceed the ability of IBD to represent spatial pattern in genetic structure (Jenkins et al., 2010). We assessed IBD at regional scales with Mantels tests in IBDWS (Bohonak, 2002) by comparing population pairwise matrices of Fst with stream distance matrices calculated in ARCMAP 10.1 (ESRI 2012, Redlands, CA, USA). Three separate regional IBD tests were carried out: across all sites, across the Upper Murray and across the lower Murray. The IBD test was not carried out for the Darling region because it was represented by only three sites.

Analysis of migration

Contemporary levels of connectivity were assessed by estimating migration rates between inferred population clusters using BAYESASS 3.0 (Wilson and Rannala, 2003). BAYESASS estimates recent (last two generations) migration rates between population pairs. The assumptions associated with this analysis (that is, no linkage between loci, Fst values >0.05) hold for our population samples. Ten million iterations were carried out for the Markov chain Monte Carlo analysis, with the initial 1 million being discarded as burn-in. Samples intervals of 100 were used to estimate posterior probabilities. Delta values for migration rates, allele frequencies and inbreeding coefficients were set at 0.15, 0.46 and 0.65, respectively. Different random seeds were used to carry out six independent runs.

Influence of genetic differentiation on genetic diversity

This was assessed following Coleman et al. (2013), who suggested that prioritizing management on differentiation alone may be detrimental to the maintenance of genetic diversity because of the strong negative relationship usually found between population differentiation and diversity. Population-specific Fst was calculated using a hierarchical Bayesian approach in Geste 2.0 (Foll and Gaggiotti, 2006) to determine the level of differentiation of each population relative to other populations. Linear regression was then used to assess the relationship between genetic uniqueness (population-specific Fst) and diversity (Ar and HE). The associated normality assumptions were tested using the Shapiro–Wilk test and independent histograms in IBM SPSS Statistics 2.1 (Chicago, IL, USA). Linear regression was carried out on Ar and HE and associated significance tested in R 3.0.1 (Vienna, Austria). As Ar and HE were highly correlated (results not shown), only Ar was used for subsequent analysis, as it is more sensitive to recent changes in population size and genetic diversity (Nei et al., 1975). The significance and strength of the relationship between Fst and Ar was also determined as above.

Landscape genetics

Landscape genetics was used to determine the impact of five habitat features on genetic variation (Ar) of demes from the 12 sites with larger sample sizes. Surrogates were used to represent meaningful habitat features based on expert opinion and field observations. These were the sum of the average proportion of grid cells used for forestry, intensive animal and plant production (that is, land use) which were used as a surrogate for sedimentation (Trimble and Mendel, 1995; Croke and Hairsine, 2006). The average difference in altitude over each upstream grid cell (that is, slope) was used as a surrogate for the frequency of deep slow-flowing pools. The distance between the nearest upstream and downstream dam (that is, barrier) was used as a surrogate for potential river area. The percentage flow contribution to mean annual discharge by the six driest months of the year (that is, perenniality) was used to represent fluctuations in available habitat. Finally, the percentage of the grid cells where industrial, agricultural, urban and point source pollution exists (that is, modified land) was used as an indicator of overall human impacts. Data were sourced from the National Catchment and Stream Environment Database 1.1.5 (Stein, 2012), available as part of the Geofabric products associated with the 9 s Digital Elevation Model (DEM)-derived streams and the National Catchment Boundaries. The segment number from the GEODATA national 9 s DEM that corresponds to sample sites was identified in ARCMAP. These were used to identify environmental characteristics linked to the stream layer in the database using ARCMAP. The mean value from all upstream grid cells was used for the predictor variables as it is likely that habitat features upstream from the sampling site have an impact on habitat in the downstream segments. Normality of all response variables was assessed as outlined above and transformations were made where necessary. Collinearity assumptions were assessed by examining the significance of pairwise correlations between all predictor variables. Perenniality was significantly correlated with both slope and land use (Supplementary Table S2) and modified land could not be transformed to meet normality. These two variables were removed from further analyses.

The log likelihoods of the linear modes for all possible combinations of the three candidate environmental variables influencing genetic variation was calculated in the all.regs function of R hier.part package (Walsh and Macnally, 2008). The best models were then selected using AICc (Akaike’s information criterion adjusted for small sample sizes; Akaike, 1973) and the hierarchical partitioning function in the hier.part package used to assess the individual contribution of each variable to the full model. Recent criticism has been made about the use of linear regression and AICc in landscape genetics because of the non-independence of explanatory variables (Legendre and Fortin, 2010; Van Strien et al., 2012). The new approaches suggested to overcome these issues relate to comparison of distance matrices (Van Strien et al., 2012) that is not suitable for our analysis based on Ar. An analysis of variance was carried out on each of our four highest ranked models to test whether models significantly explained the variation in Ar across sites. Finally, hierarchical linear regression was carried out to determine which variables significantly contributed to the predictive power of each model. We considered a model to be meaningful if it had a positive or negative regression coefficient, significantly accounted for variation in Ar (analysis of variance P<0.05) and each variable entered had a significant contribution to the predictive power of the model (P<0.05).

Results

Data quality and low levels of genetic variation

Missing data accounted for only 0.3% of the genotypes and these were spread evenly over populations and loci. No significant linkage disequilibrium was identified after sequential Bonferroni correction. There was no evidence of scoring errors. MICRO-CHECKER suggested null alleles in 8 loci, but the results were not consistent across the 12 population samples. Gama03 was the only locus with evidence for null alleles in two populations (SCC and TOO). All analyses were run with and without Gama 03 and produced similar results. Deviations from Hardy–Weinberg equilibrium were evident in 7 out of 17 populations; these were mostly because of excess of homozygotes (Table 1).

Overall genetic diversity was strikingly low with an average Ar of only 2.2 alleles (Table 1). Average number of alleles within populations was also low, ranging from 1.92 (BRN) to 3.67 (TOO). There were significant differences in Ar (P<0.001) and HE (P<0.001) between populations. Allelic richness ranged from 1.57 (KPC) to 2.9 (SCP) and HE ranged from 0.25 (SCC) to 0.66 (CRT and SCP). In particular, populations BRN and KPC showed lowest diversity, including lowest percentage of polymorphic loci (50% and 58%, respectively). High and mostly significant FIS values, suggestive of inbreeding, were found in most populations (overall average FIS=0.23; Table 1).

Population structure

The global Fst was very high (Fst =0.45, P<0.001), indicating marked population differentiation for G. marmoratus across the MDB. All pairwise population comparisons were significant, with Fst ranging from 0.13 to 0.67 (Supplementary Table S3).

AMOVA identified basin-level differentiation only between the Upper Darling and lower Murray (P=0.019), suggesting that populations are not structured mainly as a consequence of regional differences. On the other hand, both the catchment- and the regional-level AMOVAs identified that the largest amount of the total variation was accounted for by differences between populations, both between catchments and within streams within catchments (P<0.01). This indicates that populations within catchments are more related and that catchment hierarchy explains a significant level of differentiation (Table 2).

Table 2 Analysis of molecular variance (AMOVA) for Gadopsis marmoratus in the Murray–Darling Basin
Full size table

In terms of IBD results, the Mantel tests using pairwise population Fst and pairwise stream distances matrices identified no correlation between the two variables and therefore did not support a pattern of IBD. This result was consistent across all regional scales tested (Supplementary Table S4).

STRUCTURE analysis indicated a value for K of 11 populations (Figure 2). Populations were well defined, with either low or no admixture between clusters. Only one of the inferred clusters was represented by more than one sampling site (streams NAN and TOO); these are the two closest sites in our study, separated by only 11 km. Samples from these sites were pooled together for the BAYESSASS migration analysis. A few admixed individuals are also apparent in the STRUCTURE graph (Figure 2), but these were generally found between sites with smaller sample sizes (for example, AVL and MAR) that probably relates to our reduced ability to estimate allele frequency for these populations (Evanno et al., 2005). All potential cases of admixture were successfully assigned to their sampled populations via assignment tests. A high proportion of individuals in all 11 populations (97%; range of 89–100%) had multilocus genotypes that assigned to their sampled population based on GENECLASS results (Supplementary Table S5). Fish from the undersampled locations (n⩽5) were not assigned to any of the 11 population clusters (Supplementary Table S5).

Figure 2
figure 2

STRUCTURE plot depicting 11 inferred population clusters for G. marmoratus in the MDB. Only sites with ⩾8 sampled individuals were included. Each bar represents an individual and the color represents the proportion of membership to each cluster. The abbreviation of the source population follows Table 1 and black lines delineate the boundaries of the source populations.

Full size image

BAYESASS did not identify significant migration between groups at any of the hierarchical levels (that is, within catchments, between catchments and between regions; Table 3). The majority of tests indicated <2% of migrant ancestry from other populations, with confidence intervals with a lower bound overlapping zero. The only exception to this is the unidirectional gene flow of 17% from MAR to SCR, but this could be related to MAR having a small sample (n=8) (Wilson and Rannala, 2003). Overall, migration analyses strongly corroborate the results of population differentiation based on various analyses and the individual assignment, indicating that MDB river blackfish are composed of a minimum of 11 genetically distinct and currently isolated populations.

Table 3 Recent migration rates among Gadopsis marmoratus populations in the Murray–Darling Basin estimated in BAYESASS
Full size table

Influence of genetic differentiation on genetic diversity

As expected, population estimates of genetic diversity (for example, Ar and HE) were significantly correlated (P<0.001). A strong negative relationship was found between genetic diversity and population genetic differentiation (r2=0.923, P<0.001; Supplementary Figure S1). This indicates that the most genetically divergent populations have lower levels of genetic diversity.

Landscape genetics

The model that best explained the impact of environmental factors on genetic diversity included both average river slope and land use (r2= 0.312; Table 4). However, hierarchical regression showed that land use did not add significantly to the predictive power of the model after controlling for the effect of river slope. The third ranked model included only river slope, which accounted for 56.3% of the variation in genetic diversity, and was the only other model that accounted for a significant amount of variation in genetic diversity (r2=0.194; Table 4).

Table 4 Models of environmental variables potentially accounting for the distribution of genetic diversity (Ar) of Gadopsis marmoratus in the Murray–Darling Basin
Full size table

Discussion

Poor dispersal species act as important conservative baselines for biodiversity management, as they are more likely than vagile species to be negatively affected by the ecological and evolutionary processes accompanying habitat fragmentation. Large-scale habitat fragmentation, such as that observed in the MDB (Leblanc et al., 2012), initially results in increasing isolation of demes and, second, in a decrease of the total amount of habitat available for each deme (Simberloff, 1998). Understanding the genetic consequences of these impacts, both at the deme and species range levels, enables process-oriented strategies for conservation management to be implemented at appropriate spatial scales (Young and Clarke, 2000). Here we have used the poorly dispersive and threatened G. marmoratus as a surrogate indicator system for assessing the effects of fragmentation in biotas from highly modified river basins.

Our results indicate that G. marmoratus displays all the hallmarks of severe habitat fragmentation, namely scattered, small, demographically isolated, and genetically depauperate populations. Such findings make it very difficult to prioritize conservation actions throughout the basin based on patterns inferred from genetic information. Nonetheless, by combining extensive and hierarchically nested population sampling (that is, from the scale of stream, catchment, to region within an expansive river basin) with expert opinion and field observations, our study identified key associations between landscape features and genetic variation that would otherwise not be detectable in conventional population genetic surveys. As discussed below, these results can be integrated into an objective framework for prioritizing management strategies aimed at maintaining population persistence in poor dispersal species affected by habitat fragmentation.

Remarkably strong population structure in a recently fragmented ecosystem

G. marmoratus is represented by a strongly structured lineage in the MDB (for example, global Fst of 0.45). The inferred levels of population differentiation can be considered particularly high as they are based on hypervariable markers that are expected to show lower differentiation than markers with lower within-population heterozygosity (Hedrick, 1999). Irrespective of whether populations came from within the same or separate catchments, they are all essentially isolated, with no contemporary gene flow between them (Figure 2 and Supplementary Tables S3–S6). Previous mitochondrial DNA and allozyme analyses covering the entire range of Gadopsis (including many of the populations examined herein) confirm that our samples represent a single candidate species and a single ESU (that is, evolutionarily significant unit) in the MDB (Hammer et al., 2014), and thus exclude the possibility that cryptic lineages are driving this result within our data set. Many native MDB species have poor or moderate dispersal potential and a relatively high degree of habitat specialization. These characteristics, in combination with human-driven fragmentation and disturbance of the basin (Lintermans, 2007; Hammer et al., 2013), account for the high population structure and isolation typically detected in population genetic studies of MDB fishes (see, for example, Faulks et al., 2011; Brauer et al., 2013; Attard et al., 2016; Cole et al., 2016; Sasaki et al., 2016). Therefore, our microsatellite results support the proposal that G. marmoratus represents a valuable surrogate indicator species of how fragmentation might affect co-occurring low dispersal species in affected riverine ecosystems.

It is not surprising that no signal of IBD was detected at the regional scale given the low dispersal potential of G. marmoratus and large intervening sections of apparently unsuitable modified lowland river habitat (Khan et al., 2004; Koster and Crook, 2008). In highly fragmented environments this pattern can be disrupted by restricted gene flow between demes (Broquet et al., 2006). Patterns of IBD may still exist within streams or within catchments where continuous distribution of the species remains in longer river reaches (for example, the Ovens River (Lintermans, 2007)). It is expected that genetic similarity between individuals persists only over the spatial scale of a stream, as predicted for a poor disperser (Hughes et al., 2012).

It is widely understood that both contemporary and historical gene flow can influence the architecture of genetic variation in populations, and distinguishing between the two can be difficult (Poissant et al., 2005). It is likely that both the recent habitat fragmentation (past 200 years post European settlement) and species life history account for the strong contemporary population structure observed here. This hypothesis is supported by differences in the estimated widespread historic distribution compared with the current highly restricted range of G. marmoratus in the MDB (Lintermans, 2007) and by the shallow phylogeography displayed by this species in the MDB (Hammer et al., 2014). Our data also show that an unperturbed reference population from the lower Murray (Tookayerta Creek), which is by far the strongest remnant population and has experienced less decline (Hammer, 2004), had higher genetic variation than most other populations that experienced human impacts (Table 1). These findings are also consistent with coalescent-based and genome-wide analyses of a co-occurring poor dispersal species in the MDB, the southern pygmy perch Nannoperca australis. Here, isolation and basin-wide demographic declines observed for several populations are better explained by human-driven disturbances that post-date European settlement than by the historical effects of genetic drift (Attard et al., 2016; Brauer et al., 2016; Cole et al., 2016). Such analyses would probably prove useful to better delineate factors that have shaped the history of isolation of G. marmoratus populations.

The inferred overall pattern led us to reject the stream hierarchy model for G. marmoratus in the MDB and to propose either the Death Valley model or a nonequilibrium model of population structure (Meffe and Vrijenhoek, 1988; Hughes et al., 2009). The Death Valley model predicts that remnant populations across the MBD landscape have a high level of between-population structure, but no significant basin or catchment structure, no IBD, no gene flow and small local population size, perhaps driven by naturally chaotic and variable patterns of stream connectivity exacerbated by human-driven change (Walker et al., 1995). The nonequilibrium model is defined by a high level of local extinction with infrequent recolonization events (Nei and Chakravarti, 1977), indicative of a species decline in fragmented habitats. The proposed models further emphasize the extinction risk for G. marmoratus.

Conserving and managing population diversity and uniqueness across the basin

The current fragmented nature of G. marmoratus populations makes each isolated population extremely vulnerable to stochastic events (Fagan et al., 2005). This is particularly true given its low levels of genetic variation relative to other freshwater fishes (DeWoody and Avise, 2000) and values of inbreeding coefficient reported for populations across the MDB (Table 1). Conservation implications from population genetic data have traditionally focused on management units (MUs) that are demes with significantly divergent allele frequencies that occur within an ESU (Moritz, 1994). A range of analytical methods consistently indicated that each G. marmoratus population we assessed in the MDB currently represents an isolated demographic unit. Hence, we identified a minimum of 11 potential MUs, with the likely presence of others because of fish from undersampled locations not being assigned to any of the population clusters (MUs listed in Supplementary Table S1). Assessing the true number of MUs justifies further sampling effort across catchments and within streams.

The contemporary genetic architecture of G. marmoratus in the MDB makes it difficult to prioritize conservation efforts across this vast river basin. Conservation management actions should normally focus on divergent (and potentially locally adapted) populations with high levels of genetic diversity to give the species the best chance of persistence into the future (Hughey et al., 2003; Coleman et al., 2013; Brauer et al., 2016). We found a strong negative relationship between genetic diversity and population differentiation in this species. The highly divergent populations showed the lowest levels of genetic diversity (Supplementary Figure S1), possibly reflecting the influence of genetic drift in currently small and isolated populations. If management efforts were to be prioritized based on differentiation alone, the species diversity as a whole could be eroded over time, leading to an increased chance of extinction (Coleman et al., 2013). For example, focusing efforts on Browns Creek (BRN) and King Parrot Creek (KPC) would conserve the most differentiated populations, but both show relatively few alleles (1.92 and 2.08, respectively) and no or few private alleles (0.00 and 0.17, respectively). These MUs contain a subset of alleles found in other populations and their high differentiation is likely the result of recent population declines and small effective population size. On the other hand, limited conservation resources could be more efficiently directed toward MUs such as Hughes Creek (HUG) and the Tookayerta Creek (TOO/NAN), as they are genetically divergent (Supplementary Table S3) while maintaining relatively high levels of diversity (Ar=2.38 and 2.30, respectively) and a high proportion of private alleles (0.50 and 0.25, respectively). Based on this concept, in addition to HUG and TOO, the other high priority MUs that should be targeted by management actions (described below) are Molong Creek (MOL) and Catherines Creek (CAT) (Figure 3).

Figure 3
figure 3

Population genetic diversity (Ar) plotted against genetic uniqueness (population-specific Fst) for G. marmoratus in the MDB. The correlation was highly significant (P<0.001).

Full size image

The diverse and divergent MUs identified here represent appropriate conservation targets for a variety of management agencies across multiple jurisdictions of the MDB. A short-term ex situ management action recently conducted with endangered fishes in the MDB is genetically informed captive breeding followed by reintroduction to the wild (Hammer et al., 2013; Attard et al., 2016). This was implemented for four endangered species in the lower Murray and was particularly successful with two species of pygmy perches. Genetic diversity in pygmy perches was maintained after one generation of captive breeding, a few thousand offspring were reintroduced to the wild after some habitat recovery, and subsequent recapturing events demonstrated not only survival but also recruitment of captive-born individuals (Attard et al., 2016). However, the same program had very limited success in breeding G. marmoratus in captivity as a larger and territorial species. In spite of using controlled trials in a dedicated aquaculture facility, most spawning attempts proved unsuccessful and only eight captive-reared juveniles were produced (Westergaard and Ye, 2010), emphasizing the need for further research into G. marmoratus captive breeding (Hammer et al., 2013). Because of these issues, we propose that conservation agencies should instead focus on maintaining and restoring core habitat (for example, environmental water allocation, riparian restoration and enhancement of water quality; sensu Hammer et al., 2013) and that these in situ management actions should be done at the level of catchments.

A key conservation management action for species with fragmented distributions, reduced genetic diversity and lowered adaptive potential is to conduct assisted gene flow by moving individuals between populations (Tallmon et al., 2004; Frankham, 2015; Whiteley et al., 2015). Translocation of genetically diverse individuals between recently isolated populations is a form of genetic rescue, a management tool that can increase neutral and adaptive genetic diversity, population fitness and evolutionary potential (Frankham, 2015; Whiteley et al., 2015). Genetic rescue offsets the detrimental effects of drift and inbreeding in small populations, provides population resilience to environmental change and halts biodiversity loss (Weeks et al., 2011; Frankham, 2015; Whiteley et al., 2015). Population translocations remain a controversial topic among conservation biologists, in part because of the potentially negative effects of outcrossing (that is, outbreeding depression). However, recent meta-analyses have consistently indicated that concerns of outcrossing small populations are not justified (provided an appropriate systematic framework is established, that is, excluding mixing of species and ESUs) and have called for a much broader use of assisted gene flow to genetically rescue small inbred populations (Frankham, 2015; Whiteley et al., 2015). These studies also showed that outbreeding depression generally increases with genetic, geographic and environmental distance (Whiteley et al., 2015) and confirmed the effectiveness of decision trees (for example, same species, no fixed chromosomal differences, adapted to similar environments and gene flow within the past 500 years; Frankham et al., 2011) in assessing the risk of outbreeding depression (Frankham, 2015).

G. marmoratus is represented in the MDB by a single genealogical lineage (that is, one ESU) that shows very low divergence across the basin (mean within lineage mitochondrial DNA distance of 0.45%; Hammer et al., 2014). Its phylogeography is more indicative of recent range expansion (for example, late Pleistocene) across the MDB rather than relictual persistence (Hammer et al., 2014). Such a pattern is shared by many other widespread MDB fishes (see, for example, Faulks et al., 2010; Unmack et al., 2013) and appears related to the much greater flow of ancestral MDB rivers during the Pleistocene (Pels, 1964). As a specific conservation action, we propose that assisted gene flow should be performed within MDB G. marmoratus. This should be at regular intervals and involve moving individuals from four of the five targeted G. marmoratus MUs (listed above) to populations with lower genetic diversity found in catchments exposed to similar environmental conditions. Movement within the three major MDB regions may be most practical considering the potentially different climates and disease histories in these widely separated areas. In contrast, the Tookayerta Creek Catchment MU (TOO/NAN) is a genetically divergent and heterogeneous population occurring in atypical habitat (acidic peat swamp system), suggestive of a separate genetic and evolutionary trajectory (Hammer et al., 2014), and this is an example where translocation into or out of an MU should be avoided.

Pools promote diversity

Our model selection approach, hierarchical regression and hierarchical partitioning (Table 4 and Supplementary Figure S2) identified average river slope as the only major factor influencing genetic diversity in G. marmoratus in the MDB. The observed correlation is consistent with our initial hypothesis that a low average slope would increase the area of available habitat (deep slow-flowing pools), therefore supporting larger populations and higher genetic diversity. Our results may also reflect the negative influence of flow velocity on population size, given that G. marmoratus cannot sustain long periods of intense swimming (Dobson and Baldwin, 1981). The influence of river slope is supported by previous habitat use, mark recapture and radio telemetry studies that concluded that G. marmoratus have a strong affinity to slow-flowing pools and only use riffle habitat occasionally for short periods of time (Khan et al., 2004). The identification of a single proxy for genetic diversity (slope) does not necessarily imply that other habitat features are not influencing population size. Most of the few studies that identified relationships between habitat and genetic diversity had access to fine-scale habitat data and included multiple samples within a predefined habitat type (see, for example, Cena et al., 2006; Balkenhol et al., 2013). The vast geographic range but threatened status of G. marmoratus precluded replication of sampling within habitat types. In addition, detailed spatial mapping of riparian vegetation, small in-stream structures and the presence of introduced species are not available across the entire MDB. These variables are likely to influence the dynamics and demography of a range of native freshwater species (Jackson and Williams, 1980; Growns et al., 1998; Erskine and Webb, 2003). These issues, in combination with the relatively low levels of variation found in our microsatellite data set and the small size of some population samples, could potentially affect our ability to detect landscape genetic patterns.

Conclusions

We have shown G. marmoratus to be an ideal surrogate indicator species for assessing landscape fragmentation, and in the process also gleaned key data with conservation implications for this poor dispersal and declining species. Our proposals of genetic rescue relate to building resilience in the interim to landscape restoration and in the face of climate change. Restorative actions are clearly required and can build upon the spatial management framework identified herein. Indeed, ongoing temporal genetic monitoring of populations could provide empirical and applied insights on both the effectiveness and scale of influence of restoration actions such as flow restoration, habitat improvements and catered translocations, as well as the real-time and lag effects from past and present threatening processes such as fragmentation (small populations) and climate change.

Data archiving

Genotype data have been submitted to Dryad: doi:10.5061/dryad.d5842.