An investigation of genetic connectivity shines a light on the relative roles of isolation by distance and oceanic currents in three diadromous fish species

Connectivity is critical to species survival, but in many instances is compromised through human impacts (Noss et al. 2012; Haddad et al. 2015; Wilson et al. 2016). Demographic connectivity benefits relate to colonisation, population reestablishment and source and sink dynamics (Hanski 1998). Connectivity also has genetic benefits, buffering local populations against loss of genetic diversity by genetic drift, reducing inbreeding depression and promoting adaptive potential (Lowe and Allendorf 2010; Smith et al. 2014). Given its importance for species survival, conservation management practices are increasingly incorporating connectivity as a key component of long-term species management and recovery plans while actively considering the relative connectivity benefits from both a genetic and demographic perspective (Almany et al. 2009; Carroll et al. 2015; Magris et al. 2016).

What influences connectivity at a population level is often the result of complex interactions between present-day geography, the geographic history of a region, species biology and species-specific density-dependent processes, making population connectivity often highly specific to a species or area (Riginos et al. 2011; Waters et al. 2013). In many systems, large physical biogeographic barriers (e.g. land bridges, rivers and mountain ranges) can block or impede population connectivity, resulting in common patterns of population structure among diverse taxa (Pascual et al. 2017; Thacker 2017; Sánchez-Montes et al. 2018). Although physical barriers are easily observed, more subtle barriers can influence connectivity, such as resource and temperature differences in habitats limiting movement between populations by reducing the fitness of migrants arriving from different conditions (Friesen 2015). Population connectivity can also be influenced by life history traits, such as fecundity, size or either sessile or motile larval stages, which influence individual dispersal potential (Selkoe and Toonen 2011; Stevens et al. 2012; Whitmee and Orme 2013). Further, barriers and biological traits often interact, creating a relationship where stronger dispersers are less affected by barriers (Pascual et al. 2017). Ultimately, however, if the barriers to connectivity are sufficiently strong, even strong dispersers will lose population connectivity (Pascual et al. 2017). In addition to the interacting nature of barriers and biological traits, as species density increases many species have been shown to exhibit higher rates of dispersal and connectivity as a means of reducing overcrowding and resource limitations, further complicating the nature of species-level connectivity (Bitume et al. 2013, 2014).

Many studies have sought to predict connectivity based on the interactions between environment and species-specific biology (Bradbury et al. 2008; Stevens et al. 2014; Bonte and Dahirel 2017). These studies have helped identify many key factors influencing connectivity within individual species and groups of species that share similar life history traits (Stevens et al. 2014; Bonte and Dahirel 2017). However, these factors often fall short of adequately describing factors that universally influence connectivity across all scenarios (Clobert et al. 2012), highlighting the importance of species-specific connectivity studies to accurately infer population connectivity.

Riverine systems pose unique challenges to population connectivity for aquatic species. The highly linear and dendritic nature of rivers limits the connectivity of habitats and restricts the range of ways migrants can move, and therefore how populations can be connected (Jansson et al. 2007). These systems are additionally prone to complete isolation by natural barriers such as basin boundaries or human-induced barriers (e.g. dams, weirs) associated with river regulation (Puebla 2009; Birnie-Gauvin et al. 2017). In response to the more isolated and dendritic nature of freshwater habitats compared with more continuous oceanic habitats, many species have developed unique life history traits to promote connectivity, such as diadromy (McDowall 1999).

Diadromous species exhibit two distinct lifecycle stages, one in fresh water and one in the marine environment, with individuals migrating between the two environments to use resources or obtain connectivity benefits for a particular life stage (McDowall 1999). Among diadromous fishes, several specific modes of migration are recognised, including anadromy (mature adults migrate upstream from the sea to spawn in fresh water), catadromy (mature adults migrate from fresh water to spawn in the sea) and amphidromy (adults spawn in fresh water, with newly hatched larvae moving to the marine environment and returning to fresh water as juveniles; Myers 1949; McDowall 1988). Diadromous species typically show connectivity at a level intermediate to both freshwater and marine species due to their ability to use the highly connected marine environment for dispersal between otherwise isolated rivers (Ward et al. 1994). The degree of connectivity can vary greatly across diadromous species and is generally considered dependent on several life-history traits. Many life history traits have been shown to influence connectivity, with higher connectivity typically seen in diadromous species that show higher fecundities, larger body size or long-lasting pelagic larval stages, or those that actively disperse (Selkoe and Toonen 2011; Clobert et al. 2012; Feutry et al. 2013; Jones and Closs 2016). Although some economically important diadromous species such as salmonids have been studied extensively, diadromous species as a whole are a critically understudied group of fish, with many key questions surrounding how diadromy affects the overall connectivity seen within different diadromous species remaining (McDowall 1999; Delgado and Ruzzante 2020). This has ultimately resulted in insufficient knowledge of life history for many species to estimate dispersal.

Recent advances in genomic techniques have allowed for large-scale sequencing of genome-wide markers, producing thousands of markers compared with typical panels of 10–15 microsatellites traditionally used to measure genetic connectivity (Schmidt et al. 2011, 2014; Cook and Sgro 2017). With a larger number of genome-wide markers, more subtle patterns of connectivity can be resolved, allowing for finer-scale insights into population structure and connectivity. In addition, the use of large genomic datasets provides highly resolving information on important metrics of species viability, such as genetic diversity, effective population size and adaptive potential, helping highlight populations of conservation value and conservation concern (Cook and Sgro 2017).

In this study we used large-scale genomic datasets with thousands of single nucleotide polymorphic genetic markers to investigate patterns of genetic connectivity for three diadromous species: common galaxias (Galaxias maculatus), tupong (Pseudaphritis urvillii) and Australian grayling (Prototroctes maraena). These three species all use marine and freshwater environments and vary in several different life history traits commonly linked to dispersal, suggesting that connectivity may differ between the species (Miles 2007; Clobert et al. 2012). However, knowledge of the life history of each species is limited, with many key traits that can influence dispersal, such as larval mobility, unknown. Given such knowledge gaps surrounding key dispersal-linked life history traits, these species are ideal candidates for genomic-based assessments of connectivity. The aims of this study were to: (1) determine the extent of population genetic connectivity, levels of genetic diversity and effective population size estimates for each species; (2) investigate the role isolation by distance (IBD) and isolation by resistance (IBR) influenced by oceanic conditions play in population genetic connectivity for each species; and (3) determine whether differences in conservation status are reflected in the patterns of connectivity and genetic diversity observed for each species.

The three study species exhibit differences in life history traits commonly linked to dispersal, including fecundity, larval size and movement strategies, and pelagic larval duration time, all of which may contribute to connectivity differences between each species (Table S1 of the Supplementary material). However, all species studied lack information on other key life history traits associated with larval dispersal, such as mobility, growth rates and preference for ocean depth. Further, many life history traits that are known for one species are not known for others, such as initial larval length and egg size (Table S1).

Common galaxias G. maculatus is a widespread species with a distribution extending across Australia, New Zealand and South America (Gomon and Bray 2011). The form of diadromy exhibited by common galaxias is still debated, but is considered either amphidromous or semicatadromous, with adults migrating into estuaries to breed and larvae drifting into the marine environment before immigrating into rivers as juveniles (McDowall 2000; Bice et al. 2019b; Delgado and Ruzzante 2020). Adults produce thousands of eggs that hatch and undertake an ~5 month marine larval phase (McDowall et al. 1994). Otolith studies on New Zealand populations suggest high dispersal capabilities and population connectivity for the species along the New Zealand east and west coasts (Hickford and Schiel 2016). Genetic studies on populations along the Chilean coast show that the species exhibits very high diversity and connectivity (Delgado et al. 2019), but shows population genetic structuring at large spatial scales (>1000 km; Zemlak et al. 2010) as a result of past glacial events along the southern coast (González-Wevar et al. 2015). Its widespread distribution has led to common galaxias being classified as of least concern under International Union for Conservation of Nature (IUCN) status (Bice et al. 2019b).

Tupong P. urvillii is a catadromous species endemic to the coastal drainages of south-east Australia (Gomon et al. 2008). Adults spawn during late autumn–winter in the sea (Crook et al. 2010), producing up to 400 000 eggs before juveniles enter rivers from the sea after approximately 4 months (Zampatti et al. 2010; Bice et al. 2012; Walker and Humphries 2013). Previous genetic work using microsatellite markers suggests tupong exhibit weak IBD structure (Schmidt et al. 2014). Its relative abundance through south and eastern Australia has led to tupong being classified as of least concern under IUCN status (Bice et al. 2019a).

Australian grayling P. maraena is an amphidromous species. Adults mature and spawn in fresh water in autumn, producing up to 47 000 eggs, with the eggs and larvae drifting downstream to the sea and juveniles migrating back into fresh water after ~4–6 months (Berra 1982; Crook et al. 2006; Koster et al. 2013, 2021). Australian grayling has a similar distribution to tupong, but smaller estimated population sizes throughout Victoria (Allen et al. 2002), which, coupled with population declines, has led to Australian grayling being classified as a vulnerable species under the Environment Protection and Biodiversity Conservation Act 1999 of Australia and as vulnerable under IUCN status (Koster and Gilligan 2019). Previous genetic work analysing microsatellite and mitochondrial markers suggests that Australian grayling are connected between river populations around central and eastern Victoria (Schmidt et al. 2011).

Samples of each species were taken from rivers along the Victorian coast, south-east Australia (Fig. 1). The region contains two major oceanic currents, the East Australian Current (EAC), which flows south along the mainland eastern Australian coast towards Tasmania, and the South Australian Current (SAC), which flows eastward from the southern coast of Australia (Tilburg et al. 2001; Ridgway and Condie 2004). These currents vary seasonally, with the SAC flowing further into Bass Strait in winter and the EAC flowing further into Bass Strait during summer, creating a mosaic of directional flow depending on the season (Baird et al. 2006; Ridgway 2007).

Fig. 1.  Map of sampling locations for common galaxias, tupong and Australian grayling across the coastline of Victoria, Australia. Stars show the collection site in each river. Numbers indicate which species were collected along each river: 1, common galaxias; 2, tupong; 3, Australian grayling. The two-letter codes for each river are defined in Table 1. The South Australian Current (SAC) and East Australian Current (EAC) are shown, with arrows indicating the primary water flow direction.

Samples were collected from 94 common galaxias from 10 Victorian rivers (Table 1). Of these 94 individuals sampled, 69 were juveniles returning after their marine larval phase during Spring 2017 and 24 were from the same cohort sampled as young-of-the-year (YOY) during Autumn 2018 (Table 1). In addition, 94 tupong individuals were sampled from 7 rivers in Victoria. Of these 94 individuals, 28 were juveniles returning after their marine larval phase during the spring–summer of 2016 and 66 were from the same cohort sampled as YOY in February 2017 (Table 1). A total of 162 Australian grayling juveniles was sampled from 9 Victorian rivers and across three different year cohorts. Of these 162 individuals, 101 were juveniles collected after returning from their marine larval stage during the spring of 2016 and 27 were sampled from the same cohort as YOY collected in February 2017, 12 were collected as returning juveniles of a second-year cohort during the spring of 2018 and a further 22 were sampled as returning juveniles from a third-year cohort in the spring of 2019 (Table 1). All fin clips collected for common galaxias and tupong were part of the same cohort, but Australian grayling individuals were composed of two cohorts. All fin clips collected were stored in 95–100% ethanol and refrigerated at 4°C until later DNA extraction.

Table 1.  Sampling numbers taken from each study river for each species
The coordinates of each sampling site are provided in Table S2

Sampling of individuals for this study was undertaken under the ethics permits AEC15/005 and AEC18/003 (ARI Animal Ethics Committee).

For common galaxias, ~5 mm² of tissue was cut from the caudal or dorsal fins of each individual, and the tissue clips sent to Diversity Arrays Technology (DArT; Canberra, ACT, Australia), where total genomic DNA was extracted. For tupong and Australian grayling, DNA was extracted from small (~5 mm²) fin clips using a DNeasy Blood & Tissue kit (Qiagen) according to the manufacturer’s instructions. Approximately 20 µL of extracted DNA (at concentrations ranging from 5 to 20 ng µL^–1) for each sample was sent to DArT for sequencing. The DNA was then sequenced using the Dartseq platform (Diversity Arrays Technology; www.diversityarrays.com). This platform uses a form of reduced representation sequencing similar to double-digest restriction site-associated DNA sequencing to generate large numbers of single nucleotide polymorphisms (SNPs). DNA samples were digested using two different enzymes, namely PstI and SphI, with two adaptors corresponding to the two different restriction enzyme overhangs, following the digestion and ligation methods described by Kilian et al. (2012). The forward (PstI) adaptor included an Illumina flow cell attachment sequence, sequencing primer sequence and a unique barcode for multiplexing, whereas the reverse adaptor included an Illumina flow cell attachment sequence. Only DNA fragments containing both cut sites were then amplified for sequencing because fragments that lacked one or both cut sites were unable to amplify. Equimolar amounts of each amplified product were then pooled and sequenced as single reads on an Illumina Hiseq 2500 for 77 cycles. Samples were sequenced in batches of 94 per Illumina sequencing lane, with 25% of samples rerun as technical replicates for quality control.

Sequenced reads were processed using DArT proprietary analytical pipelines described in Kilian et al. (2012), with poor-quality sequences removed (based on the certainty of base calls and concordance with internal replicates run on a random subsample of 25% of individuals sequenced) and low-quality bases corrected using collapsed tags from multiple members as a template. A secondary pipeline (DArTsoft14) was then used to compile read counts into SNP loci calls with a maximum allowed difference between tags of three bases. SNPs were further filtered to remove loci whose allele read counts had a greater than fivefold difference from each other, and that scored <95% reproducibility using the sequenced technical replicates. Following the analytics pipelines, 58 997, 16 555 and 27 783 SNPs were genotyped for common galaxias, tupong and Australian grayling respectively.

SNPs processed using the DArT analytical pipelines were further filtered using the dartR package (ver. 1.1.11, CRAN.r, see https://CRAN.R-project.org/package=dartR; Gruber et al. 2018) in R (ver. 3.6.0, R Foundation for Statistical Computing, Vienna, Austria). All SNPs were filtered at a reproducibility value of 1, only retaining SNPs that were consistent in 100% of technical replicates sequenced as part of the Dartseq process (Table S3 of the Supplementary material). The data were filtered to remove any SNP that was missing in >10% of individuals, and only a single SNP per sequence tag was retained. All individuals were then filtered for a maximum allowable amount of missing data of 15%. All SNPs were then filtered to remove alleles that were present in <2% of all individuals from each species. After filtering, 8116 SNPs were retained for common galaxias, 6762 SNPs were retained for tupong, and 3410 SNPs were retained for Australian grayling.

After filtering, all loci were investigated for Hardy–Weinberg equilibrium (HWE) using the dartR package in R. All markers were found to be in HWE for each species, and therefore no markers were removed at this stage.

Two methods were used to test for neutrality within the genetic markers used. The first method used was the fixation index (F_ST) outflank method implemented in the dartR package in R, whereas the second method used was through the program Bayescan (ver. 2.1, University of Bern, see http://www.cmpg.unibe.ch/software/BayeScan/; Foll and Gaggiotti 2008). For the F_ST outflank method, left and right trims of 5% of loci were used, and an α value of 0.05 Bonferroni corrected was used to identify putatively selected loci. For Bayescan, the Markov chain was run for 50 000 steps with a burn-in of 50 000 steps, a thinning interval of 10 and a total of 40 pilot runs composed of 10 000 steps. When plotting the results of Bayescan, a false discovery rate of 0.05 was used. Both methods identified zero loci for tupong and Australian grayling, and so no loci were filtered for either species before further analysis (Fig. S1–S5 in the ‘Neutrality tests of genetic markers’ section of the Supplementary material). Two loci were identified as putatively non-neutral for common galaxias, but were found along the border of significance (see the ‘Neutrality tests of genetic markers’ section of the Supplementary material). To test the impact of each locus, a subset of analyses was repeated without the two loci present. Because all analyses returned near identical results, the two loci were not filtered out (Tables S4–S7 and Fig. S1–S5 of the ‘Neutrality tests of genetic markers’ section of the Supplementary material).

Population genetic summary statistics including observed heterozygosity (H_o) expected heterozygosity (H_e) and allelic richness (AR) were calculated for each population and as species averages using the R software package hierfstat (ver. 0.04–22, CRAN.r, see https://cran.r-project.org/web/packages/hierfstat/index.html; Goudet 2005), and percentage polymorphic loci were calculated using custom code in R.

Kinship analysis was used to identify sibling relationships between individuals, and therefore reveal connectivity patterns by identifying whether siblings were found to return to the same rivers or different rivers to each other, indicative of individual dispersal events. Further, identifying and removing sibling pairs can allow for increased accuracy when estimating population structure, dispersal and effective population sizes, because kinship can often mask true patterns of connectivity within datasets if unaccounted for (Peterman et al. 2016). Kinship analyses for each species were performed using the program Colony (ver. 2.0.6.5, Zoological Society of London, see https://www.zsl.org/science/software/colony; Jones and Wang 2010) to test for full-sibling and half-sibling pairs found across each species. Colony was run under the settings of combined pairwise and full-likelihood score (setting 2 out of 0–2), a genotyping error rate of 1%, updating allele frequencies (setting 1 out of 0/1), assumed polygamy for both males and females (setting 0 0 out of 0–1, 0–1), inbreeding, a long run (setting 3 out of 1–4) and for four runs that converged to average the most likely scenario. Using the results of these analyses, one individual from all full-sibling pairs was removed from each dataset for the analyses involving population structure, IBD, migration and oceanic resistance modelling.

To estimate the degree of population genetic structure among river populations within each species across the same geographic area, three different methods were used: (1) genetic differentiation (F_ST); (2) discriminate analysis of principal components (DAPC); and (3) and non-spatial clustering (Structure).

Pairwise F_ST values, which measure the degree of genetic differentiation between pairs of populations based on allele frequencies, were calculated using the R software package StAMPP (ver. 1.5.1, CRAN.r, https://cran.r-project.org/web/packages/StAMPP/index.html; Winter 2012). The significance of values was calculated using 100 000 bootstrap replicates of the data.

DAPC analyses were undertaken to estimate the number of genetic populations that exist within each species while also assigning individuals into their most likely genetic population. For DAPC, all individuals were assigned into optimal groups using successive K-means clustering from the find.clusters function in the R package adagenet (ver. 1.3.1, CRAN.r, see https://cran.r-project.org/web/packages/adegenet/index.html). For each species, K-means clustering suggested that in decreasing likelihood K = 1–4 were optimal based on the Bayesian inference criterion, so Clusters 1–4 were created.

Structure analysis was undertaken to complement DAPC, estimating the degree of spatial population genetic structure present in each species (Pritchard et al. 2000). For each species, Structure analysis was run with the number of genetic clusters (K) set from 1 to 10, with 14–17 (run time total of 50 CPU days per structure run) replicate runs of 200 000 Markov chain Monte Carlo iterations, after an initial burn-in period of 100 000 iterations. The results of each run were extracted using the program Structure Harvester (ver. 0.6.94, University of California, Los Angeles, CA, USA, see http://taylor0.biology.ucla.edu/structureHarvester/; Earl 2012), and the best fit for K was determined using both the Evanno method of a measure of the highest value of delta K (Evanno et al. 2005) and the highest overall likelihood of a structure output of fitting the data. Plots outside the optimal plots determined based on best fit K were also included to investigate any further structuring that may be present. Structure plots were created by averaging the result of each structure run using the program CLUMPP (ver. 1.1.2, Stanford University, see https://rosenberglab.stanford.edu/clumpp.html; Jakobsson and Rosenberg 2007) and visualised using the R package StructuRly (ver. 0.1.0, GitHub, see https://github.com/nicocriscuolo/StructuRly; Criscuolo and Angelini 2020).

To determine whether population structuring was the result of a pattern of IBD between populations, both a Mantel test and distance-based redundancy analysis (dbRDA) were undertaken comparing pairwise F_ST values with spatial distance. Estimates of IBD using Mantel tests were determined using latitude and longitude coordinates for each river converted to a log Euclidean distance and population genetic distance matrix using pairwise F_ST as the measure of genetic distance. Dissimilarity matrices were compared as a Mantel test with up to 999 permutations using the R software package dartR. Estimates of IBD using dbRDA were determined using standardised geographic distances and pairwise F_ST plotted using the dbRDA function from the R package vegan (ver. 2.5–6, CRAN.r, see https://cran.r-project.org/web/packages/vegan/index.html). The significance of the dbRDA relationships was then determined through an anova.cca in the vegan package.

To measure the potential presence of asymmetrical gene flow and migration between populations, relative population migration rates were determined using the divMigrate function and R shiny application (ver. 1.0 GitHub, see https://github.com/kkeenan02/divMigrate-online; Sundqvist et al. 2016). Estimates of asymmetrical migration between populations were calculated using Jost D statistics and confidence intervals were determined from 1000 bootstrap replicates. Relative migration rates are known to be more accurate when populations are defined based on the results of spatial structuring programs such as Structure (Sundqvist et al. 2016). As such, two analyses of divMigrate were undertaken: (1) an analysis where each river was listed as a distinct population; and (2) when population structuring was found, river populations were combined with neighbouring rivers into larger metapopulations based on the results of spatial structuring analysis (Structure program).

To estimate the role of geography and oceanic conditions on the patterns of genetic connectivity seen within each species, a series of oceanic resistance models was developed to investigate whether currents or dispersal from the coastline are influencing population genetic connectivity or whether the connectivity within each species reflects an IBD effect (Table 2). A total of six resistance models was used, the first three reflecting IBD with variation in how far larvae may travel from the coastline and the second three reflecting IBR based on current strength and direction with variation in how far larvae may travel from the coastline (Table 2).

Table 2.  Different oceanic resistance models developed to investigate how oceanographic conditions may affect population genetic connectivity
For more details of each resistance model and descriptions of how each resistance raster was produced, see File S1. For each completed resistance raster under each model, see File S2

Current strength and directional information were obtained from the Integrated Marine Observing System (IMOS; www.imos.org.au, accessed 3 September 2020). These data were then used to determine mean current strength and direction through Bass Strait between 2000 and 2016. The averages across the 16 years were then calculated for the specific periods each species was estimated to be undertaking a marine larval phase (common galaxias: April–September; tupong: June–November; Australian grayling: June–December; McDowall et al. 1994; Zampatti et al. 2010; Bice et al. 2018; Koster et al. 2018, 2021).

For all resistance models investigated, the R package rnaturalearthhires (ver. 0.2.0, GitHub, see https://github.com/ropensci/rnaturalearthhires) was used to create an accurate map of Bass Strait. From this map, the R packages sf (ver. 0.9–5, CRAN.r, see https://cran.r-project.org/web/packages/sf/index.html; Pebesma 2018) and raster (ver. 3.3–13, CRAN.r, see http://cran.stat.unipd.it/web/packages/raster/) were used to generate resistance rasters based on the resistance model being tested. After the creation of each resistance raster, ecological distance through the raster was generated using the R package gdistance (ver. 1.3–6, CRAN.r, see https://cran.r-project.org/web/packages/gdistance/index.html; van Etten 2017) using a von Neumann neighbourhood calculation of cost distance (directions = 4). Ecological distance was then plotted against pairwise F_ST for all models and against the estimated asymmetrical migration rates calculated through divMigrate for Models 4–6. Linear models of the relationship between genetic and ecological distance were then made and the correlation between variables, R², P-values and model fit through Akaike information criterion (AIC) were calculated per model.

For all analyses involving IBD and oceanic resistance modelling, the latitude and longitude coordinates of each river population were moved to the river mouth of each river. This was done to reduce the impact of Mantel tests being biased for rivers that were sampled further upstream and to allow for all populations to be sampled from an oceanic point in the model resistance rasters generated. All population structuring analyses and genetic diversity metrics were calculated on populations where a minimum of two individuals was successfully sequenced because two individuals have been shown to allow for sufficiently reliable estimates when the markers used were produced using highly resolving sequencing methods (Fumagalli 2013).

Effective population sizes are correlated with genetic diversity (Frankham 2007), which, in turn, is strongly correlated with adaptive potential within species (Frankham 2005). As such, we estimated effective population sizes to be used as a proxy for each species’ adaptive potential. However, effective population sizes are known to be downwardly biased by factors including migration and substructuring within populations (Waples and England 2011; Ryman et al. 2014). To address this, effective population sizes were calculated with individuals pooled into populations based on the results of population structuring from the program Structure. Effective population sizes were calculated using the linkage disequilibrium method implemented in the program NeEstimator (ver. 2.1, Molecular Fisheries Laboratory, see http://www.molecularfisherieslaboratory.com.au/neestimator-software/; Do et al. 2014) and calculated with the exclusion of all alleles rarer than the minor allele frequency thresholds of PCrit = 0.05 and PCrit = 0.02, where PCrit is the minimum frequency for alleles to be included in the analysis. In addition, 95% confidence intervals were generated for all effective population size estimates through a parametric (χ²-based method) within NeEstimator. When full-siblings were found, effective population sizes were additionally calculated with one sibling per pair removed, because the existence of any full siblings can disproportionately affect effective population size estimates for linkage disequilibrium measures (Peterman et al. 2016).

Mean heterozygosity was 0.138 (range 0.13–0.14) in Australian grayling, 0.257 (range 0.250–0.260) in tupong and 0.164 (range 0.160–0.165) in common galaxias. The mean percentage of polymorphic loci across all rivers was 80.4, 77.3 and 69.5% for Australian grayling, tupong and common galaxias respectively (Table 3). AR, which is standardised for sample size differences, was 1.388 in common galaxias, 1.289 in tupong and 1.237 in Australian grayling (Table 3).

Table 3.  Sample numbers and summary statistics for each population and each species
Rivers are listed from west to east across Victoria in descending order. AR, allelic richness; H_e, expected heterozygosity; H_o, observed heterozygosity; n, sample size; PIC, percentage polymorphic loci

One sibling pair was found in the common galaxias data; both siblings were sampled from the Tarwin River 1 month apart. Two sibling pairs were found in the tupong data, one pair caught from Cardinia Creek and the other pair caught from the Thomson River. Two sibling pairs were found for Australian grayling, one pair caught from Bunyip River and one pair caught from Curdies River. Because all sibling pairs were found within the same rivers, no dispersal events could be determined from kinship.

Estimates of pairwise F_ST varied greatly between species, with common galaxias showing the lowest mean F_ST (F_ST pairs ranging from 0 to 0.0037; 6.7% of comparisons were significant, with significant comparisons all between the Gellibrand River and other rivers, regardless of spatial distance). In comparison, Australian grayling exhibited a higher degree of spatial structuring (F_ST pairs ranging from 0 to 0.033; 52% of comparisons were significant, with significant values seen comparing most rivers except for spatially close rivers such as Cardinia–Bunyip, Thomson–Macalister and Glenelg–Curdies). Tupong exhibited the greatest degree of spatial structuring (F_ST pairs ranging from 0 to 0.049; 81% of comparisons were significant, with significant values seen between all rivers that were not immediately adjacent to each other and higher F_ST values as spatial distance increased; Tables S8–S10 of the Supplementary material).

DAPC analysis revealed differences in patterns of population structure among the three species. Common galaxias showed very little genetic structure across the study area, with individuals from each river population being assigned to the same genetic populations identified by DAPC regardless of the spatial distance of each river (Fig. S6–S9 of the ‘DAPC results’ section of the Supplementary material). By contrast, tupong showed clear structuring between the easternmost and westernmost rivers, with populations from the easternmost rivers classed as a separate genetic cluster to the westernmost rivers (Fig. S6–S9). Australian grayling showed a split between the westernmost rivers (Glenelg, GL; Gellibrand, GE; and Curdies, CU) and all other rivers when K = 2, whereas K = 3 suggested a slight split between the westernmost rivers (Thomson, TH; Macalister, MA) and the central rivers along the Victorian coast (Yarra, YA; Cardinia, CA; Bunyip, BU; Tarwin, TA; see the ‘DAPC results’ section of the Supplementary material).

Structure analysis results were largely consistent with those of DAPC analysis. Common galaxias exhibited little evidence of population structure across the study area, whereas tupong exhibited clear patterns of genetic structuring (Fig. 2, 3). For Australian grayling, there was evidence of moderate structuring between the most western rivers and central rivers, and some evidence of weak population structuring between the central and most eastern rivers (Fig. 4).

Fig. 2.  Common galaxias Structure and plots for K = 2–7 clusters. Rivers are arranged from west to east. The optimal K value was determined to be K = 2 using the highest likelihood method and K = 7 using the delta K method. The two non-spatial genetic clusters in the K = 2 structure plot are explored in the ‘DAPC results’ section of the Supplementary material. The two-letter codes for each river are defined in Table 1.

Fig. 3.  Tupong Structure plots for K = 2–4 clusters. Rivers are arranged from west to east. The optimal K value was determined to be K = 2 using both the highest likelihood and delta K methods. The two-letter codes for each river are defined in Table 1.

Fig. 4.  Australian grayling Structure plots for K = 2–4 clusters. Rivers are arranged from west to east. The optimal K value was determined to be K = 3 using both the highest likelihood and delta K methods. The two-letter codes for each river are defined in Table 1.

For common galaxias, the Structure outputs suggested that the most likely number of clusters was K = 2 using the highest likelihood method and K = 7 using the delta K method (Table S11 of the ‘DAPC results’ section of the Supplementary material). Based on the plot for the highest Structure likelihood (K = 2), most individuals were strongly assigned to either the yellow or blue genetic cluster (Fig. 2). However, the two clusters did not correspond to spatial location (Fig. 2). There was little genetic divergence between the two clusters (F_ST = of 0.0029; Table S14 of the ‘Additional analysis investigating population structuring for common galaxias K = 2’ section of the Supplementary material). Possible explanations for the two non-spatial genetic clusters are explored in the ‘Additional analysis investigating population structuring for common galaxias K = 2’ section of the Supplementary material, including potential cryptic subspecies, sex-based differences and potential assortative mating. Although no conclusion about the pattern can be confirmed, the genetic difference between the two clusters is very small and not the result of structuring across spatial distances. Given the two clusters detected by Structure were not detected by DAPC (Fig. 2; Fig. S7–S9), they may also be spurious. Based on the plot for the highest delta K (K = 7), no clear pattern of spatial genetic differentiation could be observed, with every river having approximately similar quantities of each cluster.

For tupong, the most likely number of clusters detected by Structure was K = 2 using both the delta K and highest likelihood methods (Table S12 of the ‘DAPC results’ section of the Supplementary material). For K = 2, Structure assigned all individuals in the westernmost rivers to one cluster and all individuals in the easternmost rivers to a second cluster, with admixture between the two clusters evident in the central rivers of Victoria (Fig. 3). Higher values of K for the Structure analyses did not reveal additional geographic patterns (Fig. 3).

For Australian grayling, the most likely number of clusters detected by Structure was K = 3 using both the delta K method and the highest average likelihood method (Table S13 of the ‘DAPC results’ section of the Supplementary material). For K = 3, individuals from the westernmost rivers (GL, GE, CU) were assigned to a mix of the yellow and blue genetic clusters, whereas <10% of individuals from any of the other rivers were assigned to the yellow cluster, suggesting moderate genetic structuring between populations from western Victoria and central–eastern Victoria (Fig. 4). More subtle genetic differentiation can be seen between the rivers found around central Victoria (YA, CA, BU, TA) and the rivers found around eastern Victoria (TH, MA). The central rivers were assigned predominately to the light blue cluster, whereas the easternmost rivers were an approximately equal mix of the light blue and dark blue clusters (Fig. 4).

A significant pattern of IBD was observed in tupong (r² = 0.3279, P = 0.003) and Australian grayling (r² = 0.2384, P = 0.025) using Mantel tests, whereas no significant pattern was observed in common galaxias (r² = 0.001375, P = 0.544; Fig. 5). The results obtained through dbRDA were consistent with the results of the Mantel tests, showing no significant relationship for common galaxias (adjusted r² = –0.0596, P = 0.79) and a clear relationship between pairwise F_ST and geographic distance in tupong (adjusted r² = 0.5028, P = 0.007) and Australian grayling (adjusted r² = 0.3265, P = 0.026; Table 4).

Fig. 5.  Isolation by distance plots for each coastal species showing (a) common galaxias, (b) tupong and (c) Australian grayling. Genetic distance is measured as fixation index (F_ST), whereas geographic distance is a log-transformed Euclidean distance. The significance of each species was calculated using 999 permutations. A higher density of pairwise comparison likelihoods (from permutations) are shown as colours moving from blue (low density) to orange and red (high density).

Table 4.  Summary results of dbRDA for each species using a pairwise F_ST matrix and standardised geographic distances
The r values show the proportion of variation in pairwise F_ST that is explained by the relationship between pairwise F_ST and geographic distance. F values are generated as the ratio of constrained and unconstrained total inertia of the dbRDA model divided by the degrees of freedom of the model

Across every species, asymmetrical migration rates were non-significant between river populations (Tables S16–S20 in the divMigrate results section of the Supplementary material). Tupong and Australian grayling were additionally analysed for migration between the genetically distinct regions found through Structure analysis (structuring found as the western, central, and eastern regions of the Victorian coast). After restructuring by these areas, divMigrate was run for analysis of asymmetrical migration between the regions, but all differing migration rates between each region for both species were non-significant (Tables S19, S20). No spatial population structuring was found for common galaxias using Structure and so no regional migration analysis was undertaken.

For each species, the optimal resistance model to describe population genetic connectivity varied (Table 5). For common galaxias, the optimal resistance model was Model 3, describing larvae that are not affected by current speed or direction and with a preference to water further from the coast. However, this model was non-significant and showed little correlation between raster resistance and pairwise F_ST (P = 0.3595, r² = 0.0196). For common galaxias, models using asymmetrical migration rates performed worse than models using pairwise F_ST, and models not including current information performed slightly better than models including these data, except for Model 1, which performed worse than Model 5 (Table 5). All models developed for common galaxias were non-significant (Table 5).

Table 5.  Summary results of each resistance model developed for each species
The optimal model based on lowest the AIC is in bold. More information on model parameters and saved resistance rasters is provided in Table 2 and Files S1 and S2 of the Supplementary material

For tupong, the optimal resistance model was Model 1, describing larvae that are not affected by current speed or direction and with no preference to water further or closer to the coastline (Table 5). This model showed a very high and significant correlation between raster resistance and pairwise F_ST (r² = 0.5722, P < 0.00007245). For tupong, models using asymmetrical migration rates performed worse than models using pairwise F_ST, and models not including current information performed slightly better than models using these data, except for Model 3, which performed worse than Models 4–6 (Table 5). All F_ST-based models for tupong were significant and each yielded a greater correlation than was observed using just Mantel statistics or dbRDA (Tables 4, 5; Fig. 5).

For Australian grayling, the optimal resistance model was Model 2, describing larvae that are not affected by current speed or direction but show a preference for waters closer to the coastline (Table 5). This model showed a very high and significant correlation between raster resistance and pairwise F_ST (r² = 0.455, P < 0.0007977). For Australian grayling models using asymmetrical migration, rates performed worse than models using pairwise F_ST, but models that included current speed and direction performed a mix of better and worse than the non-current models (in order of best model: Models 2, 4, 5, 1, 6, 3). All F_ST-based models for Australian grayling were significant and all except Models 3 and 6 yielded a greater correlation than was observed just using Mantel statistics or dbRDA (Tables 4,5, Fig. 5).

Three distinct genetic regions were identified within tupong and Australian grayling. Using these genetic regions, individuals from each species were pooled into their respective spatial cluster. The values generated for the eastern and western regions were incalculable (returned estimates which generated confidence ranges overlapping with infinity) for all species, likely due to low population sizes (Table 6). Values were generated for each species for the central Victorian region identified (through Structure), and so were used as the primary comparison between the three species. From this region, the effective population size of tupong across both allele frequency tests (full-siblings included: for PCrit = 0.05, N_e = 112.6 (CI 111.5–113.7) and for PCrit = 0.02, N_e = 118.4 (CI 117.3–119.6); siblings removed: for PCrit = 0.05, N_e = 1825.5 (CI 1583.6–2154.1) and for PCrit = 0.02, N_e = 1810.4 (CI 1583.3–2111.3) was significantly lower than for common galaxias (full-siblings included: for PCrit = 0.05, Ne = 692.6 (CI 675.6–710.6) and for PCrit = 0.02, N_e = 670.3 (CI 659.2–681.7); siblings removed: for PCrit = 0.05, N_e = 11 447.8 (CI 8102.1–19 483.7) and for PCrit = 0.02, N_e = 17  362.0 (CI 12 166.9–30 278.3), and varied between slightly higher and lower than Australian grayling (full-siblings included: for PCrit = 0.05, N_e = 1107.6 (CI 1035.4–1190.4) and for PCrit = 0.02, N_e = 1171.5 (CI 1123.7–1223.5); siblings removed: for PCrit = 0.05, N_e = 1124.3 (CI 1136.1–1326.8) and for PCrit = 0.02, N_e = 1295.8 (CI 1237.1–1360.2); Table 6).

Table 6.  Effective population sizes of each species with samples pooled into metapopulations based on the results of population genetic structuring
Three metapopulations were identified based on population genetic structuring: the Western region (GL, CU, GE), the Central region (BA, WE, YA, CA, BU, TA) and the Eastern region (TH, MA, MI). Effective population sizes were calculated with the exclusion of alleles rarer than two set minor allele frequencies (PCrit = 0.05 and 0.02) and both with all sibling pairs and with the removal of one individual from each sibling pair. Effective population sizes were calculated through the linkage disequilibrium method. Inf, values that extended into infinite estimates (this indicates that there is no clear evidence of any genetic drift driving genetic characterisation and can occur due to low sample sizes and large numbers of less informative loci; Do et al. 2014)

In this study we investigated the genetic connectivity and genetic diversity of three diadromous Australian fish species. Common galaxias exhibited no clear pattern of spatial genetic structuring and high genetic connectivity across the Victorian coastline (~700 km), whereas tupong exhibited moderate population genetic structuring and evidence of an IBD effect. Australian grayling showed some population genetic structuring and evidence of IBD, but to a lesser degree than tupong. In all three species, there was little evidence that oceanic conditions, such as distance from the coast, current strength and current direction, affected genetic connectivity. In addition, there was no evidence of oceanic conditions driving different relative migration rates between populations for any of the species studied. The lack of consistency in connectivity patterns highlights the importance of species-specific studies of population connectivity to inform the management of target species.

The results of Structure, F_ST analysis and DAPC suggested a clear pattern of population structuring for tupong and Australian grayling, and a pattern of high connectivity for common galaxias, across a large spatial scale (~700 km). These patterns are consistent with what is known for tupong and common galaxias from previous genetic studies (Zemlak et al. 2010; Schmidt et al. 2014; Hickford and Schiel 2016). Common galaxias is considered highly dispersive with strong connectivity across the scale of several hundred kilometres in New Zealand and Chile (Zemlak et al. 2010; Hickford and Schiel 2016). Although some population structuring has been observed for different diadromous common galaxias populations, it tends to be at the scale of >1000 km, or driven by glacial events (e.g. along the southern coast of Chile; Zemlak et al. 2010; González-Wevar et al. 2015). Because the scale of the present study was approximately 1000 km for common galaxias, the observation of high connectivity is not unexpected. The pattern of IBD in tupong also supports previous work measuring the genetic connectivity of the species that found a consistent moderate pattern of IBD across the Victorian coast, but provides greater resolution than reported in the previous work, which used microsatellites (Schmidt et al. 2014). Previous genetic analysis of Australian grayling revealed low genetic structuring, but this was undertaken on a smaller spatial scale and used less informative genetic (microsatellite and mitochondrial specifically) markers (Schmidt et al. 2011). Because the present study found evidence of slight population structuring between populations from central and eastern Victoria, the greater resolution of markers used is likely detecting more subtle signatures of population differentiation. The reasons for the differing connectivity between species could be the result of several factors, ranging from environmental conditions present during larval stages through to biological factors such as certain life history traits.

The Victorian coast is the point at which two strong and opposing currents, namely the EAC and SAC, end (Tilburg et al. 2001; Ridgway and Condie 2004). These currents have been shown to influence the population connectivity of species in different ways, including driving southerly asymmetrical gene flow in the estuary perch (EAC) or easterly gene flow in intertidal molluscs (SAC; Shaddick et al. 2011; Miller et al. 2013). In the present study, tupong and Australian grayling showed population structuring across Victoria. However, both species showed a lack of asymmetrical migration between populations or between regional clusters. This, coupled with resistance models, which showed lower correlation when oceanic conditions were included, suggests that contemporary oceanic conditions are not strongly influencing genetic connectivity within these species. Instead, given the clear pattern of IBD observed in all analyses performed, the pattern of spatial structuring seen within tupong and Australian grayling is likely predominantly the result of IBD.

Although we did not detect any evidence for an effect of ocean currents on patterns of genetic connectivity, our study sites are essentially at the end of the prevailing currents, and it has been demonstrated that the EAC is relatively weaker feeding into Bass Strait than when moving south along the Australian east coast (Ridgway 2007). If individuals were sampled from this wider range, a stronger influence of these currents may have been observed, as can be seen in other studies of species that extend beyond the ranges sampled here (Booth et al. 2007; Shaddick et al. 2011). In addition, although ocean currents may not be driving the patterns of connectivity observed, this does not mean that oceanic conditions play no role. Currents are known to be highly variable, fluctuating on scales as low as metres and varying over the scale of hours (Hays 2017). The lack of effect of oceanic conditions within the studied species may reflect averages failing to represent the true complexity of currents that may be influencing connectivity. Although these fine-scale oceanic patterns can be better incorporated into seascape genetics and connectivity, they require more robust knowledge of target species’ life history traits and the use of simulated dispersal events, which was not achievable with the study species investigated here (Paris et al. 2013).

Although many life history traits are poorly understood for the study species, two traits that shape dispersal and overall population connectivity within a species (fecundity and pelagic larval duration, PLD) have been studied for each species (Miles 2007; Clobert et al. 2012). For example, high fecundity has been shown to result in higher connectivity in fish species and marine dispersing organisms in general (Bradbury and Bentzen 2007; Woodings et al. 2018). Fecundity varied by a factor of 10–100 among our study species (Berra 1982; McDowall et al. 1994; Zampatti et al. 2010; Bice et al. 2012; Walker and Humphries 2013). If species connectivity patterns reflect fecundity differences, we would expect higher fecund species to be the most genetically connected. By contrast, fecundity was highest in the least genetically connected species (tupong) and lowest in the most genetically connected species (common galaxias). This suggests that the connectivity differences among species were unlikely to be related to fecundity.

Another key life history trait for diadromous species is PLD, with longer PLD being associated with greater dispersal potential (Selkoe and Toonen 2011). This trait aligns most closely with the observed patterns of connectivity, with common galaxias exhibiting the equal longest PLD along with Australian grayling (4–6 months), and tupong exhibiting the shortest (~120 days; Crook et al. 2006; Jung et al. 2009; Bice et al. 2018). Therefore, PLD, as an absolute number of days, may explain the higher connectivity seen in common galaxias and lower connectivity in tupong but does not explain the population structuring seen in Australian grayling. Further, the lengths of the PLD observed for each species are all considered very long compared with many other highly connected diadromous species (Feutry et al. 2012, 2013; Taillebois et al. 2012; Teichert et al. 2012). When comparing the PLD of each species to reviews across hundreds of studied fish (both diadromous and marine), the PLD of each species studied here suggests that these species should have the capacity to disperse up to thousands of kilometres when not limited by physical geographic factors (Shanks et al. 2003; Selkoe and Toonen 2011; Faurby and Barber 2012). When looking at examples of other small-bodied diadromous species, specifically Kuhlia sauvagii, Kuhlia rupestris (rock flagtail), Sicyopus zosterophorum (belted rockclimbing goby), Smilosicyopus chloe (bilitvat), Akihito vanuatu (Vanuatu’s emperor) and Cotylopus acutipinnis (bichique), they all have PLD between two- and fivefold shorter than any species studied here while maintaining genetic connectivity across ranges of similar or larger size than the coast of Victoria (Feutry et al. 2012, 2013; Taillebois et al. 2012; Teichert et al. 2012). Although PLD is considered a crucial component of effective connectivity for diadromous species, many other life history traits (e.g. larval size, movement strategies and initial larval length) play important roles in influencing population connectivity, making it difficult to reach any conclusions regarding the effects of life history on any connectivity patterns observed in this study.

Species differences in the extent of population structuring may reflect differences in effective population size. High effective population sizes maintain genetic homogeneity by mitigating the effects of genetic drift at individual sites, increasing the time taken before population-wide genetic differences occur following isolation (Díaz-Jaimes et al. 2010). The effective population sizes calculated for each species varied depending on the test but, after the removal of siblings to mitigate the downward bias of estimates (Peterman et al. 2016), common galaxias consistently showed effective population sizes an order of magnitude higher than the other two species. The effective population sizes align with the observed patterns of connectivity, with common galaxias exhibiting the largest effective population size, followed by Australian grayling and tupong, which exhibited similar values. It is possible that dispersal of common galaxias is restricted across the study area, but that large population sizes buffer the effect of genetic drift and maintain a signature of high genetic connectivity (Lowe and Allendorf 2010). This interpretation is further supported by the very high genetic diversity levels found for common galaxias in both the present study and in previous work on the species (González-Wevar et al. 2015).

The interpretation of effective population size masking population structuring often ignores the close relationship between dispersal and effective population sizes because highly dispersing species typically show higher effective population sizes (e.g. among highly dispersive bird species; Frankham 1997; Arguedas and Parker 2000; Hung et al. 2017). Further, there are many examples of diadromous species exhibiting high effective population sizes while also exhibiting population structuring (Waples et al. 2010; Cushman et al. 2012; Johnstone et al. 2013). Although these studies do not confirm that dispersal and gene flow were unrestricted across our study area for the common galaxias, they do suggest that over large enough time scales (thousands of years), species with high effective population sizes would show structure if they were not connected. As such, there is no direct evidence to suggest that, in this case, high effective population size is masking strong population structuring within common galaxias.

However, it is worth considering that estimating effective population sizes when the exact number of distinct populations or the exact migration rate is not known can lead to biases in the values returned (Ryman et al. 2014). Specifically, in the present study the presence of structuring and limited migration was found for tupong and Australian grayling. The method used to analyse effective population sizes for each species (linkage disequilibrium) is robust provided most genetically distinct populations across the study range have been sampled and migration rates are relatively low (<10%; Waples and England 2011). Although sampling intensity was not sufficient for accurate estimates of effective population sizes within the eastern and western Victorian genetic clusters for all species, reliable estimates for the central Victorian region were determined. Further, sampling within this central region took place over multiple rivers, and clear population structuring was found between this region and the other two genetic regions for tupong and Australian grayling (suggesting the migration rate is likely not high). These factors together suggest that the linkage disequilibrium methods used here are likely to not be so downwardly biased as to affect the overall pattern of effective population size differences, because common galaxias returned estimates up to an order of magnitude larger for the central Victorian region.

Common galaxias and tupong are currently listed as species of least concern under the IUCN status (Bice et al. 2019a, 2019b). By contrast, Australian grayling have been listed as vulnerable species under the Environment Protection and Biodiversity Conservation Act 1999 of Australia. The lower effective population size for Australian grayling is consistent with the reduced range and declining population numbers the species has experienced. In addition, the Australian grayling appears to exhibit a moderate degree of IBD across the spatial scale of the Victorian coast, suggesting population connectivity may be limited. Although the pattern of IBD is not as strong compared with tupong, local-scale management of the species to reflect the less-connected nature of distant river populations is an important consideration. Despite the species exhibiting IBD, it is important to consider there was still genetic connectivity between rivers located over 100 km apart. This suggests that rehabilitation of rivers along the Victorian coast may still be an appropriate strategy for the species to recover and repopulate additional river systems, provided source river populations are not great distances away because, over successive generations, individuals can recolonise. Further, as Australian grayling juveniles have been shown to be attracted to downstream water flow from rivers, targeted water flow offers a promising method to direct individuals to specific rivers and promote connectivity or recolonisation into new rivers (Koster et al. 2021).

Although tupong is listed as of least concern, the species was found to have similar effective population sizes to Australian grayling. Tupong also showed a moderate pattern of IBD, suggesting overall connectivity for the species is limited. Of the three species studied, tupong is likely to be the most vulnerable to human-induced disruptions to connectivity (e.g. dams and barriers blocking passage along a river, or to estuaries), with most recruitment occurring at a local scale for this species. Given these findings, as recommended for Australian grayling, conservation of tupong should consider local-scale management of the species, while further work on measuring population sizes should be undertaken as part of regular management practices. Further, because tupong is likely restricted in inter-riverine connectivity, barriers that block downstream migration and spawning in a river may have significant effects on the population size of that river: any reduction in recruitment success will not be offset by migrants from other rivers to the same degree as a more connected species. Because of this, tupong will likely benefit from local river connectivity restoration and may be an ideal indicator species for the effectiveness of such river connectivity restoration projects (e.g. fishways and environmental flows). Such potential benefits can be seen in targeted environmental releases resulting in a movement increase of almost 40% for tupong. Although tupong is currently listed as a species of least concern, active management may be required to prevent future population decline (Department for Environment and Heritage 2003). One promising avenue of active management includes the use of environmental flows to direct migration towards specific rivers, which has shown some effect on several Victorian diadromous species, although exactly how strong this relationship is for tupong is not fully known (Amtstaetter et al. 2021).

When implementing conservation strategies reflective of connectivity, it is important to consider that genetic connectivity is not the same as demographic connectivity. Although many rivers were shown to be genetically connected, comparatively few migrants are needed to move between rivers to maintain genetic connectivity (Lowe and Allendorf 2010). From a practical perspective, this means that localised recruitment can still be the norm, a possibility that is supported by all siblings identified being found within the same river in this study. As such, it is important to not consider genetically connected populations as a single demographic population that allows dispersal across hundreds of kilometres for single individuals; instead, connectivity may take place over multiple generations.

Connectivity is crucial for effective species management, particularly for riverine and diadromous species. This study shows that species with many similar life history traits, and that use a shared environment, can still vary greatly in the extent of genetic connectivity. This work reveals the complex nature of species-level connectivity, emphasising the importance of directly measuring the population connectivity of target species and using that evidence to directly inform management strategies for different species. The use of genomics to measure genetic connectivity offers an effective way to measure connectivity, providing insights into genetic connectivity that can be used to help guide conservation and management across different species and inform how connectivity takes place at a population level. The value of such analyses is particularly clear for species that disperse in marine environments, because the dispersing larvae are often numerous and difficult to track using direct tracking methods.

The raw genetic data for each species (.csv), and the R code (.R) used to filter each dataset can be found at Open at La Trobe, https://opal.latrobe.edu.au/, under doi: 10.26181/60ac9e770f631.

The authors declare that they have no conflicts of interest.

This work was funded by the La Trobe University ABC Research Funding Scheme 2017 – Securing Food Water and the Environment RFA. James O’Dwyer was supported by the Holsworth Wildlife Research Endowment, The Murray–Darling Basin joint governments Centre for Freshwater Ecosystems PhD scholarship and an Australian Postgraduate Award through La Trobe University. Katherine Harrisson was supported by an ARC Discovery Early Career Researcher Award (DECRA) Fellowship (DE190100636).