The role of anthropogenic vs. natural in-stream structures in determining connectivity and genetic diversity in an endangered freshwater fish, Macquarie perch (Macquaria australasica).

Habitat fragmentation is one of the leading causes of population declines, threatening ecosystems worldwide. Freshwater taxa may be particularly sensitive to habitat loss as connectivity between suitable patches of habitat is restricted not only by the natural stream network but also by anthropogenic factors. Using a landscape genetics approach, we assessed the impact of habitat availability on population genetic diversity and connectivity of an endangered Australian freshwater fish Macquarie perch, Macquaria australasica (Percichthyidae). The relative contribution of anthropogenic versus natural in-stream habitat structures in shaping genetic structure and diversity in M. australasica was quite striking. Genetic diversity was significantly higher in locations with a higher river slope, a correlate of the species preferred habitat - riffles. On the other hand, barriers degrade preferred habitat and impede dispersal, contributing to the degree of genetic differentiation among populations. Our results highlight the importance of landscape genetics to understanding the environmental factors affecting freshwater fish populations and the potential practical application of this approach to conservation management of other freshwater organisms.

Introduction

At large spatial and temporal scales, such as across landscapes and geological time, habitat fragmentation is a natural process that promotes the diversification of ecological niches and the evolution of species. However, contemporary habitat fragmentation is largely induced by anthropogenic activities and usually occurs at a much faster rate than the natural process. Indeed, habitat fragmentation is now recognized as one of the leading threats to ecosystems worldwide (Kingsford et al. 2009). Freshwater environments are naturally fragmented by mountain ranges and coastlines that delineate basin boundaries and are also structured as hierarchical riverine networks (Fausch et al. 2002; Magalhaes et al. 2002; Boys and Thoms 2006) with features such as waterfalls and cascades acting as natural barriers. Some of the most disruptive anthropogenic activities in freshwater environments are the construction of impoundments and agricultural or industrial practices (Cadwallader 1978; Allan and Flecker 1993; Maitland 1995). For example, the construction of impoundments can affect habitat structure, change macro-invertebrate community composition, alter natural flow regimes and impede fish dispersal (Petts 1984; Gehrke et al. 2002). Changes in water quality can also result from impoundment construction but are often a consequence of catchment land use such as clearing for agriculture and mining activities (Dudgeon et al. 2006; Growns et al. 2009). These landscape-level processes can impact immediate in-stream habitats as well as have cascading effects at larger scales within a river basin.

Similarly, habitat fragmentation can have both direct and indirect consequences at the population level. Habitat loss and fragmentation are known to result in reduced population sizes and reduced demographic and genetic connectivity (Ewers and Didham 2006; Lowe and Allendorf 2010). The latter takes place because increased isolation and novel ecological boundaries can impact on the relative contribution of net immigration to demographic connectivity – total population recruitment and also genetic connectivity – the effective number of dispersers (Lowe and Allendorf 2010). In turn, fragmented populations are supposed to experience reduced levels of genetic diversity, population growth rate and population viability (Frankham 1995; Ewers and Didham 2006) and increased risk of local extinction (Spielman et al. 2004). Ongoing habitat fragmentation also impacts negatively on evolutionary resilience – the ability of populations to persist in their current state and to undergo evolutionary adaptation in response to a changing environment (sensuBlanchet et al. 2010; Sgro et al. 2010; Ewers and Didham 2006).

The eventual impact of these demographic and genetic changes on a species can be influenced by the species biology, particularly dispersal ability, life-history characteristics and habitat requirements. Dispersal potential enables species to respond to changes in the availability of preferred habitat and can be crucial in determining the effect of habitat fragmentation on genetic diversity (Raeymakers et al. 2008; Watanabe et al. 2010). A recent review revealed a significant negative association between the genetic diversity of fish populations and human disturbance (DiBattista 2008), and there are several examples in freshwater species (e.g. Belfiore and Anderson 2001). The majority of studies of freshwater taxa have focused on the impact of fragmentation on key migratory species such as salmon; however, there is increasing recognition of whole community effects (Gehrke et al. 2002). For species that are naturally more sedentary, the genetic consequences of habitat fragmentation such as reduced diversity could be accentuated (Dehais et al. 2010; McCarney et al. 2010).

A valuable approach for exploring the influence of contemporary environmental factors on genetic diversity is landscape genetics, an area of study that integrates population genetics, environmental data and spatial statistics (Manel et al. 2003; Storfer et al. 2007; Holderegger and Wagner 2008). Landscape genetics can provide a greater understanding of the interactions between species, ecosystems and anthropogenic activities and can contribute to conservation management of endangered species and the environments in which they reside. Here, we use a landscape genetics approach to assess the impact of habitat availability on population genetic diversity and structure of an endangered (Australian Environment Protection and Biodiversity Conservation Act 1999) Australian freshwater fish, the Macquarie perch Macquaria australasica (Percichthyidae). This species is a habitat specialist and remnant populations are found in rocky gorge country with an abundance of riffles and where water quality is good (Bruce et al., 2008; Gilligan et al. 2010). Adults are thought to undergo small-scale migrations during the breeding season to spawn in riffle beds (see Methods for more details on species biology). These characteristics make M. australasica an ideal species with which to test the effect of habitat fragmentation on organisms with limited dispersal behaviour and specific habitat requirements. Specifically, we evaluated the roles of natural habitats (i.e. riffles) as well as anthropogenic in-stream structures (i.e. weirs) in driving levels of genetic diversity, population structure and contemporary patterns of genetic connectivity. We predicted that the availability of optimal habitat (riffles) would promote genetic diversity and conversely that habitat fragmentation (anthropogenic structures) would decrease genetic diversity and restrict connectivity. We anticipate that the outcomes of this study will contribute towards the conservation management of this endangered species. In addition, the results should further our understanding of the genetic consequences of habitat loss and fragmentation on other freshwater organisms that have specific habitat requirements and show limited dispersal behaviour.

Methods

Species distribution, status and biology

Remnant populations of M. australasica are restricted to the headwaters of the Lachlan, Murrumbidgee, Murray, Kiewa, Ovens, Goulburn-Broken and Campaspe Rivers in the Murray-Darling Basin (MDB) to the west of the Great Dividing Range (GDR), as well as the Hawkesbury-Nepean (HN), Georges River and Shoalhaven basins on the east of the GDR (Industry and Investment NSW Freshwater Fish Research Database; Lintermans 2007) (Fig. 1). Previous genetic and morphological studies have provided strong support for the differentiation of lineages on either side of the GDR (Dufty 1986; Faulks et al. 2010a). During the early 1900s, translocation occurred from the MDB into the Yarra River (Victoria), Mongarlowe River (New South Wales) and Cataract Dam (New South Wales) (Cadwallader 1981). However, interbasin translocations are no longer undertaken. Although both lineages co-occur in the Cataract River (because of the translocation of individuals), the extent of hybridization between them is unknown and the taxonomic status of these two groups is yet to be established (Faulks et al. 2010a). Despite the species being protected throughout its natural range, the natural population in the Shoalhaven catchment (Kangaroo River) declined rapidly during the late 1990s, and despite intensive sampling effort (118 sampling occasions), no individuals have been observed in the wild since 1998. Monitoring of M. australasica populations elsewhere indicates that the species is absent from up to 95% of its previous distribution (Bruce et al. 2008; Davies et al. 2008). Loss of habitat and the construction of impoundments are likely to be largely responsible for these declines (Lintermans 2007; Bruce et al. 2008). Other causes of population decline include increased competition from native species following stocking, introduced species (Perca fluviatilis and salmonids), overfishing and susceptibility to epizootic haematopoietic necrosis virus (EHNV) (Cadwallader 1978; McDowall, 1996; Allen et al. 2002; Lintermans 2007).

Figure 1

Map of Australia indicating the study region. The first inset shows details of sampling locations of Macquaria australasica. The second inset provides further details of sampling locations in the Hawkesbury-Nepean Basin. Dotted line indicates the Great Dividing Range (GDR).

Macquaria australasica adults from the MDB reach an average length of 350 mm, and sexual maturity is reached at 2 years for males and 3 years for females. Details of the age, growth and maturity of coastal populations are still unknown. Reproductive conditioning occurs during late winter/spring, and adults undergo small-scale migrations from pools to riffles to spawn in spring/early summer when water temperature reaches approximately 15°C (Wharton 1968; Cadwallader and Douglas 1986). Eggs and larvae are demersal and remain in the riffles and the extent of dispersal in juvenile stages is also thought to be limited, resulting in localized recruitment (Cadwallader and Rogan 1977). Even at this small scale, anthropogenic structures have the potential to interfere with the species reproductive cycle by restricting movement and habitat availability.

Sample collection and DNA extraction

A total of 299 samples from 26 locations (n = 5–40) were obtained from across the species range (Fig. 1) over several years (2003–2008). Despite the potential for temporal instability in population structure (e.g. McCairns and Bernatchez, 2008), in this case, localized recruitment and restricted dispersal mean that it is unlikely to have influence on our results. Fish were collected using nondestructive methods such as electrofishing and bait trapping. Small pieces of caudal fin tissue were taken, placed in 100% ethanol and stored at −20°C in the laboratory. Total DNA was extracted from the tissue using a salting-out method. We modified the method of Sunnucks and Hales (1996) by doubling the volume of TNES (600 μL) in the digestion step and saturated salt solution (170 μL) in the precipitation step. One sample from the South Australian Museum and a few other samples that yielded low-quality DNA using the salting-out method were extracted using a QIAGEN DNeasy kit (Qiagen, Hilden, Germany).

Microsatellite genotyping

Eight microsatellite DNA loci were amplified using PCR primers originally developed for Australian Bass (Macquaria novemaculeata) [AB009 (Schwartz et al. 2005)] and Murray Cod (Macculochella peelii) [Mpe3.B11, Mpe2.E01, Mpe1.F01, Mpe2.F07, Mpe3.G04, Mpe3.G12, Mpe1.H04 (Rourke et al. 2007)]. The forward primer for each primer pair was incorporated with a fluorescent-labelled M13 tag (Schuelke 2000). Loci were amplified and fluorescent dyes incorporated in either single (AB009 VIC, Mpe3.B11 NED, Mpe2.E01 NED, Mpe3.G12 VIC) or multiplex (Mpe3.G04 and Mpe1.H04 PET, Mpe1.F01 and Mpe2.F07 FAM) reactions. PCRs consisted of 1 μL template DNA, 2.5 mmol MgCl2, 60 mm KCl, 12 mm Tris–HCl pH 9.0, 0.12% Triton-X, 0.4 mm each dNTP, 0.78 pmol M13 forward primer, 3.8 pmol reverse primer, 3.8 pmol M13 fluorescent dye, 0.5 μg bovine serum albumin, 0.25 U Taq DNA polymerase (Promega, Madison, WI, USA) and dH2O to a final volume of 10 μL. PCR conditions were as described in Schwartz et al. (2005). PCR products for the eight loci (3–4 μL) for each sample were pooled, and 1 μL of this product was mixed with 0.1 μL LIZ500 ladder and 9.9 μL HiDi. Products were then screened on an ABI PRISM 377 with Genescan software (Applied Biosystems, Carlsbad, CA, USA) at the Macquarie University DNA Analysis Facility. Genotypes were visually inspected and scored using GeneMapper (Applied Biosystems). In cases of poor amplification, re-runs were conducted with a single locus per reaction and without pooling. Genotypes were checked for scoring errors because of null alleles, stutter and large allele dropout using Microchecker 2.2.3 (van Oosterhout et al. 2004).

Genetic diversity

Fisher's exact test of linkage disequilibrium between all pairs of loci and conformance to Hardy–Weinberg equilibrium in each population across all loci were tested using Genepop 3.4 (Raymond and Rousset 1995). Significance levels were Bonferroni-corrected to address the risk of increased Type 1 error associated with multiple tests (Rice 1989) and Benjamini-Yekutieli (B-Y)-corrected to balance the risks of Type 1 and Type 2 errors (Narum 2006). The expected (HE) and observed (HO) heterozygosities were calculated in Arlequin 3.1 (Excoffier et al. 2006), and the corrected allelic richness (AR) and inbreeding coefficient (FIS) per population were calculated in FSTAT 2.9.3.2 (Goudet 2002). To evaluate whether populations may have experienced recent bottlenecks, we calculated the M ratio [Arlequin 3.1 (Garza and Williamson 2001)] and levels of heterozygosity [Bottleneck 1.2.02 (Piry et al. 1999)]. The M ratio is the mean ratio of the number of alleles compared to the range in allele size, and values smaller than 0.68 indicate that the population has experienced a recent and severe reduction in population size (Garza and Williamson 2001). Hybridization can invalidate the M ratio; therefore the Cataract River population, where translocation of distinct lineages has occurred, was not tested. Excess heterozygosity, owing to a faster rate of the number of alleles lost compared to the reduction in gene diversity, was assessed using both the stepwise (SMM) and two-phased models of mutation (TPM) and a Wilcoxon test in Bottleneck (Piry et al. 1999).

Population genetic structure

Population differentiation was evaluated using global and pairwise FST tests, and levels of population subdivision were assessed (AMOVA) among basins and among locations within basins. Weir and Cockerhams (Weir and Cockerham 1984) FST values and significance levels were calculated in Arlequin 3.1 (Excoffier et al. 2006). We used FST rather than RST because of our relatively small number of loci (<20) and recently diverged populations (Gaggiotti et al. 1999), as shown in our mtDNA phylogeographical study of this species (Faulks et al. 2010a).

Preliminary data analyses using both STRUCTURE (nonspatial) (Pritchard et al. 2000) and Geneland (spatial) (Guillot et al. 2005) to estimate population structure obtained similar results. Geneland was chosen for final analyses as the incorporation of spatial information provides a more realistic representation of the study system, optimizing the ability to determine population structure (sensuDionne et al. 2008; Latch et al. 2008). All translocated populations were removed from this analysis (Yarra River, Mongarlowe River, Cataract Dam and River and Loddon Ck). Geneland uses a Markov chain Monte Carlo (MCMC) approach to identify genetic spatial discontinuities. Five replicate runs were performed, each with 1 × 106 MCMC iterations, the Dirichlet model of allele frequencies, no uncertainty in the spatial coordinates, maximum rate of Poisson process 100 and the maximum number of nuclei in the Poisson-Voronoi tessellation set to 200. Plots of MCMC runs were assessed for mixing. The most probable number of population clusters (K) was the same for each run. The final output displays maps of posterior probabilities of locations belonging to a particular genetic cluster (K) or population. Although these clusters (classes) may include areas of suboptimal habitat, the important consideration is that a class contains populations that are grouped together with a higher probability than with populations from any other class.

Gene flow

Estimates of recent gene flow across the GDR as well as among populations within the MDB and east coast were assessed in BayesAss+ 1.3 (Wilson and Rannala 2003). The program was run for 3 × 106 iterations including a burn-in of 1 × 106 iterations. Delta values for allele frequencies, inbreeding coefficients and migration rate were set at 0.05, 0.10 and 0.02 for the entire data set; 0.45, 0.45 and 0.25 in the MDB; and 0.35, 0.35 and 0.15 on the east coast. These values achieved the recommended acceptance rates of changes of 40–60% (Wilson and Rannala 2003). Convergence was assessed by plotting the cumulative log likelihoods of the iterations. Estimates were considered ‘real’ if they were consistent in at least seven of 10 replicate runs, and the run with the best acceptance rates and convergence was chosen for the parameter estimates.

Landscape genetics

We used a general linear model selection approach (e.g. Banks et al. 2007) to identify the environmental variables influencing genetic diversity [allelic richness (AR)]. AR was used as the measure of diversity as it is a sensitive indicator of changes in population demography (Nei et al. 1975). The following environmental variables were considered: distance of the site from the source (headwater), river slope, total number of anthropogenic barriers within the entire reach (including all upstream tributaries and downstream to the junction with the next larger waterway in the stream hierarchy), riparian land use immediately surrounding the location (ranked from 1 to 4 in increasing order of disturbance) and elevation and latitude of the sampling location. Distance from the source was measured in ArcMap. River slope was calculated as the change in altitude in m/km from 5 km upstream to 5 km downstream of each sampling location. This measure was used as a surrogate of potential riffle habitat within the reach. Although river slope may also indicate flow velocity and valley shape, this approach is considered appropriate in the absence of detailed habitat mapping across the species distribution (Jowett 1993; Thompson et al. 2006). The number of barriers was calculated using the I&I NSW database (I&I NSW, 2008). These data include all structures of all ages, ranging from culverts on road crossings to weirs and large dams on impoundments (I&I NSW, 2008). Unfortunately, equivalent data from catchments in Victoria were unavailable, so data from those locations were removed from this analysis. Land use was inferred from Global Map 2001 (Geoscience Australia) and was categorized as follows: 1 = Forests, 2 = Grasslands, 3 = Water, 4 = Agriculture. For the model selection analyses, all variables were standardized to a mean of zero and a standard deviation of one. The assumptions of general linear models, including normality and colinearity of variables, were assessed prior to performing the modelling analyses. Only two of the variables were slightly skewed to lower values: river slope and distance from the headwater source. In addition, only five of the 15 pairwise comparisons were significantly correlated, with none of the variables consistently correlated with all others (Appendix 2). Therefore, all variables were retained in the modelling analyses.

The best models were selected using Akaike's Information Criterion for small sample sizes (AICc) (Burnham and Anderson 2002) and the independent effect of each variable to the model estimated using hierarchical partitioning (Walsh and MacNally 2008). In addition, we conducted a re-sampling-based regression analysis to estimate the effect of our environmental variables on AR while controlling for unequal sample sizes among sites. In each of 10 000 iterations, we re-sampled an equal number of individuals (5) from each population, recalculated AR over loci within a population and fitted a simple linear regression model of the best environmental variable to the AR data. We obtained the distribution (95% CI) of the regression coefficient and considered our conclusions robust if >97.5% of the coefficients were either consistently positive or consistently negative.

Finally, we performed simple and partial Mantel tests in FSTAT 2.9.3.2 (Goudet 2002) to determine the effect of anthropogenic barriers and riverine distance on population differentiation (FST). These tests were conducted on the coastal populations only, because of the small number of sites for which data were available in the MDB. Riverine distances between pairs of populations were calculated using ArcMap, and the number of barriers between populations was calculated using the I&I NSW database (I&I NSW, 2008). To investigate the effect of variable samples sizes, we performed two tests, the first with FST values calculated from the entire data set and the second with FST values calculated from a single random selection of five individuals from each population. Each test involved 2000 permutations.

Results

We amplified eight polymorphic loci with no consistent evidence for stutter, large allele dropout or linkage disequilibrium. Some loci showed evidence of null alleles (E01, G04 and G12). However, these were not consistent across populations, and when analyses were run without E01, G04 and G12, similar results were obtained. Thirteen populations (of 26) deviated significantly from Hardy–Weinberg equilibrium when averaged over all loci (15 of 26 using the B-Y correction) (Table 1). These deviations are most likely due to local inbreeding effects as shown by the high inbreeding coefficients (FIS > 0, Table 1) and the lack of loci specific patterns of disequilibrium. Diversity indices for each population are also shown in Table 1. Overall, the species showed very low genetic diversity. Mean allelic richness was only 1.4, and mean observed heterozygosity was 0.28. The M ratio results for TPM and SMM were the same and indicated that most populations have undergone recent bottlenecks. However, results from the Wilcoxon test of excess heterozygosity were more conservative and only eight populations showed evidence of bottlenecks (Table 1): Broken, Buffalo, Cotter, Bowen, Cataract Dam, Cordeaux, Loddon and Wheeny. Four of these populations are in or above major impoundments.

Table 1

Genetic diversity indices characterized at eight microsatellite loci for Macquaria australasica populations

Location	N	Mean AR	Mean HO	Mean HE	FIS	P value	M ratio
MDB
AB	40	1.5	0.41	0.50	0.19	<0.001*	0.56
BR	5	1.45	0.28	0.45	0.42	0.12	0.54*
BUFF	6	1.3	0.30	0.30	0.00	0.69	0.62*
CO	10	1.28	0.12	0.28	0.60	<0.001*	0.60*
DD	7	1.46	0.33	0.46	0.30	0.14	0.52
HU	9	1.42	0.33	0.42	0.22	0.10	0.60
KPC	9	1.44	0.23	0.44	0.50	<0.001*	0.49
LN	14	1.54	0.42	0.54	0.23	<0.001*	0.39
MB	20	1.37	0.15	0.37	0.59	<0.001*	0.40
YR	12	1.5	0.33	0.50	0.36	<0.001*	0.40
COAST
BOC	6	1.37	0.25	0.37	0.34	0.07	0.60*
CD	6	1.41	0.29	0.41	0.31	0.41	0.53*
CR	16	1.59	0.33	0.59	0.44	<0.001*	NA
COR	6	1.4	0.22	0.40	0.47	0.03	0.65*
CX	12	1.57	0.40	0.57	0.32	<0.001*	0.44
EC	7	1.54	0.43	0.54	0.25	0.11	0.42
GB	12	1.35	0.18	0.35	0.50	<0.001*	0.51
KD	12	1.59	0.41	0.59	0.31	<0.001*	0.45
KW	24	1.63	0.42	0.63	0.34	<0.001*	0.46
LD	8	1.38	0.18	0.38	0.55	0.01*	0.50*
LT	9	1.57	0.42	0.57	0.30	0.003*	0.34
MON	6	1.34	0.21	0.34	0.42	0.02	0.53
WB	5	1.21	0.05	0.21	0.79	0.07	0.83
WC	12	1.35	0.09	0.35	0.75	<0.001*	0.51*
WL	18	1.39	0.30	0.39	0.22	<0.001*	0.45
WOC	8	1.26	0.14	0.26	0.46	0.02	0.69

N is the sample size, AR is allelic richness standardized for sample size, HO and HE are the observed and expected heterozygotes, FIS is the inbreeding coefficient (FIS > 0 evidence of inbreeding), the P value relates to test of Hardy–Weinberg equilibrium with significant values after Bonferroni correction indicated in bold (P < 0.0019) and after B-Y correction indicated with an asterisk (P < 0.013), M ratios indicating bottlenecks (Arlequin; M < 0.68) are highlighted in bold, and asterisks indicate significant excess heterozygosity and/or shifted mode as calculated in BOTTLENECK. NA indicates the test was not performed (see Methods).

AB, Abercrombie R; BR, Broken R; BUFF, Buffalo Ck; CO, Cotter Dam; DD, Dartmouth Dam; HU, Hughes Ck; KPC, King Parrot Ck; LN, Lachlan R; MB, Murrumbidgee R; YR, Yarra R; BOC, Bowen Ck; CD, Cataract Dam; CR, Cataract R; COR, Cordeaux Dam; CX, Coxs R; EC, Erskine Ck; GB, Glenbrook Ck; KD, Kedumba R; KW, Kowmung R; LD, Loddon Ck; LT, Little R; MON, Mongarlowe R; WB, Webbs Ck; WC, Wheeny Ck; WL, Wollemi Ck; WOC, Wongawilli Ck.

Population structure

There was a high degree of population structure across the species range (FST = 0.43, P < 0.001). Pairwise fixation indices (Appendix 1) had a large range (FST = 0.00–0.68) and were highest between the Webbs and Wongawilli populations on the coast. Geneland separated samples into four major groups, which correspond broadly with catchment structure: Class 1. upper-mid-HN (Kowmung, Coxs, Kedumba, Little, Erskine, Glenbrook), Class 2. lower HN (Bowens, Cordeaux, Wollemi, Wongawilli, Webbs), Class 3. Murrumbidgee and Class 4. Lachlan and Murray. The inclusion of Cordeaux Dam and Wongawilli with the lower HN group (Class 2) may indicate that some long-distance dispersals or localized translocations have occurred in the HN. However, there were no reliable estimates of migration among the coastal populations as BayesAss+ results were inconsistent between runs (<7 of 10 runs). There was very little recent gene flow across the GDR, with the proportion of migrant individuals per generation being 0.004 from the HN to the MDB (CI 0.000–0.016) and 0.010 from the MDB to the HN (CI 0.000–0.036). In contrast, estimates of recent gene flow among populations of the MDB consistently (10 of 10 runs) indicated significant amounts of migration from the Abercrombie to the Lachlan [0.231 (CI 0.121–0.314)]. However, estimates of migration among the remaining MDB populations were inconsistent between runs (<7 of 10 runs).

The model that best explained levels of genetic diversity within populations included just one environmental variable: river slope (Table 2). Other highly ranked models retained river slope but also included distance from the headwater source (ΔAICc = 0.2892) and land use (ΔAICc = 0.8075). Hierarchical partitioning of the data indicated that the majority of variation in genetic diversity among populations could be explained by variation in the river slope (58.9%) and distance from the source (21.5%). However, the only significant variable in all of the models was river slope, which was positively correlated with genetic diversity (AR) (Fig. 2). The results of our regression analyses were shown to be robust to variation in sample size. Our re-sampling procedure indicated that the coefficient for river slope was always positive [0.041 CI (0.029, 0.052)]. Mantel tests revealed that both riverine distance and the number of barriers between sites are correlated with genetic differentiation (FST) among the coastal populations (Appendix 3). The number of barriers that separate localities explains a significant portion of genetic divergence across localities (R2 = 0.25, P = 0.0245), even when riverine distance is factored out in partial Mantel correlation analysis. The same results were obtained when we recalculated FST using a random sample of five individuals from each population.

Table 2

Models of environmental variables to explain the distribution of genetic variation at eight microsatellite loci (AR) among populations of Macquaria australasica

Model rank	River slope	km from source	Barriers	Land use	Elevation	Latitude	Intercept	AICc	deltaAICc	Rsqu (adj)	ANOVA	Sig variables
1	0.0592						1.43	−29.3577	0	0.228	0.019	rs
2	0.0551	0.0357					1.43	−29.0685	0.2892	0.278	0.024	rs
3	0.0657			0.0436			1.44	−28.5502	0.8075	0.26	0.03	rs
Null (12)								−25.5947	3.7631
Full (63)	0.0679	0.0359	−0.0014	0.0343	−0.0206	−0.0456	1.46	−15.2042	14.1535	0.145	0.242	rs
HP%IC	58.9	21.5	4.8	8.5	2.6	3.7

All variables are standardized to mean of zero and standard deviation of one. HP %IC, independent contribution of each environmental variable to explaining the variation in AR; AICc, Akaike's Information Criterion adjusted for small sample sizes; R2 (adj), R2 adjusted for the number of variables in the multiple regression analysis; ANOVA, P value for the analysis of variation of environmental variables in the multiple regression (significance at P < 0.01 following Bonferroni correction); Sig., variables making a significant contribution to the model according to the ANOVA; rs, river slope.

Figure 2

Relationship between genetic diversity (AR) and river slope (m/km).

Discussion

Our study of an endangered freshwater fish shows that recent habitat fragmentation leads to negative genetic consequences. Genetic diversity in M. australasica was significantly higher in locations with a higher river slope, a surrogate for riffles, the species preferred habitat structure. In contrast, barriers both fragment and degrade preferred habitat and impede dispersal, contributing to the genetic differentiation of populations. We highlight the importance of the landscape genetics approach to understand environmental factors affecting freshwater fish populations and its potential practical application to conservation management of other freshwater organisms with specific habitat requirements and limited dispersal ability.

Preferred habitat promotes genetic diversity

Macquaria australasica is known to be a habitat specialist, as extensive habitat mapping and surveying within some of the most abundant remnant populations (the upper Lachlan and Abercrombie Rivers and the HN) suggest that the area of riffles within a reach is one of the best predictors of the species presence (Bruce et al., 2008; Gilligan et al. 2010). Species-habitat association modelling indicates that at least one ha and an optimum of three ha of riffle habitat per kilometre of stream are required for M. australasica to be present (Gilligan et al. 2010). Riffles are known to be a vital factor for successful recruitment of M. australasica populations (Cadwallader and Rogan 1977). These areas of the stream habitat are where demersal eggs are laid and remain until hatching and where larvae shelter during development. Therefore, a greater area of riffles provides a greater area of suitable spawning substratum and may also improve the survival of larval fish. Gilligan et al. (2010) indicated that the Abercrombie River had more optimal habitat for the species, with 56% of reaches containing suitable habitats, including riffles, compared with 20% in the Lachlan River. Our measure of riffle habitat, river slope, was also higher in the Abercrombie than in the Lachlan River (6.4 and 4.3 m/km, respectively), supporting our use of this measure as a surrogate for the species preferred habitat.

Although natural river structures and river slope have been observed to influence the levels of genetic structure in other species (Wofford et al. 2005; Boizard et al. 2009; Blanchet et al. 2010; Cook et al. 2010), we believe there is also a clear need to consider the patterns and processes that influence genetic diversity of populations. This is especially true in endangered species where management resources need to be directed to populations of high priority, such as those with reduced levels of genetic diversity (Amos and Balmford 2001). Maintaining the resilience and evolutionary potential of such populations could be facilitated by enhancing the environmental processes identified by landscape genetics studies to be important drivers of genetic diversity.

Anthropogenic structures contribute to population genetic structure and restricted gene flow

Overall, broad-scale population genetic structure reflected drainage basin divisions with the clear differentiation of the MDB and HN, as well as the major catchments within these drainage basins. There was no evidence of significant amounts of recent gene flow across the GDR. This provides further support to our finding of historically isolated lineages based on mtDNA data (Faulks et al. 2010a) as well as studies that showed allozyme and morphological differentiation among populations across the GDR (Dufty 1986). In addition, many characteristics of the species biology and ecology are thought to differ on either side of the GDR (A. Bruce, personal communication). Therefore, support for the designation of the coastal and MDB forms as separate evolutionary significant units (ESUs) (sensuCrandall et al., 2000) is strong, and further investigation into the taxonomic distinction of these forms should be considered a priority.

Although the movement patterns of M. australasica are not well documented, the species has been observed undertaking small migrations within a localized home range (Katie Ryan, University of Canberra, personal communication). Despite relatively limited dispersal behaviour, we have demonstrated that habitat fragmentation, in the form of anthropogenic structures, has led to increased genetic differentiation of populations. We hypothesize that M. australasica may undertake rare long-distance dispersal events that are hindered by anthropogenic structures. Alternatively, habitat fragmentation may have a cumulative negative effect on the species genetic architecture. Until further studies clarify patterns of dispersal behaviour in M. australasica, the specific effect that anthropogenic structures exert on the species will remain unclear.

The separation of populations within the MDB is likely driven by a combination of extensive geographical distances, unsuitable habitat conditions as well as the presence of major anthropogenic structures, e.g. Cotter Dam, Dartmouth Dam and Wyangala Dam. The only evidence of recent migration was from the Abercrombie River into the Lachlan River, adjacent waterways in the same catchment. Migration between the two populations appears unidirectional with the Abercrombie potentially acting as a source population for the Lachlan. It is possible that greater habitat quality and quantity within the Abercrombie River (Gilligan et al. 2010) result in a larger population size and consequent Ne within the Abercrombie River population. Within the HN basin, the two population groups (upper-mid-HN and lower HN) are separated not only by geographical distance, in some cases the saline tidal waters in the Hawkesbury River estuary, but also by anthropogenic barriers. For example, the upper-mid-HN group contains populations above Warragamba and Cordeaux Dams; these are large impoundments that likely act as effective barriers to dispersal in both downstream and upstream directions.

Management

Macquaria australasica has experienced widespread population declines (Bruce et al., 2008; Davies et al. 2008), and despite extensive sampling effort, many of the populations in this study are represented by relatively small sample sizes. We acknowledge that sample sizes may have contributed to the Hardy–Weinberg disequilibrium observed in some populations as well as reduced the power of BayesAss+ to detect recent migration. Nonetheless, simulation studies with BayesAss+ have shown that increasing sample size has only a moderate effect on the accuracy of estimates, particularly when the number of loci and migration rates are low (Wilson and Rannala 2003). Importantly, the major findings of our study – the influence of habitat fragmentation (anthropogenic and natural in-stream structures) on genetic structure and diversity – are supported by re-sampling procedures designed to account for variable and small sample sizes. Overall, we believe that our results have significant biological meaning and can contribute to the conservation management of this species.

Levels of genetic variation in M. australasica were very low, both in absolute values and in comparison with the congeneric M. ambigua (Faulks et al. 2010b). Our analyses also provided strong evidence for population bottlenecks and inbreeding in M. australasica. The tests performed are sensitive to the detection of inbreeding within the past three generations (Piry et al. 1999; Goudet 2002), a period of at least 10 years for M. australasica. The widespread construction of barriers to movement as well as declines in population abundances across the species range has occurred within this timeframe (Cadwallader and Rogan 1977) and has most likely contributed to the genetic deterioration of populations. We suggest that in order to help boost genetic variation, there is a need to mitigate the effects of habitat fragmentation and barrier construction (Hughes 2007). The installation of fishways to assist dispersal of fish over dams and weirs has proven beneficial for rehabilitating fish communities in Australia and overseas (Calles and Greenberg 2007; Stuart et al. 2008) and should be considered as an important component in the management of M. australasica. Allowing passage over barriers could provide M. australasica populations with access to additional spawning habitat and enhance recruitment success. The rehabilitation of preferred habitat, e.g. removing sand or sedimentation from riffles or enhancing environmental flows to scour sediment loads, could also help increase population sizes. Other management options being considered in the current species recovery planning process include reducing illegal fishing, reducing disease risk, establishing new populations, increasing community awareness and establishing a long-term monitoring programme (I&I NSW, 2010). In addition, we recommend that the following distinct catchment-scale groups be recognized to help maintain the genetic diversity, the evolutionary resilience and the evolutionary potential of the species as a whole: Lachlan, Murrumbidgee, Murray, upper HN and lower HN. In conclusion, we highlight the value of incorporating landscape genetics studies into established frameworks such as the ESU, particularly in cases where conservation outcomes can be achieved and assessed through adaptive management strategies.

	AB	BUFF	DD	BR	CO	HU	KPC	LN	MB	YR	CD	LD	MON	BOC	COR	CR	CX	EC	GB	KD	KW	LT	WB	WC	WL	WOC
AB	0.00
BUFF	0.14	0.00
DD	0.15	0.01	0.00
BR	0.11	0.08	0.07	0.00
CO	0.28	0.49	0.38	0.44	0.00
HU	0.16	0.12	0.08	0.12	0.35	0.00
KPC	0.09	0.16	0.18	0.13	0.37	0.13	0.00
LN	0.06	0.05	0.07	0.09	0.35	0.12	0.12	0.00
MB	0.19	0.29	0.25	0.30	0.30	0.22	0.21	0.26	0.00
YR	0.09	0.16	0.11	0.15	0.33	0.14	0.15	0.13	0.22	0.00
CD	0.12	0.18	0.13	0.15	0.34	0.16	0.09	0.12	0.22	0.15	0.00
LD	0.24	0.31	0.23	0.32	0.37	0.27	0.24	0.14	0.32	0.31	0.16	0.00
MON	0.19	0.22	0.22	0.25	0.36	0.19	0.12	0.21	0.02	0.16	0.16	0.29	0.00
BOC	0.44	0.58	0.47	0.53	0.61	0.47	0.52	0.42	0.53	0.46	0.53	0.53	0.55	0.00
COR	0.43	0.54	0.44	0.50	0.58	0.44	0.49	0.40	0.52	0.44	0.50	0.50	0.53	0.28	0.00
CR	0.16	0.14	0.10	0.17	0.28	0.15	0.15	0.09	0.18	0.14	0.09	0.13	0.17	0.35	0.33	0.00
CX	0.34	0.45	0.37	0.40	0.48	0.38	0.40	0.33	0.44	0.34	0.40	0.44	0.44	0.32	0.27	0.23	0.00
EC	0.38	0.48	0.40	0.40	0.51	0.40	0.43	0.37	0.50	0.32	0.41	0.48	0.47	0.43	0.43	0.28	0.15	0.00
GB	0.47	0.60	0.55	0.55	0.62	0.53	0.55	0.48	0.60	0.49	0.55	0.58	0.59	0.51	0.57	0.39	0.32	0.17	0.00
KD	0.36	0.39	0.35	0.34	0.46	0.36	0.37	0.32	0.43	0.34	0.38	0.41	0.41	0.29	0.25	0.26	0.05	0.17	0.29	0.00
KW	0.33	0.37	0.32	0.33	0.40	0.33	0.34	0.30	0.36	0.30	0.35	0.38	0.37	0.27	0.23	0.22	0.00	0.14	0.27	0.03	0.00
LT	0.35	0.52	0.44	0.41	0.56	0.44	0.42	0.35	0.52	0.40	0.46	0.47	0.49	0.44	0.40	0.26	0.16	0.25	0.41	0.09	0.10	0.00
WB	0.50	0.67	0.60	0.60	0.71	0.59	0.59	0.50	0.66	0.56	0.62	0.62	0.65	0.46	0.52	0.44	0.40	0.44	0.55	0.34	0.32	0.42	0.00
WC	0.44	0.53	0.48	0.51	0.61	0.50	0.51	0.42	0.55	0.48	0.53	0.51	0.54	0.54	0.54	0.37	0.46	0.51	0.58	0.42	0.38	0.52	0.64	0.00
WL	0.46	0.57	0.53	0.52	0.61	0.51	0.52	0.46	0.57	0.49	0.55	0.55	0.57	0.45	0.54	0.40	0.39	0.34	0.34	0.35	0.35	0.36	0.46	0.55	0.00
WOC	0.48	0.65	0.53	0.60	0.65	0.51	0.57	0.45	0.60	0.53	0.59	0.54	0.63	0.55	0.34	0.40	0.37	0.52	0.64	0.36	0.30	0.49	0.68	0.63	0.60	0.00

	km	Barriers	rs	Latitude	Landuse
Barriers	0.496
	0.026
rs	0.121	0.288
	0.611	0.218
Latitude	−0.129	0.052	0.483
	0.588	0.829	0.031
Land use	0.191	0.169	−0.211	−0.084
	0.419	0.476	0.372	0.725
Elevation	0.466	0.571	−0.193	−0.649	0.063
	0.038	0.009	0.415	0.002	0.793

Variable	P value	Coefficient	R2
Riverine km	0.0005	0.429	0.184
Barriers	0.0005	0.492	0.243
Barriers (km controlled)	0.0245	0.429	0.250
km (barriers controlled)	0.0005	0.257	0.250