Phylogeography of the freshwater glassfishes (Ambassidae) of eastern Australia: cryptic species and hybrid zones

Sample collection

All collections were obtained under permit from Queensland and New South Wales Fisheries and research was conducted under approval from the University of Canberra Institutional Animal Care and Use Committee (approvals: CEAE 13-06, CEAE 14-17, CEAE 15-05, IACUC 0447, IACUC 9371). Snap-frozen (liquid nitrogen) or ethanol-fixed tissues were deposited in the South Australian Museum’s Australian Biological Tissues Collection. We sampled A. agassizii across its range, except for populations north of the Mulgrave River, which we were unable to obtain. We also included fish from the Bulloo River (Lake Eyre Basin). These have never been formally identified but are clearly genetically related to the lineage containing A. agassizii on the basis of our broader unpublished data. We also sampled other Ambassis populations from all adjacent river basins to the west of the main study area and further north through Cape York. Sample sizes per site were usually low (typically two individuals per population, range: n = 1–10), as is appropriate for studies that primarily focus on species boundaries and broader phylogeography (Adams et al. 2014). Although each molecular dataset employed a different number of populations and individuals, there was considerable site overlap between marker types, thus ensuring that all key lineages were always included. In total, 115 sites (Fig. 2) were surveyed for at least one dataset, including 239 individuals for SNPs, 146 for allozymes and 101 for mtDNA (Table 2).

Fig. 2.

Map showing the geographic location of all sites and their final taxon designation by using the symbols provided. Sites numbers follow those used in Table 2.

Allozyme analyses

Because this study was initiated before low-cost genomic screening became routinely available, our initial genetic assessment involved a comprehensive allozyme screen of 146 individuals from 72 sites, with these encompassing the entire geographic range for which frozen tissues were available prior to 2017. Allozyme electrophoresis of muscle homogenates was conducted on cellulose acetate gels as detailed by Richardson et al. (1986). All laboratory protocols, enzyme details and the methodological procedures for undertaking principal co-ordinates analysis (PCA) and calculating genetic distances or fixed differences have been published previously (Richardson et al. 1986; Hammer et al. 2007; Adams et al. 2014).

SNP genotyping and filtering

Considering the allozyme and cytb results, 239 individuals (from 103 sites) were ultimately selected for SNP genotyping (Table 2). DNA was extracted and sequenced by diversity arrays technology (DArT Pty Ltd, Canberra, ACT, Australia; see www.diversityarrays.com) on the basis of a combination of complexity reduction by using two restriction enzymes (PstI, recognition sequence 5′-CTGCA|G-3′, and SphI, 5′-GCATG|C-3′), implicit fragment size selection, and next-generation sequencing (Kilian et al. 2012). All other details of the sequencing methods are presented in Georges et al. (2018).

The SNP data and associated metadata were read into the package dartR (ver. 2.9.4, see https://cran.r-project.org/package=dartR; Gruber et al. 2018) as a genlight object (sensu adegenet, ver. 2.1.10, see https://cran.r-project.org/package=adegenet; Jombart 2008) for subsequent processing with dartR. As detailed elsewhere (Unmack et al. 2023), preliminary analyses to identify individuals that had either excessive missing data (likely owing to poor quality DNA) or excessively high heterozygosity (likely reflecting some tissue or DNA cross-contamination) resulted in the removal of four individuals from eastern coast basins (high heterozygosity counts). We then generated the ‘primary’ dataset (32,010 SNPs for n = 239 fish) by filtering for reproducibility ≥0.99 and removing any monomorphic loci (this latter filter being applied at every level of data reduction). This ‘primary’ dataset was then further filtered to create two initial working datasets, namely (a) the ‘PCA’ dataset (4057 SNPs and 3.71% missing data after filtering for a per-locus missing data threshold ≥0.9 and removing secondaries using the ‘best’ option) for the analyses most sensitive to the presence of missing data and linked loci (PCA and fixed difference calculations), and (b) the ‘Phylo’ dataset (204 individuals, 14,209 SNPs and 12.42% missing data after filtering by locus ≥0.7 and removing individuals that were considered highly introgressed) for the tree-based analyses. We also repeated each relevant filtering step with the same thresholds for any subset of individuals derived from the working ‘PCA’ and ‘Phylo’ datasets.

Detecting admixture

Given that the inclusion of even a few admixed individuals can distort phylogenetic trees and confound assessments of diagnosability (Unmack et al. 2022, 2023), we followed our standard analytical framework to distinguish primary lineages from admixed individuals and sites. As detailed in Unmack et al. (2022), primary lineages can be distinguished from ‘hybrid’ individuals, sites or subpopulations by differences in their genetic profiles when subjected to the combined application of stepwise ordination by PCA, observed heterozygosity counts for both groups and individuals (using the primary SNP dataset), phylogenetic trees, and assessments of diagnosability based on pairwise fixed differences (FDs). Support for the presence of any putative sites or zone of admixture was then further assessed by comparing their geographic location with those of their putative parent lineages.

We also used STRUCTURE (ver. 2.3.4, see https://web.stanford.edu/group/pritchardlab/structure.html; Pritchard et al. 2000) to offer an independent perspective on lineage diagnosability and the extent and nature of admixture among all relevant lineages. On the basis of stepwise PCA and the phylogenetic trees, two subsets of the ‘PCA’ dataset were assembled, the first for all individuals referable to A. agassizii and the second for all other primary lineages plus the ‘CEQ’ phylogroup of A. agassizii (implicated by PCA as a source lineage for some northern admixed sites) plus all putatively admixed fish. Each of these subsets was further filtered using a per-locus missing data threshold of ≥0.98, read depth range of 8–50, and global per-locus minor allele frequencies of <0.015 (gl.filter.maf(gl,threshold = 0.015)), and then converted to a STRUCTURE input file in dartR. STRUCTURE was run using the following parameters: 5 replicates each for K = 1–7, admixture model, correlated allele frequencies, burn-in = 50,000, and iterations = 200,000. We then used KFinder (ver. 1, see https://www.zsl.org/about-zsl/resources/software/kfinder; Wang 2019) to provide estimates of the optimum K by using both the ΔK method (Evanno et al. 2005) and Wang’s parsimony method. Barplots for various values of K were created using CLUMPP (ver. 1.1.2, see https://rosenberglab.stanford.edu/clumpp.html; Jakobsson and Rosenberg 2007).

Ordination of SNPs

Ordination was achieved using PCA, with individuals as entities and loci as attributes as implemented by the gl.pcoa and gl.pcoa.plot functions of dartR. Follow-up PCAs were performed on various major groups and subgroups on reduced and re-filtered datasets by using the same protocols.

Phylogenetic trees

Genomic sequences were constructed by concatenating only the variable bases from each SNP locus into a single partition for each individual. Loci for which a SNP was not called were designated as an ‘N’ to maintain sequence alignment. Phylogenetic analyses were conducted by maximum likelihood (ML) using IQ-TREE (ver. 1.6.12, see http://www.iqtree.org/; Minh et al. 2020) run on the CIPRES server (Miller et al. 2010), incorporating IQ-TREE’s model selection procedure (-m TEST; Kalyaanamoorthy et al. 2017), resulting in the selection of the TVM + F + G4 model followed by 10,000 replicates of ultrafast bootstrapping with the SH-like approximate likelihood-ratio test (-aLRT 1000; Hoang et al. 2018).

Fixed difference analyses for SNPs

The diagnosability of primary lineages or phylogroups and any putative admixed populations was assessed on the ‘PCA’ dataset by counting the pairwise number of SNPs displaying fixed or nearly fixed allelic differences (Unmack et al. 2022). This was achieved using the gl.fixed.diff function in dartR (t = 0 for an absolute difference and t = 0.025 or t = 0.05 to allow either 5 or 10% tolerance for shared alleles at a locus). We routinely use nearly fixed differences in recognition of the fact that low levels of hybridisation are very common among sister species of freshwater fishes (Verspoor and Hammart 1991) and to echo the widespread use of partially diagnostic characters in morphological taxonomy (and indeed for Ambassis in particular; Allen and Burgess 1990).

mtDNA genotyping and analyses

We ultimately sequenced 101 individuals from 53 sites for the mitochondrial cytb gene following the polymerase chain reaction (PCR) protocols in Hammer et al. (2014), except for the following primer details. Most samples were amplified for the entire gene by using the primers Glu31 TGRCTTGAAAAACCACCGTTGT and AM.Thr.28 CCGTAACTGGATTACAAGACCG or AM.Thr.32 CCTTCGACATTCGGCTTACAAG. When amplification was insufficient, the gene was amplified in two halves by using Glu31 and HD.Hyps GGGTTGTTGGAGCCSGTTTCGT and PP.505 TCAGTAGACAACGCCACCCT and AM.Thr.28 or AM.Thr.32. Sequences were edited using Chromas (ver. 2.6.5, Technelysium, Tewantin, Qld, Australia, see https://technelysium.com.au/wp/chromas/) and imported into BioEdit (ver. 7.2.5, see https://bioedit.software.informer.com/; Hall 1999). Sequences were aligned by eye and checked by amino acid coding in MEGA (ver. 11.0.13, see https://www.megasoftware.net) to test for unexpected frame shift errors or stop codons. Data were analysed phylogenetically by using IQ-TREE (ver. 1.6.12) as per the SNP analysis, by using the model HKY + F + G4.

Lineage, taxon and phylogroup nomenclature

Individuals, groups and lineages are designated throughout by a unique symbol matching their final operational taxonomic unit (OTU) membership, matched with an informal nomenclature, whereby each primary lineage or phylogroup is designated a three-letter code and putative admixed groups by a combination of the codes for their putative parents. Thus, for A. agassizii we delineate four primary phylogroups: central-eastern Queensland encompassing the type locality (CEQ), southern Queensland ‘northern phylogroup’ (SQN), southern Queensland ’southern phylogroup’ (SQS), and Murray–Darling Basin (MDB), plus one admixed site, SQN × SQS. The remaining primary lineages (all genetically distinct from A. agassizii) were Bulloo and Burdekin basins (BUL), northern coastal Queensland (NCQ), and Cape York and Cooper Creek (ASP), with each found to be involved in admixture events involving one or more of the abovementioned OTUs.

Allozyme dataset

We were able to assign allozyme profiles at 57 putative loci for the 146 fish screened. Allozyme profiles for all key OTUs and admixed groups are presented in Supplementary Table S1. An initial PCA (Fig. 3a) showed a complicated pattern of genetic affinities among the sites sampled, with much but not all reflecting the geographic arrangement of sites (Fig. 2, Table 2). Sites south, west and south-west of Site 37 fell into one of three broad but discrete clusters, namely a widespread A. agassizii cluster, a cluster of three inland sites (25, 30 and 94) spanning ~900 km from the western edge of Burdekin Basin into the Bulloo Basin (lineage BUL), and a cluster that included widely dispersed sites from west and north of the Bulloo (lineage ASP) plus a number of coastal sites from the Ross Basin northward (Sites 3–8, and 19, lineage NCQ). There was also qualitative evidence of admixture between the BUL and ASP–NCQ clusters, with six individuals placed as intermediate between these two clusters. Three additional individuals from Site 35, ultimately shown by their SNP profiles to be admixed, were placed in the position consistent with their final taxon allocation (NCQ(CEQ)).

Fig. 3.

Scatterplots of PCA scores in the first two dimensions for the allozyme dataset (relative contribution of each dimension in parentheses). (a) Initial PCA on all 146 individuals. (b) Follow-up PCA on the 99 individuals referable to A. agassizii. (c) Follow-up PCA on the 33 individuals belonging to either NCQ or ASP. Symbols and taxon labels are as given in Fig. 1.

Given that all but five of the allozyme sites were also included in the SNP dataset (Table 2), we conducted only two additional explorations of genetic structure in deeper PCA dimensions. Surprisingly, even though a follow-up PCA restricted to the A. agassizii cluster (Fig. 3b) largely supported the presence of four geographically based phylogroups (CEQ, SQN, SQS and MDB) plus limited admixture between phylogroups SQN and SQS (at Site 70), it also inferred that the primary phylogeographic dichotomy in this species involves phylogroup SQS (occupying just ~10% of its geographic range) v. CEQ–SQN–MDB combined. This primary dichotomy is also diagnosable by two fixed differences (Table S2). The other follow-up PCA on the ASP–NCQ cluster (Fig. 3c) found modest support for distinguishing these two lineages (which are also diagnosable by 1 FD; Table S2) and further inferred a phylogeographic split in ASP between the three most-northerly sites (95–97) v. those from Site 101 south.

Ordination of SNP dataset

While numerous targeted PCAs were ultimately undertaken to identify the putative parental lineages for all admixed OTUs, their main outcomes can be summarised by the four presented herein (Fig. 4a–d). The initial PCA of all 239 fish (Fig. 4a) unequivocally found the same three primary clusters inferred for the allozyme dataset, i.e. A. agassizii, the inland taxon BUL, and a ‘composite’ cluster of all other sites to the west or north. Within taxon BUL, a primary phylogeographic break was evident between the Bulloo and Burdekin sites in the first dimension (Fig. 4b), with additional structure detected within the Burdekin (Belyando v. upper Burdekin).

Fig. 4.

Scatterplots of PCA scores in the first two dimensions for the SNP dataset (relative contribution of each dimension in parentheses). (a) Initial PCA on all 239 individuals. (b) Follow-up PCA on taxon BUL (43 fish, 1760 SNPs, 1.92% missing data). (c) Follow-up PCA on A. agassizii (98 fish, 1544 SNPs, 4.3% missing data). (d) Follow-up PCA after removing A. agassizii and BUL (98 fish, 3335 SNPs, 4.2% missing data). Symbols and taxon labels are as given in Fig. 1.

Stepwise PCA on the A. agassizii cluster (Fig. 3c) again supported the allozyme results, namely four distinct phylogroups (CEQ, MDB, SQN and SQS) plus a single admixed SQN × SQS site, but was less emphatic about a primary split between SQS and the other three phylogroups. Interestingly, while most river catchments contained a single phylogroup (Table 2), the Burnett harboured three (all except SQS), including a geographic-outlier population of CEQ (Site 52).

Last, a follow-up PCA on the composite cluster (Fig. 4d) confirmed the taxa ASP and NCQ as the two most distinctive groups in the first dimension and displayed a diverse assortment of other geographically intermediate clusters in both dimensions. The putative parental affiliations of all these likely admixed OTUs were initially assigned by comparing their scatterplot positions in both the initial and follow-up PCAs with the PCA-placement of their parental lineages. Ultimately, we postulated five likely admixed OTUs in this composite cluster, two from a single site (ASP × NCQ and NCQ × CEQ) and three across multiple sites (BUL × NCQ × ASP at four sites, NCQ × BUL from eight sites, and NCQ(CEQ) from two sites).

Distinguishing primary lineages from admixed OTUs

Reflecting both nuclear genetic datasets, ordination resulted in a working hypothesis of seven pure lineages or phylogroups, nested within four candidate taxa, and six putatively admixed OTUs for the Ambassis sites surveyed herein (summarised in Table 3). Several additional lines of evidence support this primary taxonomic scaffold. First, diagnosability assessments among all OTUs show the presence of numerous fixed and nearly fixed differences among all four candidate taxa (e.g. range 77–463 for SNPs allowing 10% tolerance; Table 4), including some involving allozyme loci (Table S2). This contrasts markedly with the comparatively low number of FDs distinguishing the two admixed OTUs with sufficient sample sizes (sites > 3 and N > 10; BUL × NCQ × ASP, NCQ × BUL) to provide reasonably robust FD comparisons with their parental taxa (range 0–28 across all tolerance levels; Table 4). Second, admixed OTUs mostly displayed considerably higher heterozygosity counts than those found for every parental lineage (Tables 3 and S1 for SNPs and allozymes respectively), the one exception being SQN × SQS, which had higher levels of heterozygosity for allozymes but not for SNPs than did lineage SQN. Third, the geographic placement of every admixed OTU is broadly consistent with the presence of a ‘hybrid’ zone between its putative parental lineages (Fig. 2). Fourth, individuals assigned to these admixed OTUs all occupied early branching positions in the preliminary ML tree constructed from the primary SNP dataset (Supplementary Fig. S1). Finally, these same SNP-based diagnosability assessments did offer some additional support for a primary split within A. agassizii between lineage SQS v. the rest (Table 4), as flagged by the allozyme analyses.

Exploring further instances of admixture between lineages and phylogroups

The results of the STRUCTURE analyses on A. agassizii are presented in Fig. 5a for K = 4. Given the disparity between the optimum value of K suggested by the ΔK and parsimony methods (3 v. 6), we chose K = 4 to match the number of primary groups shown by ordination. In addition to correctly classifying the two SQN × SQS individuals as an ~50:50 mix of its two parental phylogroups, there was evidence of modest to moderate admixture between all abutting phylogroups. With one exception, the sites involved were geographically close to their nominated ‘donor’ phylogroup (Fig. 2, 5a), with the most prominent examples being CEQ Sites 45–48 (close to and admixed with SQN), SQN Sites 56–58 (adjacent to CEQ, close to MDB, and admixed with both), SQS Sites 64–69 (surrounded by and modestly admixed with SQN), plus MDB Sites 49–51, 53, 54, 81 and 82 (modestly admixed with and close to CEQ or SQN). Even the single exception, MDB Site 84 in the Warrego, was plausibly located much closer to its CEQ admixture parent than to any other phylogroup.

Fig. 5.

STRUCTURE plots of K = 4 for each of two subsets of the Ambassis SNP dataset. Taxon names, symbols and lineage colours match those used in Fig. 1. Each individual in the STRUCTURE plot is labelled by site number (following Table 2, Fig. 1). (a) Plot for all A. agassizii (98 individuals, 384 SNPs after stringent filtering). (b) Plot for individuals referable to lineages BUL (excluding Site 94 in the Bulloo Basin), CEQ (excluding Site 52), ASP, NCQ and their putative admixed sites (149 fish, 1011 SNPs).

Fig. 5b presents the STRUCTURE plot for K = 4, focusing on the other primary Ambassis lineages plus CEQ. This value of K matches the number of lineages shown by ordination and was also the lower optimum value for the ΔK and parsimony methods (4 v. 5). Note that to better focus on the geographic region where all the targeted taxonomic complexity occurs, we excluded the genetically divergent BUL population from outside the region (Site 94, Fig. 2, 4b). As predicted by other analyses, the five admixed OTUs NCQ × CEQ, NCQ × ASP, NCQ (CEQ), NCQ × BUL, and NCQ × BUL × ASP were each successfully classified consistent with their expected parental lineages. Of the four primary candidate taxa, all except NCQ showed modest evidence of admixture at some sites with other taxa in the region. Generally, the sites involved (particularly BUL Sites 27, 33 and 34 plus ASP Sites 111 and 114) were either geographically proximate or located in adjacent river catchments, thus further reinforcing the notion that the Ambassis species of eastern Australia are likely to have semi-permeable species boundaries.

Phylogenetic affinities of primary lineages and phylogroups

The removal of individuals initially classified as belonging to admixed OTUs followed by refiltering produced a final ‘phylo’ dataset of 204 individuals and 14,209 SNPs. ML recovered one tree with a –ln score of −263471.385 (Fig. 6; full tree for all 239 individuals in Fig. S1). Support across most of the deeper nodes of the tree was strong, with A. agassizii being the sister group to BUL and NCQ being the sister group to ASP, with these two pairs being deeper sister groups. As would be expected from their broad distribution, A. agassizii displayed deep phylogeographic structure. However, none of its four PCA phylogroups was resolved as a monophyletic lineage, a result that is consistent with the high degree of among-phylogroup admixture observed in the STRUCTURE analysis (23 of 53 sites displayed >15% admixture involving two or more phylogroups; Fig. 5a). Phylogenetic relationships within A. agassizii appeared particularly confounded by the presence of admixed populations primarily centred on the Burnett Basin and the rivers immediately north (Site 48) and south (Site 58). As such high levels of reticulate evolution are known to effectively render such datasets unsuitable for producing bifurcating trees (Morrison 2016; Unmack et al. 2023), we generated a second ML tree for the 55 individual A. agassizii (from 30 sites) that were classified by STRUCTURE as ‘pure’ for one of the four phylogroups. However, even this restricted tree (Fig. 6b) supported monophyly for only two of the four ‘phylogroups’ (MDB and SQS).

Fig. 6.

Phylogenetic trees for the SNP dataset. Individuals are labelled by site code + River Basin (as per Table 2) + unique fish number. Bootstrap proportions of 100% are indicated by an asterisk for the earlier-branching nodes. (a) Final ML tree, rooted at the midpoint, for 204 Ambassis individuals following removal of the 35 fish identified as obviously admixed (Fig. 4, 5). (b) ML tree for the 55 A. agassizii individuals classified by STRUCTURE as ‘pure’ for one of the four phylogroups (represented by the symbols used in Fig. 2).

As for A. agassizii, the presence of numerous admixed individuals in the ML tree precludes any detailed assessment of the admittedly shallow phylogeographic structure within both BUL and NCQ. By contrast, taxon ASP displayed a well-supported phylogeographic dichotomy between the Cooper sites and those in more-northerly catchments and a clear geographically aligned pattern of genetic relatedness in the latter phylogroup.

Phylogenetic analysis of mtDNA

The cytb dataset contained a smaller subset of 101 individuals for 1141 base pairs. ML recovered one tree with a -ln score of −3119.278 (Fig. 7). The cytb dataset presented quite a different topology from and lower support at deeper nodes than the SNP and allozyme perspectives. The first branching lineage consisted of a subset of ASP from populations on the northern portion of Cape York (Sites 95–97 from the Archer River north). The next sequential branching lineage was NCQ, followed by a second distinctive ASP lineage, representing all other sites sampled. The remaining two taxa, A. agassizii and BUL, were sisters, although only separated by very short branch lengths. Apart from the presence of two distinctive paraphyletic mtDNA lineages with ASP, representing a major north–south phylogeographic dichotomy based between Sites 97 and 98, there was minimal clear genetic structure present within any other taxon.

Fig. 7.

Maximum likelihood tree for the 101 cytb sequences of eastern Australian Ambassis. The tree was rooted using A. elongata as an outgroup (GenBank sequences KF414820.1 and KF414813.1). Individuals are labelled as in Fig. 6. Bootstrap proportions are indicated for the earlier-branching nodes by asterisks (**100%, *98–99%).

Candidate species

Both nuclear datasets were broadly concordant in delineating four primary Ambassis taxa across the 115 sites sampled in eastern Australia. The great majority of sites (n = 99; Table 1) appeared to harbour only non-introgressed populations of one of these four taxa, each genetically diagnosable from its congeners by numerous absolute fixed differences for SNPs (range = 37–265; Table 4) plus supported by fixed allozyme differences (Table S2). Such levels of diagnosability compare favourably with the total number (~10) of morphological characters available for the taxonomic revision of Allen and Burgess (1990), particularly given that few of these characters were shown to be unequivocally diagnostic for members of the A. ‘agrammus’ complex. Non-introgressed populations of each taxon occupied allopatric distributions (Fig. 8), with geographic ranges varying from widespread and broad (A. agassizii and ASP), to widespread and narrow (BUL), to restricted (NCQ). However, disrupting this simple taxonomic framework was our finding of obvious among-taxon admixture involving various combinations of either two or three taxa at 15 sites, all occupying geographically intermediate positions between their primary parental taxa. Three putative hybrid zones were evident, one between A. agassizii and NCQ, a second between BUL and NCQ, and a third involving BUL, NCQ and ASP. In addition, a single F₁ hybrid between NCQ and ASP was detected at a site broadly consistent with being in an overlap zone between these two taxa.

Fig. 8.

Map depicting the distribution by river catchment of non-introgressed populations for each of the four taxa, plus the general location of all hybrid zones. Also shown are the names of the river basins (black text) and individual rivers (grey or white text) referred to in the Discussion.

Prior to the ready availability of genomic data, it was widely believed that hybridisation and introgression were rare phenomena in most animal groups (Dowling and Secor 1997; Taylor and Larson 2019), although even then fish were regarded as one of the vertebrate groups most likely to display between-species hybridisation (Verspoor and Hammart 1991; Scribner et al. 2000). However, it now well established that reticulate evolution is pervasive in most animal groups to the point where many have argued that species boundaries ought to be considered semi-permeable and that speciation can occur even in the presence of ongoing gene flow (Edwards et al. 2016; Harrison and Larson 2016; Mallet et al. 2016; Campbell et al. 2018; Taylor and Larson 2019; Bangs et al. 2020; Walker et al. 2024), particularly when habitats have undergone rapid alteration as a result of anthropogenic modification or climate change (Scribner et al. 2000; Banerjee et al. 2023; Brauer et al. 2023; Zbinden et al. 2023).

Acceptance of these realities does of course not remove the key dilemma confronting those attempting to construct taxonomic frameworks for a group of related organisms. Put simply, how do we distinguish gene flow between distinctive phylogroups within a single species from between-species hybridisation accompanied by levels of introgression that remain under the threshold required to maintain species boundaries? Clearly, both alternate scenarios are potentially in play for the Ambassis taxa studied herein and neither can be unequivocally chosen without more detailed biological investigations assessing the morphological and genomic integrity of both primary parental lineages and admixed individuals in and around these hybrid zones. Numerous studies of hybrid zones have demonstrated that, outside of and even within their zone(s) of hybridisation, parental species can maintain their morphological and genomic integrity, provided ‘differential’ introgression occurs, i.e. some key regions of the genome remain immune from reticulate evolution (Harrison and Larson 2016). The persistence of such ‘islands of genomic divergence’, despite the occurrence of widespread introgression involving other genomic regions, has been observed in fish (e.g. Souissi et al. 2018; Sotola et al. 2019; Du et al. 2024) as well as in many other animal and plant groups (e.g. Harrison and Larson 2014; Walsh et al. 2016; Yang et al. 2020; Hagberg et al. 2022; Bruni et al. 2023).

Although we hesitate to conclude herein that all four Ambassis taxa are definitely valid candidate species, several observations based on the available molecular data do support their recognition as likely candidates. In addition to being readily diagnosable at numerous independent genomic markers, the combined geographic distribution of the four primary taxa vastly overshadows the total area occupied by all observed hybrid zones (>97 v. <3%; Fig. 8). Second, two of the four taxa are directly (A. agassizii and A. agrammus = ASP) referable to an existing described species, although the perceived ranges of the following two warrant a reappraisal of their morphological distinctiveness on the basis of the updated geographic ranges presented here (BUL may or may not represent Ambassis sp. NW, and there would be no available name for NEQ). Third, all four taxa were characterised by distinct matrilineal clades (two paraphyletic clades for ASP and one per taxon for the others). Moreover, all admixed sites displayed only the single cytb clade expected for the dominant parental ancestor, suggesting that hybridisation is non-random with respect to maternal ancestry, as is often found for species that only partially introgress (Toews and Brelsford 2012). Under these perspectives, we suggest that the hypothesis of four candidate species is a sensible starting point for all future phylogenetic and morphological assessments.

Patterns of introgression

Most of the introgression found occurred in the Burdekin Basin and adjacent basins to the north and south (see Fig. 8 for all references to catchment and river names). This region represents the intersection among three candidate taxa, which are likely to have a complex history of exchanges. The Herbert Basin was dominated by NCQ, but with differing admixture levels of BUL present at all sites except Stone River (Site 16), a lowland tributary (Fig. 3, 4, 6). Three patterns of introgression are present in the Burdekin Basin (Fig. 3, 4, 6). Meadowbank Station (Site 16) is immediately adjacent to Herbert Basin and had a similar genetic signature of admixture between NCQ and BUL. Sites in or close to the mainstem Burdekin River above Burdekin Falls were primarily admixed between BUL and ASP, but with some NCQ influence as well. Upper Burdekin River tributary sites and those near to Dalrymple Dam (Site 34, where the Burdekin and Belyando rivers meet) were primarily BUL, but with low levels of introgression from ASP and NCQ. Locations in the Belyando River (Sites 30, 31) were ‘pure’ BUL. The two populations from the lower Burdekin Basin below Burdekin Falls were primarily NCQ, with some influence from CEQ. This influence extended south of the Burdekin Basin, at least to the Don Basin (Site 37). Given the large gap in our sampling between the Don and Pioneer basins, it is unclear how far south this introgression extends. No introgression was evident immediately north of the Burdekin Basin in the Haughton Basin (Site 21).

Our sampling on the northern range limits of NCQ was extremely limited, but the one fish genotyped from Jeannie catchment (Site 1) appeared to be an F₁ hybrid between NCQ and ASP (Fig. 4, 6). Further sampling is necessary to better understand what the patterns might occur in that region.

One perplexing finding was an outlier population of SQN geographically nested within the distribution of SQS (Fig. 3, 4, 6). Tweed River (Site 71) was ‘pure’ SQN, whereas the next population sampled to the north, the Nerang River (Site 70) was well admixed between SQN and SQS (Fig. 6). Further sampling is required to determine the extent of SQN in this region and what the broader patterns of introgression are. It is unclear whether this is a natural pattern or the result of translocations. Unfortunately, we lack mtDNA data for Site 71, but the one fish genotyped from Site 70 was identical to one of the SQN haplotypes (Fig. 7).

It cannot be determined here whether the regions with introgression are stable or changing over longer time frames. Glacial-scale climatic changes result in different degrees of aridity for inland populations, or higher potential connectivity for coastal populations during periods of lowered sea levels (Unmack 2001), which allow populations to mix to a greater extent than today. Although it is clear that these zones of introgression are limited to the close proximity of adjacent taxa or lineages, our data cannot determine what drives the lack of gene flow beyond these narrow contact zones. One possibility is that there are a sufficient number of adaptive genetic differences among these taxa to provide definite selective advantages for each taxon in its own unique climatic region, similar to those proposed for rainbowfishes (Attard et al. 2022; Gates et al. 2023). It is important to note that qualitatively none of the introgressed populations appears to be less abundant than their non-introgressed counterparts.

Biogeography

The Ambassis populations in this study present an eclectic mix of biogeographic patterns, some of which are difficult to explain and are often idiosyncratic relative to other species in the same region.

The most unusual biogeographic pattern is the disjunct population of BUL in the Bulloo Basin. The easiest way for fish to be exchanged between Bulloo and Burdekin basins is via Cooper Creek. Several fish species as well as eastern short-necked turtle (Emydura macquarii) have a close relationship indicative of recent exchanges between Cooper Creek and Burdekin Basin (Unmack 2013; Georges et al. 2018). However, Cooper Creek has a different Ambassis lineage present (taxon ASP), which has presumably replaced BUL within that system, with no evidence of introgression. Given the generally low divergence between BUL populations in the Bulloo and Burdekin basins (e.g. 1 base pair of 1141 for cytb), this replacement is likely to be recent (perhaps in the past few hundred thousand years). Disjunct populations such as this in widespread genera are uncommon in Australian freshwater fishes (and other aquatic organisms). One species, the softspine catfish (Neosilurus mollespiculum), likely shows the same biogeographic pattern. Although its primary distribution is the Burdekin Basin, specimens of this catfish, or one closely related to it, were also recorded from Bulloo River (I.27071-001, Australian Museum Sydney), but have never been found in Cooper Creek.

Most evidence for disjunct patterns at the river basin scale come from genetic studies. It is also important to note that in some cases species may show only a disjunct pattern with their mitochondrial DNA, but not their nuclear DNA. Examples in other groups include patterns based on mtDNA for three Craterocephalus lineages (Unmack and Dowling 2010) within the ‘Disjunct’ group, Lentigenosus clade from north-western Australia, with a sister relationship to Dalhousie clade from central Australia and the ‘Atherton’ clade present in north-eastern Queensland. In addition, there are disjunct populations of Craterocephalus fulvus and Melanotaenia duboulayi in the Fitzroy Basin in eastern Queensland. In rainbowfishes the mitochondrial lineage within M. splendida separates monophyletic western and eastern rainbowfish species such as M. australis and M. eachamensis, resulting in a major disjunction between drainages to the west of Darwin and eastern Australia (McGuigan et al. 2000; Unmack et al. 2013).

Although speculative, we suggest that BUL may have some advantage over the other Ambassis taxa in the more arid parts of the Burdekin Basin above Burdekin Falls, given its broad persistence there, especially in the southern Belyando River, which has no evidence of introgression despite a lack of any clear physical barriers to gene flow. BUL has experienced introgression with NCQ and ASP in parts of the Burdekin River. There is evidence for exchanges between the upper Burdekin and Herbert rivers, given that BUL has introgressed into Herbert Basin populations of NCQ and NCQ has introgressed into BUL in the Burdekin Basin (Fig. 2), with at least one exchange being quite recent. For instance, Wyoming Creek (Site 22) in the upper Burdekin River has clear stronger affinities with other Herbert Basin populations (Fig. 2, 5). This is adjacent to an area with a very low drainage divide, which could allow exchanges between basins during overland flow events. No doubt exchanges are likely to have occurred between these two basins over a considerable period of time, given the general geological stability of the region and a flat topography in some parts of the drainage divide. Also likely recent is an invasion by ASP from western drainages into the mid-Burdekin River, which has spread downstream, being typically found in the main Burdekin channel and the lower reaches of some tributaries (Fig. 2). We have previously found evidence of such an exchange in purple spotted gudgeons (Mogurnda spp.) from either the Flinders or Einsleigh River into Gregory River, a mid-tributary of the Burdekin River (S. Amini, P. J. Unmack and M. Adams, unpubl. data). It is unclear where this exchange is likely to have occurred, although there are numerous swamps and small lakes very close to the drainage divide in this region (e.g. Barkers Swamp, Agnes Lake). It is possible that Ambassis entered the system by the same pathway, a supposition that would need to be confirmed by genotyping Ambassis from Gregory River. Despite these other taxa dispersing into the Burdekin Basin, they have not overwhelmed BUL genotypes, unlike the hypothesised situation in Cooper Creek with ASP replacing BUL. However, it is also possible that ASP have only recently invaded the Burdekin Basin and may replace BUL over time.

There is a lack of any clear physical explanations for the broad geographic boundary between the deep phylogenetic separation of A. agassizii and NCQ, which occurs from the lower Burdekin River (below Burdekin Falls, Sites 35 and 36) south to at least the Don River (Site 37, Fig. 2). Molecular data for two other freshwater fishes show a phylogeographic boundary in the same general region, namely purple spotted gudgeon (Mogurnda spp., nuclear DNA, Adams et al. 2013; S. Amini, P. J. Unmack and M. Adams, unpubl. data) and Pacific blue-eye (Pseudomugil signifer) for mtDNA (Wong et al. 2004). Clearly, some gene flow is ongoing at the boundary, but it cannot be determined whether A. agassizii is expanding north, NCQ is expanding south, or the boundary is stable. Gene flow at the boundary of NCQ and admixed populations may be quite recent, given that introgression is widespread in the lower Burdekin but absent from the Haughton River (site 21) immediately north, which is only minimally separated from the Burdekin, especially during recent lowered sea levels. The northern range limits of NCQ are unclear owing to a lack of samples north of the Russell–Mulgrave Basin (Site 2), with the introgressed population from Muck River (Site 1) being based on a single sample involving a distance of ~330 km separating these two locations. This population occurs immediately east of Normanby Basin, which contains non-introgressed ASP (Fig. 2).

Populations of A. agassizii separated into the following four phylogroups: SQS from the Clarence River (Site 76) north to the Caboolture River (Site 64), SQN from the Mooloolah River (Site 63) north to the lower portion of the Burnett River (Site 55), CEQ from the Kolan River north (Site 48) and MDB from the Murray–Darling Basin and upper portion of the Burnett River. A central melting pot for three of these groups is the Burnett River. Biotic exchanges between the Burnett and the Murray–Darling Basin (Condamine River specifically) are common and include a number of fishes such as carp gudgeons (genus Hypseleotris, Thacker et al. 2007; Unmack et al. 2023), purple spotted gudgeons (genus Mogurnda, Sasaki et al. 2016), flathead gudgeons (genus Philypnodon, Thacker et al. 2008; Hammer et al. 2019b), and hardyheads (genus Craterocephalus, Unmack and Dowling 2010). In most cases, exchanges have been recent (probably within the past 200,000 years) between the Burnett and Murray–Darling Basin (on the basis of very small genetic divergences). It has been hypothesised these exchanges have occurred in both directions across shallow drainage divides (Unmack 2013).

An examination based on mitochondrial sequence data for 33 freshwater species in south-eastern Queensland (from the Mary River south, Page and Hughes 2014) found that a number of species displayed deep divergences between the Brisbane and Mary basins. However, there was considerable variation in where exactly the separation occurred relative to the smaller river basins that separate the Brisbane and Mary. Coastal populations of A. agassizii separated between the Caboolture and Mooloolah rivers, which matched one smaller divergence in Page and Hughes (2014), but other than this divergence it showed low levels of geographic structure across their study area. The partial biogeographic break between the Burnett and Kolan basins (rendered incomplete by the finding a single CEQ outlier in the Burnett at Site 52) also appears to be different from other freshwater fish species examined in this region (e.g. Unmack and Dowling 2010; Sasaki et al. 2016; Hammer et al. 2019b; Unmack et al. 2023).

Conclusions

Our comprehensive molecular data of the A. ‘agrammus’ complex in eastern Australia have demonstrated the likely presence of four candidate species, one referable to A. agassizii, one to A. sp (sensuAllen et al. 2002), one that may represent A. agrammus, and one putative sp. nov. Given that both A. sp and A. agrammus as currently described (Allen and Burgess 1990) are widespread outside of eastern Australia, any detailed assessment of the likely taxonomic affinities of candidate species BUL and ASP will be deferred to a future companion study of the other glassfish populations and species of central, northern and western Australia.

All four candidate species appear to readily hybridise in regions where their distributions abut, although the area occupied by these zones of admixture are a tiny proportion of the total range sampled for ‘pure’ exemplars of the four taxa. Intriguing questions such as the age of these hybrid zones, how far the observed admixture bleeds into parental taxa, whether there is assortative mating, etc. will require considerably more detailed geographic sampling for both morphological and genomic markers. Nevertheless, the systematic, biogeographic and phylogeographic framework outlined here for all candidate species will serve as the ideal scaffold for any such future explorations of this intriguing group.