Patterns of genetic structuring at the northern limits of the Australian smelt (Retropinna semoni) cryptic species complex

Sampling strategy

A total of 389 individual samples were collected from 15 locations in south-east Queensland, Australia (Fig. 1; Table 1A). Samples were collected using a hand—held seine net from an upstream and a downstream site from each of eight river systems (except the Noosa River—downstream only). We aimed to collect at least 30 individuals per site but this was not always possible, as the species was not abundant in all rivers. Fin clips or entire individuals were placed in 95% ethanol in the field and stored prior to preparation for analysis. All procedures were carried out according to Australian Ethics Commission protocol number ENV/23/14/AEC.

Molecular methods

Genomic DNA was extracted from fin tissue using the DNeasy Blood and Tissue kits (Qiagen, Germantown, MD, USA) following the manufacturer’s directions. Microsatellite markers developed for R. semoni were amplified and genotyped using primers developed by Islam, Schmidt & Hughes (2017). Ten loci were screened across all individuals. The ten loci were BS18, BS3, BS4, BS5, BS20, BS21, BS22, BS24, BS8 and MS24. Microsatellite screening was carried out in 10 µl polymerase chain reaction (PCR) consisting of 0.5 µl of genomic DNA, 0.2 mM reverse primer, 0.05 mM tailed forward primer, 0.2 mM tailed fluorescent tag (either FAM, VIC, NED or PET, Applied Biosystems), 1 × PCR buffer (Astral Scientific, Caringbah, Australia) and 0.02 units of taq polymerase (Astral Scientific, Caringbah, Australia). The following basic thermocycler settings for PCR were performed: initial denaturation at 94 °C for 4 min, followed by 35 cycles at 94 °C for 1 min, 57 °C for 30 s, 72 °C for 1 min and a final extension at 72 °C for 7 min. Fluorescently labelled PCR products were pooled and added to 10 µl of Hi-Di™ formamide with 0.1 µl of GeneScan™ 500 LIZ size standard. Fragment analysis was conducted on an ABI PRISM 3130 Genetic Analyzer (Applied Biosystems, Foster City, CA, USA) according to the manufacturer’s instructions. Data were scored using GENEMAPPER version 3.1 software (Applied Biosystems, Foster City, CA, USA).

Two individuals from each of the 15 populations represented in the microsatellite study were randomly selected for mtDNA analysis except the Mooloolah_upper (MLU) site from which one additional sample was included for analysis. Samples from four additional sites not included in microsatellite analysis were also sequenced—two sites from the Mary River and two sites from the Brisbane River (Fig. 1 and Table 1B). In total 68 individuals from nineteen sites were sequenced. A 666 bp fragment of the cytochrome b region of the mtDNA genome was selected for sequence analysis. The primers HYPSLA and HYPSHD (Thacker et al., 2007) were used to amplify the region in 10 µL reaction mixtures. PCR conditions were 4 min at 95 °C, followed by 45 cycles of 30 s at 95 °C, 45 s at 53 °C, 45  s at 72 °C and a final extension cycle of 7 min at 72  °C. MtDNA sequences were edited and aligned using Geneious version 9.1.5 (Kearse et al., 2012).

Data analysis

Genetic variability and levels of differentiation

After Bonferroni correction, 3 out of 15 populations exhibited deviations from HWE in only one or two loci (Table S1). All loci were kept for further analyses since deviations were not consistent across populations. Instances of null alleles estimated using MICRO-CHECKER were rare and not consistently associated with specific loci or populations (Table S2). No evidence for genotypic linkage disequilibrium was observed between any pair of loci across all populations.

Population genetic diversity indices are shown in Table 1A. Microsatellite genetic diversity was high. Mean number of alleles per population ranged from 4.60 (MLD) to 14.70 (CMU) and heterozygosity averaged across loci ranged from 0.566(MLU) to 0.887 (CMD and CMU). Most sites exhibited positive F_IS values, indicating that most of the populations had slight heterozygote deficits.

Most of the pairwise F_ST values between the 15 populations were significant and ranged from −0.018 to 0.404. The SEQ-N populations were more diverged from one another than the populations in the SEQ-S group. The lowest pairwise F_ST value (F_ST =  − 0.018; P < 0.05) was observed between populations NRD and NRU. The highest genetic divergence (F_ST = 0.404; P < 0.05) was observed between populations NSD and MLU. Out of 105 comparisons, only two comparisons were non-significant (P > 0.05) and each of these pairs was from within the same river (Logan-Albert and Nerang). However, one site of each of these pairs contains a small number of samples which probably reduces the power to detect the significant F_ST values. Generally F_ST comparisons revealed much less divergence among populations within the same river than between populations from different rivers (Table 2).

STRUCTURE analysis revealed the highest likelihood was at K = 8 clusters (Average log probability of data Ln[P(DK)] ± SD = − 15125.4 ± 0.584618) (Fig. 2, Table S3) indicating this as the best estimate of the true number of genetic clusters. The height of ΔK was used as an indicator of the strength of the signal detected by STRUCTURE (Evanno, Regnaut & Goudet, 2005). ΔK also showed the highest peak at K = 8, suggesting eight genetically homogeneous clusters across the sampled populations (Fig. S1) and negligible migration was observed among rivers (Fig. 2).

Strong genetic structure was observed by AMOVA. The AMOVA showed significant genetic differentiation between the two groups (SEQ-N and SEQ-S) (F_CT = 0.05), but also among populations within groups (F_SC = 0.19) (Table 3A). There were similar patterns observed in the two groups when each group (SEQ-N and SEQ-S) was analysed separately, with the F_CT value (among rivers) higher than the F_SC value (among sites within rivers) in both groups (Table 3Bi and 3Bii). However, the overall F_ST values, and each of the other F statistics in the hierarchy were slightly higher within the SEQ-N group than the SEQ-S group.

Isolation-by-distance

There was a significant correlation between genetic differentiation and stream distance among populations in the SEQ-S group (R² =0.3687, p = 0.001; BRD, BRU, LAD, LAU, CMD, CMU, NRD, NRU, CRD, and CRU) (Fig. 3A), but not for the SEQ-N group (R² =0.0355, p = 0.302; MRD, NSD, MRU, MLD and MLU) (Fig. 3B).

Contemporary migration

Very little contemporary migration was observed among the coastal river populations. Only six sampled populations contained individuals that were identified as potential immigrants from the BAYESASS analysis. In all cases, the putative source population was the paired site within the same catchment (Table 4). In most cases, dispersal was from the upstream to the downstream site. Only individuals from Currumbin Creek were estimated to have dispersed in an upstream direction. Only thirty-six (<10%) of 389 individuals across all sites were identified as F₀ migrants using the “detection of first generation migrants” option in GENECLASS2 (Table 5).

MtDNA sequences analysis

The edited alignment for the cyt b gene was 575 bp and included 121 variable positions. All sequences are lodged under GenBank accession numbers MG867590 –MG867657. The neighbour-joining tree revealed two strongly supported clades (bootstrap 89% SEQ-S; 96% SEQ-N; Fig. 4). Phylogeographic structure was also clearly evident within clades. All individuals from four sites in the Brisbane River formed a distinct clade, and all three rivers sampled for the SEQ-N lineage formed shallow clades (i.e., Mary, Noosa and Mooloolah rivers; Fig. 4). Genetic distance was high between northern smelt lineages and the southern smelt sequences used as outgroups (uncorrected mean nucleotide distance 0.15–0.17). The mean nucleotide distance between the two northern lineages (SEQ-N and SEQ-S) was 0.04 (SE = 0.007).

Population structure and dispersal

Based on previous studies of Australian smelt (Woods et al., 2010; Hughes et al., 2014), we had hypothesized that R. semoni in south-east Queensland would exhibit limited genetic connectivity among river systems. Our findings of strong genetic differentiation among rivers support this hypothesis. For both groups (SEQ-N and SEQ-S), there were highly significant F_ST values, which indicated that populations were not panmictic. The results of the significant genetic differentiation among rivers were also consistent with the results of population structure revealed from the STRUCTURE analysis, suggesting restricted gene flow and limited dispersal among populations in both groups. Limited dispersal was supported by our first-generation migrant detection analysis in Geneclass2, which demonstrated that less than 10% of individuals in each population were immigrants.

The sample from Tinana Creek (MRD site), was differentiated from the rest of the populations in the SEQ-N group (Table 2). This might be the result of a barrier which separates Tinana Creek from the rest of the Mary River system despite their close proximity to one another. Tinana Creek runs into the Mary River not far from the mouth, with both drainages having tidal estuarine reaches in the lower sections. The differentiation of the Tinana Creek population from the main stem of the Mary River is also observed in a number of other freshwater species including Mary River Cod, Maccullochella mariensis (Huey, Espinoza & Hughes, 2013), Mary River Turtle, Elusor macrurus (Schmidt et al., 2018), freshwater crayfish Cherax dispar (Bentley, Schmidt & Hughes, 2010) and Australian lungfish Neoceratodus fosteri (Hughes et al., 2015).

Although populations in the SEQ-N group were slightly more highly structured than those in the SEQ-S group, fishes in both groups exhibited restricted gene flow. The high genetic structuring in SEQ-N populations might be the result of genetic divergence within the lineage that occurs over a very small geographical range. The lower levels of genetic structure of R. semoni populations in the SEQ-S group than the SEQ-N group suggests that more current or recent gene flow occurs within and among the catchments in SEQ-S group than for the SEQ-N group. This could be because the catchments in the SEQ-S group have been connected more recently and/or more often than those in the SEQ-N group. This could be because during extreme flood events, SEQ-S catchments might be connected by freshwater plumes flowing into Southern Moreton Bay. Another plausible reason might be that these cryptic species (SEQ-N and SEQ-S) have different microhabitat preferences, with SEQ-N being less tolerant of brackish conditions (Hughes, 2015).

An important model for stream dwelling species is isolation by distance (IBD). In this model, equilibrium between genetic drift and gene flow may be reached in species where the life time dispersal distance is less than the range. Here, a relationship between stream distance and genetic differentiation should be evident (Wright, 1943). In this study, a strong IBD relationship was identified among the SEQ-S populations, but not among SEQ-N populations. This suggests that for SEQ-S populations, dispersal, when it occurs, is more likely between nearby catchments. Similar IBD relationships have been reported for other coastline restricted species (Keenan, 1994; Jerry & Baverstock, 1998; Shaddick et al., 2011; Schmidt et al., 2014). Lack of IBD for the SEQ-N group may be attributed to an insufficient number of population samples available for comparison and/or the greater degree of population isolation within this group relative to the SEQ-S group, consistent with the overall higher F_ST estimates among SEQ-N populations. Hughes et al. (2014) observed similarly contrasting patterns of population genetic structure between cryptic species groups of southern Australian smelt. In that study, two informal species groups (MTV and SEC) with adjacent distributions along the western and eastern coast of southern Victoria had microsatellite-based F_ST values of 0.19 and 0.07 respectively (Hughes et al., 2014). Using otolith microchemistry, Hughes et al. (2014) also showed that the more structured western group (MTV) had a greater proportion of nondiadromous populations relative to the weaker structured eastern group (SEC). The similar pattern of contrasting structure observed here between northern groups in the Australian smelt complex (SEQ-N and SEQ-S), is probably not due to differences in diadromous behaviour because preliminary evidence from otolith chemistry suggests all of these populations are nondiadromous (R Islam, 2017, unpublished data). Higher structuring of the SEQ-N group could possibly be due to genetic drift if these populations have been established for a longer period of time at the northern-most limit of Australian smelt distribution relative to the SEQ-S populations.

The complementary pattern of population structuring in both microsatellite and mtDNA data between the SEQ-N and SEQ-S groups could reflect long-term phylogeographic isolation or perhaps further evidence for cryptic species within Australian smelts as first highlighted by Hammer et al. (2007). Mean cyt b divergence of 4% between SEQ-N and SEQ-S samples is close to the 3.6% divergence observed for the full mitochondrial molecule reported by Schmidt, Islam & Hughes (2016), and is within the range of lineage divergence reported for R. semoni in southern Queensland (Page & Hughes, 2014). It should be noted that Schmidt, Islam & Hughes (2016) incorrectly attributed mitochondrial sequence data for the SEQ-N group to the CEQ group (i.e., Genbank accession KX421785 corresponds to SEQ-N and KX421784 corresponds to SEQ-S). The level of cyt b divergence between the SEQ groups and southern lineages of R. semoni is very large (15–17%) and adds to previous studies that have highlighted the likely existence of a cryptic species complex within the taxon currently referred to as R. semoni (Hammer et al., 2007; Hughes et al., 2014).

Contemporary migration

The Bayesian assignment analysis detected contemporary movement of individuals only between proximate sites within rivers (Table 4). Contemporary dispersal was not observed between rivers. Although, most of the sites that we sampled were within 10–60 km stream distance of another sampled site, there was no contemporary dispersal among the majority of those rivers in either group. The present results coincide with the findings of some southern smelt populations, where contemporary movement was observed among populations within a catchment (Woods et al., 2010). The species appears to occur in a wide range of freshwater habitats, many of which are isolated by long stretches of unfavourable habitat. Our data therefore suggest that if local extinctions occur in one or more of these streams within a reach of the river, then recolonization from elsewhere is unlikely to occur rapidly.