
==== Front
Syst Biol
Syst Biol
sysbio
Systematic Biology
1063-5157
1076-836X
Oxford University Press US

38577768
10.1093/sysbio/syae016
syae016
Regular Manuscripts
AcademicSubjects/SCI00960
AcademicSubjects/SCI01130
AcademicSubjects/SCI01130
Museum Skins Enable Identification of Introgression Associated with Cytonuclear Discordance
Potter Sally School of Natural Sciences, 14 Eastern Road, Macquarie University, Macquarie Park, NSW 2109, Australia
Division of Ecology and Evolution, Research School of Biology, 134 Linnaeus Way, The Australian National University, Acton, ACT 2601, Australia
Australian Museum Research Institute, Australian Museum, 1 William St, Sydney, NSW 2010, Australia

Moritz Craig Division of Ecology and Evolution, Research School of Biology, 134 Linnaeus Way, The Australian National University, Acton, ACT 2601, Australia

Piggott Maxine P Division of Ecology and Evolution, Research School of Biology, 134 Linnaeus Way, The Australian National University, Acton, ACT 2601, Australia
Research Institute for the Environment and Livelihoods, Charles Darwin University, Casuarina, NT 0811, Australia

Bragg Jason G National Herbarium of New South Wales, The Royal Botanical Gardens and Domain Trust, Mrs Macquaries Road, Sydney, NSW 2000, Australia

Afonso Silva Ana C Univ. Lille, CNRS, UMR 8198 - Evo-Eco-Paleo, F-59000 Lille, France

Bi Ke Museum of Vertebrate Zoology and Department of Integrative Biology, University of California Berkeley, Berkeley, CA 94720, USA

McDonald-Spicer Christiana Division of Ecology and Evolution, Research School of Biology, 134 Linnaeus Way, The Australian National University, Acton, ACT 2601, Australia

Turakulov Rustamzhon Australian Genome Research Facility, Victorian Comprehensive Cancer Centre, 305 Grattan Street, Melbourne, VIC 3000, Australia
Earth Sciences, College of Science and Engineering, Flinders University GPO Box 2100, Adelaide, SA 5001, Australia

Eldridge Mark D B Australian Museum Research Institute, Australian Museum, 1 William St, Sydney, NSW 2010, Australia

Barrow Lisa Associate Editor
Correspondence to be sent to: School of Natural Sciences, 14 Eastern Road, Macquarie University, Macquarie Park, NSW 2019, Australia; E-mail: sally.potter@mq.edu.au.
5 2024
05 4 2024
05 4 2024
73 3 579593
10 3 2022
14 3 2024
03 4 2024
20 4 2024
© The Author(s) 2024. Published by Oxford University Press on behalf of the Society of Systematic Biologists.
2024
https://creativecommons.org/licenses/by-nc/4.0/ This is an Open Access article distributed under the terms of the Creative Commons Attribution-NonCommercial License (https://creativecommons.org/licenses/by-nc/4.0/), which permits non-commercial re-use, distribution, and reproduction in any medium, provided the original work is properly cited. For commercial re-use, please contact reprints@oup.com for reprints and translation rights for reprints. All other permissions can be obtained through our RightsLink service via the Permissions link on the article page on our site—for further information please contact journals.permissions@oup.com.

Abstract

Increased sampling of genomes and populations across closely related species has revealed that levels of genetic exchange during and after speciation are higher than previously thought. One obvious manifestation of such exchange is strong cytonuclear discordance, where the divergence in mitochondrial DNA (mtDNA) differs from that for nuclear genes more (or less) than expected from differences between mtDNA and nuclear DNA (nDNA) in population size and mutation rate. Given genome-scale data sets and coalescent modeling, we can now confidently identify cases of strong discordance and test specifically for historical or recent introgression as the cause. Using population sampling, combining exon capture data from historical museum specimens and recently collected tissues we showcase how genomic tools can resolve complex evolutionary histories in the brachyotis group of rock-wallabies (Petrogale). In particular, applying population and phylogenomic approaches we can assess the role of demographic processes in driving complex evolutionary patterns and assess a role of ancient introgression and hybridization. We find that described species are well supported as monophyletic taxa for nDNA genes, but not for mtDNA, with cytonuclear discordance involving at least 4 operational taxonomic units across 4 species which diverged 183–278 kya. ABC modeling of nDNA gene trees supports introgression during or after speciation for some taxon pairs with cytonuclear discordance. Given substantial differences in body size between the species involved, this evidence for gene flow is surprising. Heterogenous patterns of introgression were identified but do not appear to be associated with chromosome differences between species. These and previous results suggest that dynamic past climates across the monsoonal tropics could have promoted reticulation among related species.

Cytonuclear discordance
exon capture
introgression
rock-wallabies
museum skins
==== Body
pmcComplex evolutionary histories are now being resolved using genomic data sets and advanced computational analyses (e.g., Fontaine et al. 2015; Figueiró et al. 2017; Layton et al. 2020; Rivera et al. 2022). Empirical studies are increasingly revealing evidence of reticulate evolutionary histories across populations and species (see Edwards et al. 2016; Mallet et al. 2016; Aguillon et al. 2022). However, our understanding of the drivers of introgression and its impacts on divergence and speciation across diverse organisms is still limited (e.g., Cruickshank and Hahn 2014; Harrison and Larson 2014; Wolf and Ellegren 2017). By expanding sampling of populations, historical museum specimens provide an important source of genomic information to resolve longstanding questions in systematics, including the role of ancient introgression and hybridization in the speciation process (Card et al. 2021).

Reticulate processes can present as incongruent phylogenetic signals in data (Mallet et al. 2016). Most obviously, strong “cytonuclear discordance” where the divergence in mitochondrial DNA (mtDNA) differs from that for nuclear genes more (or less) than expected given the higher mutation rate and lower population size of mtDNA, is now evident across diverse organisms (see review by Toews and Brelsford 2012; Phuong et al. 2017; Firneno Jr et al. 2020; Sarver et al. 2021). Among the many factors that could contribute to cytonuclear discordance (e.g., Rheindt and Edwards 2011; Bonnet et al. 2017), gene flow among species is important to consider. Our ability to generate population genomic data sets now enables more powerful tests to detect introgression separate from incomplete lineage sorting (e.g., Seixas et al. 2018; Sarver et al. 2021). Generating empirical data across diverse species and ecosystems enables us to understand the role of intrinsic species’ traits (e.g., body size) versus extrinsic, biome-specific dynamics (e.g., environment and climatic stability) in driving these processes (see Singhal et al. 2021). Cytonuclear discordance has been associated with climatic instability and its effects on population range dynamics (e.g., Linnen and Farrell 2007; Good et al. 2008; Cahill et al. 2013; Chavez et al. 2013). Range instability in response to past climate change is expected to increase the opportunity for hybridization and introgression through neutral demographic processes (e.g., Currat et al. 2008; Excoffier et al. 2009; Phuong et al. 2017).

The growing body of molecular research evaluating the phylogenetic and phylogeographic diversity of organisms across the Australian monsoonal tropics (e.g., Eldridge et al. 2011; Pepper and Keogh 2014; Rosauer et al. 2016; Potter et al. 2018) has identified the Kimberley region as a diversification hotspot in vertebrates, plants, and other organisms (Crisp et al. 2001; Powney et al. 2010; Köhler and Criscione 2015; Oliver et al. 2016; Shelley et al. 2018). Along with high phylogenetic endemism in the relatively mesic northwest Kimberley (Rosauer et al. 2018), there is also evidence of interspecific introgression across the whole region for diverse groups (e.g., freshwater fish – Shelley et al. 2020; lizards – Moritz et al. 2016, 2018; Laver et al. 2018; birds – Kearns et al. 2014 and frogs – Jaya et al. 2022). The Kimberley and Top End regions of northern Australia have shown contrasting climatic histories with the Kimberley experiencing more extreme aridity than the Top End over the last 80 Kya (see Reeves et al. 2013; Potter et al. 2018). Thus, it is possible that greater instability of species’ distributions in the Kimberley has increased opportunity for introgression among closely related species.

However, comprehensive sampling of taxa from important but poorly surveyed regions can be problematic. Fortunately, technological advances now enable sequencing of museum specimens which are often our only source of important spatial and/or temporal data (see Card et al. 2021) in remote areas. The ability to unleash this genomic information might be the only evidence to resolve taxonomic questions and complex historical processes of divergence. Expanding geographic coverage within species is important for 2 reasons. First, this will reduce bias in estimates of population statistics that can result from clustered sampling in species with strong isolation-by-distance (Battey et al. 2020). Second, it can avoid the “ghost population” problem, where failure to sample a population involved in past reticulations can mislead analyses of introgression (Beerli 2004; Slatkin 2005; Hey et al. 2018; Linck et al. 2019; Hibbins and Hahn 2022). Although a single individual can be important in systematic research, population sampling provides a foundation for a better understanding of spatial distribution and historical evolutionary processes, particularly where introgression may be restricted in space and time.

Here we incorporate sampling of museum skins and genome-scale sequencing in the brachyotis group of rock-wallabies (Petrogale) to disentangle their complex evolutionary history. Rock-wallabies are a unique and valuable empirical system for understanding evolutionary processes driving divergence and speciation (Potter et al. 2017). They are one of the most speciose living Australo-Papuan marsupial genera with 17 distinct species, 12 subspecies, and 16 distinct karyotypes, and have been used to evaluate the role of chromosome rearrangements in speciation (see Potter et al. 2017, 2022). The genus represents a series of recent and rapid radiations (Potter et al. 2012a) which, previous to the use of genomic methods, has created difficulties in resolving phylogenetic relationships and determining the evolutionary histories of species. Many questions remain in relation to species delimitation, introgression between species, the role of demographic population history, and the interaction of chromosome rearrangements and genic divergence in this genus (see Potter et al. 2017).

The brachyotis chromosomal group, distributed across the Australian monsoonal tropics, represents the oldest radiation in the genus (Potter et al. 2012a). It currently comprises four species, sympatric over parts of their ranges, containing 10 operational taxonomic units (OTUs; identified in Potter et al. 2012b, 2014 on the basis of initial genetic data and morphology): the western short-eared rock-wallaby (Petrogale brachyotis) which includes two subspecies (P. brachyotis brachyotis and P. b. victoriae) that are distinct for mtDNA and morphology and with 2 ESUs defined by mtDNA divergence within the former; Monjon (Petrogale burbidgei) with 2 mtDNA-defined ESUs; Nabarlek (Petrogale concinna) including 3 morphological subspecies (P. c. canescens, P. c. concinna and P. c. monastria); and the eastern short-eared rock-wallaby (Petrogale wilkinsi) including 3 candidate ESUs based on mtDNA divergence and morphology (see Fig. 1). The brachyotis species group has the most rearranged chromosomes amongst Petrogale (2n = 16, 18) compared with the ancestral (2n = 22) macropodid karyotype (Fig. 1; Sharman et al. 1990; Eldridge et al. 1992; Potter et al. 2017). However, there is less chromosomal variation among these species than within the penicillata group, in which genome reorganization is associated with strong reproductive isolation and suppression of gene flow (Potter et al. 2022). There is also substantial body size variation within the brachyotis group, with both P. burbidgei and P. concinna weighing 1–2 kg, the smallest species in the genus, whereas P. brachyotis and P. wilkinsi at 3–5 kg are significantly larger (Potter et al. 2014; Baker and Gynther 2023). Of particular emphasis here are the body size differences which can be expected to influence competitive ability and physiology. In addition, P. concinna is the only marsupial species to have continually erupting molars, a unique morphological feature that previously placed it in the monotypic genus Peradorcas (Thomas 1904). In previous studies of this group, mtDNA was highly discordant with morphological and cytological evidence in that each of the smaller species (P. burbidgei and P. concinna) had polyphyletic mtDNA interspersed amongst the variation within P. brachyotis (Potter et al. 2012b, 2014). Initial analyses of nuclear genes were uninformative (Potter et al. 2012b, 2014), but genomic data (~1800 nuclear genes) for 2 individuals indicated evidence of historical introgression (Potter et al. 2017).

Figure 1. Sampling and karyotypes of the brachyotis group of rock-wallabies from northern Australia. Sample locations are shown for modern (closed symbol) and historical (open symbol) samples for each OTU, which are designated by a unique symbol and the distribution of OTUs are outlined by dashed lines. The location for the unsampled subspecies P. concinna concinna is included for reference in black. Three karyotypes differing by at least one fusion and two centric shifts are known amongst the four species, in addition to some inter and intraspecific variation of the X chromosome. Robertsonian fusions are highlighted by chromosome numbers in reference to the ancestral karyotype, (a) denotes acrocentric chromosome, (sm) denotes submetacentric chromosome, (s) denotes submetacentric chromosome, (m) denotes metacentric chromosome, and polymorphic karyotypes are shown for the X chromosome for P. brachyotis and P. wilkinsi.

In this study we use a targeted exon capture approach, which enables us to incorporate high-coverage genomic data for museum specimens with modern tissue samples, to fill important sampling gaps for the brachyotis group. This approach also allows us to capture mitochondrial sequence data as by-catch from the experiments and assemble more extensive mitochondrial genomic data. Here, we 1) evaluate the relationships and divergence of the previously identified 10 OTUs in this chromosomal group using phylogenomic and population genomic approaches, 2) compare mitochondrial and nuclear divergence histories for discordance, 3) test for introgression between OTUs, and 4) evaluate processes causing cytonuclear discordance between OTUs. The increased genomic coverage, together with more robust sampling, allows us to test models of isolation versus migration between OTUs and statistically assess the divergence histories in this species complex in the context of body size differences and regional climatic instability.

Materials and Methods

Sampling and DNA extraction

DNA was extracted from the ear and liver tissue (modern samples) stored at the Australian Museum, as well as from samples of museum specimens (skull or skin; historical samples). Modern samples were extracted using a salting-out method (Sunnucks and Hales 1996), whereas historical samples were extracted using the DNeasy Blood and Tissue Kit (Qiagen GmbH, Hilden, Germany) in a dedicated environmental and trace DNA laboratory separate from where modern samples were extracted. All DNA extractions were conducted using aerosol barrier pipette tips and all working surfaces and equipment were wiped down with Lookout DNA Erase (Sigma–Aldrich) before each use. A total of 126 samples were analyzed from across the distributions of the 5 species/10 OTUs (excluding the subspecies P. concinna concinna due to lack of sample availability), however due to samples failing to pass assembly and genotype calling filters, the final data set was based on 79 samples (Supplementary Table 1; Fig. 1). This included the following OTUs: P. b. brachyotis East Kimberley—EK ESU (n = 12), P. b. brachyotis West Kimberley—WK ESU (n = 15), P. b. victoriae (n = 6; including one skin replicate), P. burbidgei northern lineage—N (n = 8), P. burbidgei southern lineage—S (n = 3), P. concinna canescens (n = 4), P. c. monastria (n = 5), P. wilkinsi Top End ESU (n = 16), P. wilkinsi Gulf ESU (n = 5), and P. wilkinsi Groote ESU (n = 3). Two P. penicillata individuals were used as an outgroup.

Exon Capture Approach and Bioinformatics

We used a custom in-solution exon capture approach (SeqCap EZ Developer Library; Roche NimbleGen) using target sequences (3960 target exons) designed from a yellow-footed rock-wallaby (Petrogale xanthopus) transcriptome (Bragg et al. 2016; Potter et al. 2017). Orthologs, with targets >200bp in length from BLAST hits to the Tasmanian devil (Sarcophilus harrisii) and tammar wallaby (Notamacropus eugenii) genomes were included. Samples included in the in-solution exon capture experiment had genomic libraries prepared following the protocol of Meyer and Kircher (2010), including modifications made by Bi et al. (2013). Each individual had a unique barcode and samples were pooled in equimolar ratios (1.2μg total) (see Supplementary Table S1). A total of 56 individuals were pooled in any one experiment and samples were prepared across 3 different experiments (SP12, SP13, and SP14). The pooled genomic libraries were combined with 5 μg of mouse Cot-1 DNA (Life Technologies Corporation) and 56 barcode-specific blocking oligos (1000 pmol) designed to block the unique barcodes and adapters used in the Meyer and Kircher (2010) protocol. The hybridization mix was added to the target probes and hybridized for ~72 h following the SeqCap EZ Developer Library protocol. Hybridization reactions were amplified in 2 independent enrichment PCRs (17 cycles) and then cleaned up using the QIAquick PCR purification kit (Qiagen). Quality control checks were made using the DyNAmo Flash SYBR green qPCR kit (Thermo Fisher Scientific Inc.; see Bi et al. 2012) to assess global enrichment of the target exons by comparing precapture pooled genomic libraries to the postcapture cleaned hybridization reaction and specifically designed to hit targets of the hybridization probes. After passing these quality control checks, the enriched hybridization samples were run on a BioAnalyzer (2100; Agilent Technologies, Inc.) to check the quality and quantity of the libraries and then sequenced on a single lane of an Illumina HiSeq 2500 (100 bp paired-end run) at the ACRF Biomolecular Resource Facility.

The raw sequencing data was processed following the workflow of Singhal (2013) which removed duplicate, contaminant, and low complexity reads (scripts for this pipeline are archived in Dryad Repository doi:10.5061/dryad.931zcrjrr). Each locus from modern samples was assembled de novo from the cleaned sequencing reads using the workflow described in Bragg et al. (2016). For historical samples, we used an alternate approach using custom scripts (https://github.com/CGRL-QB3-UCBerkeley/denovoTargetCapturePhylogenomics) and published methods (Bi et al. 2012; Portik et al. 2016). The best-assembled haplotype from each individual (h0) and the diplotype sequences obtained from the previous assembly steps were aligned and filtered using the EAPhy (v1.2; Blom 2015) pipeline. The assembled diplotype alignments were mapped against the Petrogale penicillata genome using ncbi-blast-2.13.0 + (Camacho et al. 2009), which was used to create a BLAST database and BLAST-QC (Torkian et al. 2020) to determine their chromosomal location. Refer to Supplementary Materials for more details.

A by-product of exon capture experiments is the ability to recover some mitochondrial sequence data as by-catch (Guo et al. 2012). We extracted mitochondrial genome data from the raw sequence data following the methods of Hahn et al. (2013) using a docker version for MITObim (https://github.com/chrishah/MITObim) (see Supplementary Materials for details). Data were generated in 1 of 2 approaches: 1) mapped to the common wallaroo (Osphranter robustus) mitochondrial genome (Janke et al. 1997), or 2) reconstructed from mitochondrial ND2, COI, and Cytb seeds from previously sequenced individuals (Potter et al. 2012b, 2014). Final mitochondrial genome assemblies were then aligned using Geneious Prime as well as mafft (Katoh et al. 2002) and edited to protein-coding regions of the genome (~11,500 bp).

Phylogeographic Structure and Divergence

We estimated the divergence histories of individuals using the concatenated nuclear haplotype (h0) dataset. Given the low phylogenetic signal within each locus, coalescent species tree approaches such as ASTRAL-III that form species trees from gene trees may be inconsistent (Wascher and Kubatko 2021). Therefore, we used SVDquartets (Chifman and Kubatko 2014) executed in PAUP* (v 4.0a; Swofford 2003) to estimate genetic relationships of nuclear loci, where sites sampled from the concatenated sequence data are analyzed under the Multi Species Coalescent (MSC) model. SVDquartets computes quartet scores from a decomposition matrix of site pattern frequencies to infer a phylogeny. First, we evaluated 100,000 random quartets and 1000 bootstrap replicates to infer the relationships of individuals. We then estimated divergence times of OTUs using qage in SVDquartets by assigning individuals to their OTU and implementing the MSC approach using 100,000 random quartets and 1000 bootstrap replicates. For the divergence dating, we used a standard bootstrapping approach, the exact JC model, an outwidth of 140, and rooted the tree with P. penicillata. The estimated time in coalescent units from the analysis was converted to years by the following equation: divergence time = coalescent units * generation time * 2N, where N = theta/(4 * u). The qage analysis calculates an average theta for the run and we used the mutation rate (u) as 1.45 × 10−8 (based on those used in another marsupial, the koala, Phascolarctos cinereus; Johnson et al. 2018). We note this mutation rate is likely too fast, given our target data being nuclear exons. We also estimated the relationships of individuals using a phylogenetic network approach (NeighborNet) in SplitsTree4 (Bryant and Moulton 2002; Huson and Bryant 2006). The approach clusters individuals based on an uncorrected P distance matrix to examine the patterns of reticulation. Although this does not provide bootstrap support for phylogenetic relationships, it provides a visual assessment of complex evolutionary processes of reticulation.

For mtDNA, we applied the maximum likelihood approach in RAxML (v8; Stamatakis 2014), using the rapid bootstrap analysis with 100 bootstrap replicates and the GTRGAMMA model to estimate the relationships of individuals for both concatenated mitochondrial data sets (all assembled mitogenome data and a dataset where individuals had at least 50% data). Osphranter robustus was used to root the tree. To assess the support of the mitochondrial data for the nuclear topology, we compared the likelihood of the mitochondrial phylogeny (unconstrained) to a phylogeny constrained to the nuclear topology in IQtree2 (Minh et al. 2020). We compared likelihoods using the RELL approximation and the approximately unbiased (AU) test with 10,000 replicates and the GTR + G model.

Lastly, we estimated average genetic differentiation (Dxy) between OTUs and species for both mtDNA and nDNA, as well as in relation to chromosome rearrangements. Genetic differentiation was first estimated for the concatenated nuclear h0 dataset (including missing data) in PopGenome in R (Pfeifer et al. 2014). We then calculated Dxy for individual exons and compared average Dxy for exons on rearranged (R) chromosomes to those from nonrearranged (NR) chromosomes for each pairwise comparison between OTUs with varied karyotypes (e.g., all interspecific comparisons excluding those between P. wilkinsi and P. brachyotis OTUs). We averaged across loci specific to differences in karyotype. So, for P. brachyotis and P. wilkinsi OTU comparisons with P. burbidgei OTUs (R: chromosomes 5,8,9; NR: chromosomes 1,2,3,4,6,7,10); P. brachyotis and P. wilkinsi OTU comparisons with P. concinna OTUs (R: 4,5,8,9; NR:1,2,3,6,7,10); P. burbidgei OTUs comparisons with P. concinna OTUs (R: 4; NR:1,2,3,5,6,7,8,9,10). We excluded the X chromosome from this analysis as there were only four loci mapped.

Historical Demographic Analyses

As introgression events can be associated with range fluctuations, we tested for population expansion. We estimated Tajima’s D on individual OTUs for the concatenated h0 alignment in the PopGenome R package (Pfeifer et al. 2014). Coalescent simulations were run for 1000 simulations using Hudson’s MS (Hudson 2002) to evaluate if the Tajima’s D value was significant. Significance was calculated by comparing the observed Tajima’s D to the P = 0.05 threshold from the simulated data. We also tested for demographic expansion using the rangeExpansion R package which applies the approach outlined in Peter and Slatkin (2013, 2015). This approach estimates the directionality index (ψ) which assesses asymmetries in two-dimensional allele frequency spectrum between populations, assuming that asymmetries are caused by founder events that occur due to population expansion. Low-frequency alleles are lost from founder events creating clines in the allele frequency SFS (Peter & Slatkin 2013). Each colonization event from an expansion origin is expected to be accompanied by a founder event and shifts in allele frequencies following successive colonization events. In an equilibrium isolation-by-distance model ψ = 0, but under models of population expansion ψ > 0. Evaluating the SFS is more sensitive to patterns of range expansion than Tajima’s D which uses pairwise differences and segregating sites. This analysis was only performed on a subset of OTUs with appropriate sample sizes (> 5 individuals), including P. wilkinsi TE, and the 2 P. b. brachyotis OTUs. This approach requires the ancestral state of alleles to be known and uses SNP data to generate a site frequency spectrum. We gathered SNP data sets for this analysis from the exon sequence alignments using the EAPhy pipeline (Blom 2015), including P. penicillata as the outgroup to call derived SNPs. For each OTU, SNPs (a single SNP per locus) were output into a single file and SNPs heterozygous for the outgroup were removed.

Coalescent Models of Divergence Histories

We applied 2 approaches that exploit having extensive sequence data to either test for the presence of introgression, or to estimate the amount of introgression, in each case separate from incomplete lineage sorting. First, we used an Approximate Bayesian Computation framework to statistically test for isolation with migration using the demographic inference with linked selection program (DILS; Fraïsse et al. 2021). Using 2-population models, we focused on testing divergence among species where there is evidence of substantial cytonuclear discordance (i.e., between interspecific OTUs that are sympatric), but also included comparisons that were most informative for understanding divergence between intraspecific OTUs (23 comparisons total). Using the h0 haplotype dataset we estimated (1) the best demographic model of divergence (e.g., isolation—strict isolation or ancient migration versus migration via secondary contact versus isolation with migration, Roux et al. 2016) and (2) the best model of genomic divergence, allowing for population size and migration rates that are either homogeneous or heterogeneous across loci and so incorporating potential effects of linked selection. Analyses were run with the same demographic priors for all pairwise comparisons, including filtering to a maximum 10% of missing data in loci with a minimum length of 30bp for which at least 2 samples had sequence data. After optimization of the priors, the models were performed with the following priors: changes in population size, an assumed mutation rate of 3 × 10−9 mutations per generation and per base pair (which we note is slower than the previous rate used for dating but is optimized for the best model fit), 0.1 of ratio of recombination over mutation, a time of split varying between 100 and 10 Ma, population sizes between 100 and 1.2 million, and finally modeling the barrier loci with the Bimodel distribution, although taking migration rates as a minimum of 0.4 and maximum of 20. Five replicate runs per analysis were compared for convergence and the run with the highest posterior probability was used. We also used the computed Dxy (and Da) from DILS and compared these with the probability of migration from the best demographic model. Here we discuss only the model estimation rather than the parameter estimates from DILS, due to limitations in parameter estimation described in Fraïsse et al. (2021).

For the second analysis we used a Bayesian coalescent-based approach to estimate mutation-scaled population sizes (Θ) and migration rates averaged over the whole divergence history (M) in MIGRATE 4.0 (Beerli and Palczewski 2010; Beerli et al. 2019). We estimated Θ and M for 3 separate comparative analyses as motivated by comparisons of mtDNA and nDNA divergences. The first examined the divergence histories between the 3 P. wilkinsi OTUs, the second examined the divergence histories between all OTUs of P. brachyotis, P. burbidgei, and P. concinna, and the third examined divergence histories between parapatric P. wilkinsi Top End and P. b. victoriae OTUs. These analyses were run using the diplotype nuclear data set as this program can evaluate ambiguity-coded sequence data, the Jukes-Cantor sequence model with base frequency 0.25, a uniform prior distribution for Θ (0,0.006–0.0005, 0.0006–0.00005; minimum, maximum, delta, respectively) and M (0, 100,000, 10,000). Random starting parameters from the prior distribution were used for estimation of Θ and M, and analyses were run starting from a random number seed, a constant mutation rate estimated from the data. Two independent analyses were run for each comparison, each with one long chain and four heated chains. A static heating scheme was used with temperatures set to 1, 1.5, 3, and 106 ordered from cold to hot, with sampling every 100 steps and run for one million steps after a burn-in of 100,000.

Results

Genomic Data sets

We generated 4 different data sets from the exon capture experiments: 1) a phased haplotype dataset (h0) from the nuclear loci (396,135 bp; 841 exons); 2); a diplotype dataset from the nuclear loci (589,797 bp; 1094 exons); 3) a mitochondrial dataset encompassing protein-coding portions of the mitochondrial genome (~11,500 bp); and 4) a reduced mitochondrial dataset where individuals had at least 50% data. A total of 26 skins (46%) were excluded from our experiment, largely due to either complete failure of recovering alignments from a sample, or insufficient data for robust analyses. The final nuclear data sets contained a complete data matrix of individuals for each locus with 1.6% missing data overall per individual. The success in recovering mitochondrial sequence data as by-catch was extremely variable, with missing data ranging from 0% to 81% (mean 35%; Supplementary Table S2). Allowing for up to 50% missing data, we were able to recover mitochondrial sequence data for 40 individuals which represent the geographic spread of all 10 OTUs.

Phylogeographic Structure and Divergence

The coalescent species analysis and phylogenetic network based on the nuclear haplotype data resolved species as monophyletic groups (Fig. 2). Not all the designated OTUs were recovered as monophyletic, with the subspecies of P. brachyotis and P. concinna being paraphyletic (Supplementary Figs. S1 and S2). The 2 OTUs of P. b. brachyotis formed separate monophyletic lineages but the P. b. brachyotis EK clade included P. b. victoriae. There was strong support for two monophyletic OTUs in P. burbidgei (N and S lineages, separated by the Prince Regent River). With addition of skin samples to fill key sampling gaps (Fig. 1), our results support 3 divergent ESUs within P. wilkinsi (P. wilkinsi from the Top End, P. wilkinsi GULF from Sir Edward Pellew Islands and the southern Gulf of Carpentaria region of the mainland, and P. wilkinsi GROOTE from Groote Eylandt).

Figure 2. Comparison of mtDNA and nDNA phylogenetic relationships within sampled brachyotis group taxa. (a) Mitochondrial phylogenetic tree based on concatenated alignments with <50% missing data. Sample symbols from (2a) match OTUs outlined in (2b) and Figure 1. (b) Nuclear phylogenetic tree of the brachyotis group based on the multispecies coalescent analysis in SVDquartets. All bootstrap values (in 2a and 2b) are 100 unless indicated otherwise for major nodes in the tree.

The estimated divergence times between intraspecific OTUs ranged from 29 thousand years ago (kya) (between P. brachyotis OTUs) to 136 kya (between P. concinna OTUs), whereas interspecific divergence estimates ranged from 183 kya to 278 kya. Divergence time estimates between P. b. brachyotis EK and P. brachyotis victoriae were effectively zero. Interspecific divergence between P. brachyotis and P. wilkinsi was 183 kya, and between P. burbidgei and P. concinna was 256 kya, while the deepest split in the clade is 278 kya (Supplementary Table S3). Average Dxy for nDNA between intraspecific comparisons ranged from 0.0008 between P. b. victoriae and P. b. brachyotis EK to 0.0017 between the P. concinna OTUs. For interspecific comparisons average Dxy ranged from 0.0021 to 0.0027 and showed similar trends to divergence times. The mtDNA showed different patterns with average Dxy for interspecific comparisons lower than intraspecific comparisons (see Supplementary Table S3).

In contrast to the nuclear phylogenies, the mitochondrial data sets supported an alternate relationship of OTUs. All species except P. wilkinsi were nonmonophyletic in the mtDNA phylogeny (Fig. 2; Supplementary Fig. S3). Of major discord were the placements of P. burbidgei S, P. b. victoriae and the 2 P. concinna subspecies. Consistent with the nuclear data, the samples for P. wilkinsi GULF and GROOTE OTUs were each monophyletic and distinct from the P. wilkinsi Top End OTU (Supplementary Fig. S3), but several samples have >50% missing data and so are not in Fig. 2. P. wilkinsi forms the sister group to P. brachyotis, P. burbidgei and P. concinna, but one sample groups with P. b. victoriae, although the data set consisting of more missing data (Supplementary Fig. S3). Given the patchiness and large amount of missing data (see Supplementary Table S2), particularly for the first dataset, we are cautious not to infer too much based on these results. However, P. c. monastria was the only lineage to show a discrepancy between trees from the 2 mtDNA data sets. It did not form a monophyletic lineage based on the dataset with more missing data which could be an artefact, as only a single individual remained when allowing for <50% missing data. We therefore draw more heavily from the more complete data with fewer individuals (Fig. 2).

When testing the constrained nuclear topology with the mitochondrial dataset, we found the constrained (i.e., nuclear) topology was rejected, P-AU < 0.05 and a difference in the log-likelihood of 754.34 between the trees. The other associated tests also supported the unconstrained tree (Supplementary Table S4) indicating that the nuclear topology is not supported by the mitochondrial data.

Tests for Introgression Where There is Cytonuclear Discordance

We tested for introgression between OTUs to determine if cytonuclear discordance was caused by introgression. We find higher migration probabilities are largely linked to lower Dxy (e.g., intraspecific comparisons) compared with those between species (Fig. 3). The ‘outliers’ here, with high migration probability and high Dxy, are between P. b. victoriae and P. concinna subspecies, P. b. victoriae and the P. wilkinsi Top End OTU, and between P. b. brachyotis WK and P. burbidgei N. The MIGRATE results between OTUs from distinct species indicated low levels of introgression, although greatest from P. burbidgei N into P. c. monastria (0.001–0.19; Supplementary Table S5). There were no significant differences in average Dxy between interspecific OTUs for rearranged versus nonrearranged chromosomes (Supplementary Fig. S4).

Figure 3 Probability of migration results from DILS analysis for both intraspecific and interspecific lineage pairwise comparisons (where there is evidence of cytonuclear discordance) on the y-axis plotted against the average Dxy between lineages on the x-axis. These results include whether the loci were homogeneous (homo) or heterogeneous (hetero) for migration (M) and population size (N) estimates. Open symbols correspond to intraspecific (within species—W) comparisons and closed symbols correspond to interspecific (between species—B) comparisons. Unique symbols correspond to the best model result from the DILS analysis: square (IM) = isolation-with-migration model with heterogeneous M and N; circle (IM) = isolation-with-migration model with homogeneous M but heterogeneous N; hexagon (AM) = ancient migration model with heterogeneous N; cross (AM) = ancient migration model with homogeneous N; triangle (SC) = secondary contact model with heterogeneous M but homogeneous N; star (SC) = secondary contact model with homogeneous M but heterogeneous N. OTUs include: bbE (P. b. brachyotis EK), bbW (P. b. brachyotis WK), bv (P. b. victoriae), bN (P. burbidgei N), bS (P. burbidgei S), cc (P. c. canescens), cm (P. c. monastria), w (P. wilkinsi Top End), wGR (P. wilkinsi GROOTE), and wGU (P. wilkinsi GULF).

P. brachyotis and P. burbidgei: P. brachyotis, and P. burbidgei were estimated to have diverged ~278 kya (Supplementary Table S3). There was evidence of recent introgression from secondary contact between P. b. brachyotis WK and P. burbidgei N (posterior probability—Pp = 0.7; Supplementary Fig. S5). All other comparisons supported models of ancient introgression (pp = 0.5–0.7). The MIGRATE results revealed low levels of introgression, but with some asymmetric gene flow from P. brachyotis subspecies into P. burbidgei S (0.2–0.3 vs. 0.001) and higher values than those for P. b. brachyotis WK and P. burbidgei N (Supplementary Table S5).

P. brachyotis and P. concinna: P. brachyotis and P. concinna were estimated to have diverged ~278 kya (Supplementary Table S3). Isolation-with-migration models were supported for analyses of P. b. victoriae with both P. concinna OTUs (pp = 0.6–0.9), compared with models of ancient migration for OTUs of P. b. brachyotis and P. concinna (pp = 0.5–0.9; see Supplementary Fig. S4). Migration estimates from MIGRATE were low (~0.1 migrants per generation; Supplementary Table S5).

P. b. victoriae and P. b. brachyotis/P. wilkinsi: Phylogenetic estimates of divergence time between P. brachyotis and P. wilkinsi (~183 kya) was 6-fold greater than between the P. b. brachyotis WK and other OTUs in this species (Supplementary Table S3). Between P. b. victoriae and P. b. brachyotis EK the divergence time estimates from species-trees were effectively zero, in contrast to the divergence with secondary contact model from DILS. Analyses of divergence histories between P. brachyotis subspecies and between P. b. victoriae and P. wilkinsi Top End OTU revealed support for introgression. Between P. b. victoriae and each of the P. b. brachyotis OTUs, a model of gene flow upon secondary contact best supported the data (Fig. 3; pp ~0.65, Supplementary Fig. S5). In contrast, a model of isolation-with-migration was supported for P. b. victoriae and P. wilkinsi Top End OTU (pp = 0.6; Supplementary Fig. S5). MIGRATE results indicated higher gene flow from both P. b. brachyotis OTUs into P. b. victoriae, despite still being low in general (0.22–0.31; Supplementary Table S5) and similar to levels detected between P. wilkinsi Top End OTU and P. b. victoriae (0.21 and 0.32).

Intraspecific Divergence and Tests for Recent Expansion

Generally, the DILS analyses between intraspecific OTUs supported models of migration (isolation-with-migration and secondary contact) (Fig. 3). The clear “outlier” here was between the island P. wilkinsi GROOTE OTU and the 2 other P. wilkinsi OTUs, where there was a low probability of migration despite low Dxy (Fig. 3). A model of ancient migration was supported for these comparisons (pp ~0.55–0.85, Supplementary Fig. S5). Estimated divergence times between these OTUs range from ~95 to 109 kya (Supplementary Table S3).

MIGRATE analyses estimated very low population sizes for P. wilkinsi OTUs (Θ = 4 Neµ = 2 × 10−6 to 8 × 10−5; see Supplementary Table S5 for Θ estimates) and very limited gene flow between OTUs (0.0003–0.04 genomes per generation). The population sizes for P. brachyotis, P. burbidgei and P. concinna OTUs were higher than those of P. wilkinsi (~0.0001) but the gene flow between intraspecific OTUs was variable (Supplementary Table S5). Most estimates were low (0.001–0.1), however some indicated higher levels of gene flow (0.21–1.36 genomes per generation) and included P. b. brachyotis OTUs and P. b. victoriae, and P. burbidgei OTUs.

None of the OTUs indicated deviations from mutation-drift equilibrium in the form of significantly negative Tajima’s D but the largest deviations were observed in P. brachyotis and P. wilkinsi (Supplementary Table S6). The ψ statistic was estimated only for populations with large sample size (P. brachyotis OTUs and P. wilkinsi Top End OTU). Among these, significant demographic expansion was only detected for P. b. brachyotis EK OTU (P < 0.001).

Discussion

Historical Museum Specimens Resolve Complex Evolutionary Patterns

The increase in power from genomic data sets to differentiate between patterns of recent versus historical introgression is increasing our knowledge of the processes driving divergence and speciation (e.g., Sarver et al. 2021; Guo et al. 2022; Hibbins and Hahn 2022). Cytonuclear discordance is often detected, and our ability to now disentangle ILS from introgression is enabling unique insights into the speciation process, exposing how permeable species boundaries can be (Harrison and Larson 2014; Payseur and Rieseberg 2016). Here, with extensive geographic and genomic coverage through incorporation of museum skins to represent remote locations across the brachyotis group of rock-wallabies we find evidence of both historical (between OTUs which diverged 183–278 kya) and recent introgression. Incorporation of historical specimens was essential in this study to apply robust statistical approaches to disentangle their complex evolutionary history. Expanding geographic sampling enabled evaluation of the role of climatic instability in driving processes of cytonuclear discordance. Our results are consistent with growing evidence for links between geographic proximity of species, climatic fluctuations and range stability with introgression (see Singhal et al. 2021).

The more extensive nuclear data set available from museum specimens has enabled population sampling to resolve the species phylogeny and the complex evolutionary processes driving discordant patterns of evolution from the mitochondria and nuclear data sets. For the first time, we find strong support for the monophyly of each of the brachyotis group species, grouping species of similar body size and karyotype (Fig. 2). P. burbidgei and P. concinna which share a 5–9 fusion and a metacentric chromosome 8 formed a sister relationship (diverged 256 kya), whereas P. brachyotis and P. wilkinsi which share the alternative karyotype, larger morphology and ear size, formed a separate sister pair (diverged 183 kya). We note our molecular divergence estimates are orders of magnitude lower than previous fossil-calibrated dates (Potter et al. 2012a). The mutation rate is likely an overestimate because our data is focused on exons which are likely under purifying selection and therefore divergence estimates here would be a minimum. However, estimates do still imply Pleistocene divergence times between species.

Introgression, A Driver of Cytonuclear Discordance

This study substantiates previous suggestions of cytonuclear discordance, which was detected for OTUs of P. brachyotis, P. burbidgei and P. concinna (Potter et al. 2012b, 2014; refer to Fig. 2; Supplementary Figs. 1–3). We find evidence of historical and recent introgression (secondary contact and speciation with gene flow) between these species despite their differences in body size. There is no evidence that the few chromosome rearrangements are suppressing gene flow, which contrasts with a strong role of chromosome change in the related penicillata group (Potter et al. 2022). Our results add to the growing number of studies finding introgression rather than ILS driving cytonuclear discordance (e.g., Lin et al. 2019; Taylor et al. 2021; Wang et al. 2022; Ji et al. 2023). The incorporation of historical samples gave us the population sampling to apply ABC modeling to distinguish between models of ongoing migration, secondary contact, and ancient migration. It also enabled population demographic analyses to understand at a lineage level the extent of gene flow in relation to past demographic histories. A single individual and phylogenetic analysis would not have enabled us to tease apart intrinsic and extrinsic factors in driving cytonuclear discordance.

Interestingly, the cytonuclear discordance we see in this study is largely concentrated in the Kimberley and Victoria River regions of the monsoonal tropics, areas suggested to be less climatically suitable for much of the last glacial cycle, rather than the Top End which has been variable but consistently mesic (e.g., Reeves et al. 2013; Potter et al. 2018). Heterogeneity in climatic stability in the Kimberley has been linked to populations contracting and expanding (e.g., Potter et al. 2016, 2018; Afonso Silva et al. 2017; Fenker et al. 2021; Jaya et al. 2022) and a foundation for gene flow to occur between populations as they reconnect (e.g., Catullo and Keogh 2014; Eldridge et al. 2014; Kearns et al. 2014; Potter et al. 2018). Range instability has been associated with cytonuclear discordance across many diverse species (e.g., Krosby and Rohwer 2009; Singhal and Moritz 2012; Phuong et al. 2017) and has been associated more broadly with macroevolutionary patterns of introgression (see Singhal et al. 2021). We find evidence of spatial expansion for P. b. brachyotis EK in the Kimberley which we predict may contribute to some of the cytonuclear discordance detected.

ABC modeling of interspecific divergence histories support a genic view of speciation (Wu 2001). Models of ancient migration indicate initial divergence could have been accompanied by gene flow and support previous results using phylogenetic network analysis based on a small number of individuals for the brachyotis group (Potter et al. 2017). Many studies have found similar processes of speciation with heterogeneous gene flow (e.g., Roux et al. 2016; Peñalba et al. 2019). We do note however, although DILS has been reported to accurately discriminate between models of isolation and ongoing migration, it is poorer at discriminating between strict isolation and ancient migration, with a higher tendency to support the latter (Fraïsse et al. 2021).

Despite many interspecific divergence models supporting ancient migration, we also find instances of gene flow transcending species boundaries and substantial size differences (1–2 kg vs. 5 kg). Evidence of secondary contact and recent gene flow were detected between P. burbidgei N and P. b. brachyotis WK, and P. concinna OTUs and P. b. victoriae and likely contributed to patterns of cytonuclear discordance. P. burbidgei and P. b. brachyotis WK diverged ~278 kya but are broadly sympatric in parts of the Kimberley region. Divergence models supported recent introgression from secondary contact between P. b. brachyotis WK and P. burbidgei N which could well have resulted in the discordant relationships of P. burbidgei N mtDNA (Figs. 2 and 3). No introgression was detected involving P. burbidgei S, although low levels of gene flow were detected from P. b. brachyotis WK into P. burbidgei S from MIGRATE analysis (Supplementary Table S4).

Recent introgression was also detected between P. b. victoriae and the two P. concinna OTUs which diverged ~278 kya (Supplementary Fig. S5). Given the polyphyly of P. concinna with P. b. brachyotis in the mtDNA phylogeny, not P. b. victoriae (Fig. 2), it is difficult to disentangle how this migration could be a driver for this pattern of incongruence. Gene flow between P. b. brachyotis EK and P. b. victoriae may be involved. Analysis of the currently unsampled P. c. concinna subspecies which is sympatric with P. b. victoriae may shed light on this, however only a single specimen of this taxon exists, collected in 1839 (Eldridge 1997). It is also unresolved if P. c. monastria mtDNA forms a monophyletic lineage based on current low support and missing data in the mtDNA data set (Supplementary Fig. S3). In general, even including skins, we have a very limited sampling of P. concinna relative to its recorded range.

Evidence of recent introgression was also detected between P. b. victoriae and P. wilkinsi Top End OTU which diverged ~183 kya. (Fig. 3). Unlike the nDNA phylogeny, in the mtDNA phylogeny P. b. victoriae was basal to the P. brachyotis/P. burbidgei/P. concinna clade. Introgression with P. wilkinsi could result in the shift of the P. b. victoriae lineage to the base of the mtDNA tree as it mixed with the more divergent P. wilkinsi. Given that these two interspecific OTUs are parapatrically distributed, opportunity for introgression is possible and their past distributions may even have overlapped. The amount of gene flow between species is often linked to molecular divergence, however, speciation occurs across a continuum and many organisms do not fit this trend (see Roux et al. 2016).

At an intraspecific level, there was evidence of cytonuclear discordance between P. brachyotis OTUs (Fig. 2). P. b. victoriae has a highly divergent and unique mitochondrial lineage, as well as pelage differences (Potter et al. 2014), yet clustered with P. b. brachyotis EK in the nuclear results. The SVDquartet analysis could not estimate a divergence time between P. b. brachyotis EK and P. b. victoriae, however, the ABC modeling of nuclear genes, supported a history of divergence with recent introgression upon secondary contact between these subspecies (Fig. 3). This model would appear more parsimonious with the mitochondrial results. Demographic expansion was detected for P. b. brachyotis EK and a scenario of allele surfing of the nuclear genome from P. b. brachyotis EK into P. b. victoriae could explain the cytonuclear discordance and the divergent mitochondrial lineage for P. b. victoriae. Demographic expansion has been associated with cytonuclear discordance in other systems (e.g., Wilson and Bernatchez 1998; Melo-Ferreira et al. 2005; Cahill et al. 2013; Phuong et al. 2017). Allele surfing can target standing variation and more commonly occurs in small populations and in rapidly growing populations (e.g., recent expansion) (see Excoffier et al. 2009). Further data are required to test the demographic expansion hypothesis as a driver of introgression, as signals of spatial expansion can result from other factors (e.g., gene flow between divergent lineages and allele surfing; Marchi and Excoffier 2020) which have been detected in this system. Whole genome sequencing, or a clinal study at the contact zone between subspecies is required to understand the exact drivers of divergence. This would enable assessment and processes of gene flow across the genome, as adaptation or divergence may only affect a small proportion of the genome.

Unlike the divergence and secondary contact detected between P. b. brachyotis EK and P. b. victoriae, we also find evidence of divergence at an intraspecific level within P. wilkinsi (Fig. 3). Previous studies had limited power to assess how distinct OTUs from GROOTE and GULF were, with only a single representative sample for each (Potter et al. 2012b, 2014). Here, with the inclusion of museum samples and increased genomic coverage, our results support the three distinct ESUs within P. wilkinsi which may represent additional taxa (diverged 95–109 kya). Only the 2 mainland P. wilkinsi OTUs (Top End and GULF) show evidence of gene flow. The divergence of the 3 P. wilkinsi OTUs is likely driven by strong biogeographic barriers involving a break in rocky habitat east of the Roper River and the Arafura Sea isolating Groote Eylandt from the mainland. This pattern of divergence has also been seen for other organisms, including lizards (Smith et al. 2011; Rosauer et al. 2016; Oliver et al. 2020), rodents (Kitchener 1989), and birds (Ford 1978).

Heterogeneity in Genome Divergence and Evolution

We detected heterogeneous patterns of migration at both an interspecific and intraspecific level which could suggest the presence of barrier loci (Fig. 3; Supplementary Fig. 5). There is growing evidence that introgression does not occur evenly across a genome (see Aguillon et al. 2022). Both recombination landscapes and density of genes, as well as selection, are likely linked to which introgressed regions are retained or lost (e.g., Martin et al. 2019; see Harrison and Larson 2016). Interestingly, the average Dxy between OTUs showed no significant difference between exons on rearranged versus nonrearranged chromosomes (Supplementary Fig. 4). The exons targeted in this study are dispersed sparsely across the chromosomes and may not be suitable to detect fine-scale divergences associated with chromosome rearrangements (e.g., recombination suppression model of divergence; Rieseberg 2001). Chromosome level, whole genome data is required to determine the landscape of genomic divergence between species and OTUs in the brachyotis group, particularly in relation to structural variation and introgression. However, the modest chromosomal differences amongst these 4 brachyotis group species (e.g., single fusions) may not impede gene flow (see Rieseberg 2001) unless they are associated with adaptive loci under strong selection (see Guerrero and Kirkpatrick 2014). Heterogeneous migration, as seen here with the exon data, may just reflect variation in their strength of selection against introgression. Positive selection driving mitochondrial introgression has been supported across diverse systems (e.g., Melo-Ferreira et al. 2014; Morales et al. 2018) and could be an important mechanism driving cytonuclear discordance in this system (see Bonnet et al. 2017).

Disentangling Factors Ivolved in Introgressive Hybridisation

Our analysis of the brachyotis group of rock-wallabies adds to a growing number of genomic studies that find evidence of historical and recent introgression (e.g., Park and Park 2020; Ferreira et al. 2021) contributing to cytonuclear discordance and complex evolutionary histories. Despite evidence of introgression, the nuclear phylogeny at a species level is well resolved, suggesting that gene flow might be localized across the genome and potentially driven by strong selective pressures (either mitochondrial or chromosomal). Both intrinsic and extrinsic factors can predispose taxa to introgressive hybridization. Intrinsic factors involve changes to the genome, either through genetic drift or selection, including genic or structural variation (e.g., chromosome rearrangements), or phenotypic changes (e.g., life-history effects). In the brachyotis system, intrinsic barriers did not appear to impede historical introgression despite stark differences in body size and minor chromosome differences. Instead, our results highlight the over-riding role that extrinsic climatic instability and range fluctuations can have in driving introgressive hybridization. Hybridization can occur differently in space and time (Harrison and Larson 2014), and it is only with extensive geographic sampling, here enabled by using museum specimens, that we can disentangle historical introgression from incomplete lineage sorting.

supplementary material

Data available from the Dryad Digital Repository: http://doi.org/10.5061/dryad.931zcrjrr.

Acknowledgements

We thank the numerous people, groups and organisations who provided samples or assisted with sample collection. We especially thank Miriuwung Gajerrong rangers and Miriuwung Gajerrong Corporation, as well as staff from the WA Department of Biodiversity, Conservation and Attractions and Australian Wildlife Conservancy. Dunkeld Pastoral Co. Pty. Ltd. and the Water Authority of Western Australia are thanked for access to their properties, financial and sampling support. We are grateful to Leo Joseph of the Australian National Wildlife Collection (https://ror.org/059mabc80), for assistance in undertaking this research. We also thank the staff at the Australian Museum (Sandy Ingleby), Western Australian Museum (Mark Harvey), and Museum Victoria (Kevin Rowe) for permission to sample specimens in their care, as well as the South Australian Museum Australian Biological Tissue Collection and Australian Museum for access to tissue samples and data. We also thank Matt Morgan and Leo Joseph for helpful discussions about the project.

Conflict of Interest

The authors declare no conflict of interest.

Funding

This research was funded by the Australian Research Council Discovery Project Grant (DP160100187), a Centre for Biodiversity Analysis Ignition Grant and the Australian Museum Research Institute. S. Potter is a recipient of an Australian Research Council Future Fellowship (FT210100715). M. Piggott was a recipient of a Discovery Early Career Research Award from the Australian Research Council (DE130100777).

Data Availability

The scripts and data underlying this article are available in Dryad Digital Repository doi:10.5061/dryad.931zcrjrr. In addition, sequence data is available at the NCBI Sequence Read Archive (SRA) under Bioproject PRJNA# (TBA; refer to Supplementary Table 1 for Biosample IDs—TBA).
==== Refs
References

Afonso Silva A.C. , BraggJ.G., PotterS., FernandesC., Manuela CoelhoM., MoritzC. 2017. Tropical specialist vs. climate generalist: diversification and demographic history of sister species of Carlia skinks from northwestern Australia. Mol. Ecol. 26 :4045–4058.28543871
Aguillon S.M. , DodgeT.O., PreisingG.A., SchumerM. 2022. Introgression. Curr. Biol. 32 :R865–R868.35998591
Baker A.M. , GyntherI.C., editor 2023. The mammals of Australia. 4th ed. Sydney: New Holland. p. 317–340.
Battey C.J. , RalphP.L., KernA.D. 2020. Space is the place: effects of continuous spatial structure on analysis of population genetic data. Genetics. 215 :193–214.32209569
Beerli P. 2004. Effect of unsampled populations on the estimation of population sizes and migration rates between sampled populations. Mol. Ecol. 13 :827–836.15012758
Beerli P. , MashayekhiS., SadeghiM., KhodaeiM., ShawK. 2019. Population genetic inference with MIGRATE. Curr Protoc Bioinformatics. 68 :e87.31756024
Beerli P. , PalczewskiM. 2010. Unified framework to evaluate panmixia and migration direction among multiple sampling locations. Genetics. 185 :313–326.20176979
Bi K. , LinderothT., VanderpoolD., GoodJ.M., NielsenR., MoritzC. 2013. Unlocking the vault: next-generation museum population genomics. Mol. Ecol. 22 :6018–6032.24118668
Bi K. , VanderpoolD., SinghalS., LinderothT., MoritzC., GoodJ.M. 2012. Transcriptome-based exon capture enables highly cost-effective comparative genomic data collection at moderate evolutionary scales. BMC Genomics. 13 :403.22900609
Blom M. 2015. EAPhy: a flexible tool for high-throughput quality filtering of exon-alignments and data processing for phylogenetic methods. PLoS Curr. 7 . doi:10.1371/currents.tol.75134257bd389c04bc1d26d42aa9089f
Bonnet T. , LebloisR., RoussetF., CrochetP. 2017. A reassessment of explanations for discordant introgressions of mitochondrial and nuclear genomes. Evolution. 71 :2140–2158.28703292
Bragg J.G. , PotterS., BiK., MoritzC. 2016. Exon capture phylogenomics: efficacy across scales of divergence. Mol. Ecol. Resour. 16 :1059–1068.26215687
Bryant D. , MoultonV. 2002. NeighborNet: an agglomerative method for the construction of planar phylogenetic networks. In: GuigóR., GusfieldD., editors. International Workshop on algorithms in bioinformatics. Heidelberg: Springer. p. 375–391.
Cahill J.A. , GreenR.E., FultonT.L., StillerM., JayF., OvsyanikovN., SalamzadeR., St. JohnJ., StirlingI., SlatkinM., ShapiroB. 2013. Genomic evidence for island population conversion resolves conflicting theories of polar bear evolution. PLoS Genet. 9 :e1003345.23516372
Camacho C. , CoulourisG., AvagyanV., MaN., PapadopoulosJ., BealerK., MaddenT.L. 2009. BLAST+: architecture and applications. BMC Bioinf. 10 :421.
Card D.C. , ShapiroB., GiribetG., MoritzC., EdwardsS.V. 2021. Museum genomics. Annu. Rev. Genet. 55 :633–659.34555285
Catullo R.A. , KeoghJ.S. 2014. Aridification drove repeated episodes of diversification between Australian biomes: evidence from multi-locus phylogeny of Australian toadlets (Uperoleia: Myobatrachidae). Mol. Phylogenet. Evol. 79 :106–117.24971737
Chavez A.S. , MaherS.P., ArbogastB.S., KenagyG.J. 2013. Diversification and gene flow in nascent lineages of island and mainland North American tree squirrels (Tamiasciurus). Evolution. 68 :1094–1109.
Chifman J. , KubatkoL. 2014. Quartet inference from SNP data under the coalescent model. Bioinformatics. 30 :3317–3324.25104814
Crisp M.D. , LaffanS., LinderH.P., MonroA. 2001. Endemism in the Australian flora. J. Biogeogr. 28 :183–198.
Cruickshank T.E. , HahnM.W. 2014. Reanalysis suggests that genomic islands of speciation are due to reduced diversity, not reduced gene flow. Mol. Ecol. 23 :3133–3157.24845075
Currat M. , RuediM., PetitR.J., ExcoffierL. 2008. The hidden side of invasions: massive introgression by local genes. Evolution 62 :1908–1920.18452573
Edwards S.V. , PotterS., SchmittJ.C., BraggJ.G., MoritzC. 2016. Reticulation, divergence, and the phylogeography-phylogenetics continuum. Proc. Natl. Acad. Sci. U.S.A. 113 :8025–8032.27432956
Eldridge M.D.B. 1997. Taxonomy of rock-wallabies, Petrogale (Marsupialia: Macropodidae). II. An historical review. Aust. Mammal. 19 :113–122.
Eldridge M.D.B. , JohnstonP.G., LowryP.S. 1992. Chromosomal rearrangements in rock wallabies, Petrogale (Marsupialia: Macropodidae). VII. G-banding analysis of Petrogale Petrogale (Marsupialia: Macropodidae). VII. G-banding analysis of Petrogale brachyotis and P. concinna: species with dramatically altered karyotypes. Cytogenet. Cell Genet. 61 :34–39.1505229
Eldridge M.D.B. , PotterS., CooperS.J.B. 2011. Biogeographic barriers in north-western Australia: an overview and standardisation of nomenclature. Aust. J. Zool. 59 :270–272.
Eldridge M.D.B. , PotterS., JohnsonC.N., RitchieE.G. 2014. Differing impact of a major biogeographic barrier on genetic structure in two large kangaroos from the monsoon tropics of Northern Australia. Ecol. Evol. 4 :554–567.25035797
Excoffier L. , FollM., PetitR.J. 2009. Genetic consequences of range expansions. Annu. Rev. Ecol. Evol. Syst. 40 :481–501.
Fenker J. , TedeschiL.G., MelvilleJ., MoritzC. 2021. Predictors of phylogeographic structure among codistributed taxa across the complex Australian monsoonal tropics. Mol. Ecol. 30 :4276–4291.34216506
Ferreira M.S. , JonesM.R., CallahanC.M., FareloL., TolesaZ., SuchentrunkF., BoursotP., MillsL.S., AlvesP.C., GoodJ.M., Melo-FerreiraJ. 2021. The legacy of recurrent introgression during the radiation of hares. Syst. Biol. 70 :593–607.33263746
Figueiró H.V. , LiG., TrindadeF.J., AssisJ., PaisF., FernandesG., SantosS.H.D., HughesG.M., KomissarovA., AntunesA., TrincaC.S., RodriguesM.R., LinderothT., BiK., SilveiraL., AzevedoF.C.C., KantekD., RamalhoE., BrassalotiR.A., VillelaP.M.S., NunesA.L.V., TeixeiraR.H.F., MoratoR.G., LoskaD., SaragüetaP., GabaldónT., TeelingE.C., O’BrienS.J., NielsenR., CoutinhoL.L., OliveiraG., MurphyW.J., EizirikE. 2017. Genome-wide signatures of complex introgression and adaptive evolution in the big cats. Sci. Adv. 3 :e1700299.28776029
Firneno T.J. Jr , O’NeillJ.R., PortikD.M., EmeryA.H., TownsendJ.H., FujitaM.K. 2020. Finding complexity in complexes: assessing the causes of cytonuclear discordance in a problematic species complex of Mesoamerican toads. Mol. Ecol. 29 :3543–3559.32500624
Fontaine M.C. , PeaseJ.B., SteeleA., WaterhouseR.M., NeafseyD.E., SharakhovI.V., JiangX., HallA.B., CatterucciaF., KakaniE., MitchellS.N., WuY., SmithH.A., LoveR.R., LawniczakM.K., SlotmanM.A., RichS., HahnM.W., BesanskyN.J. 2015. Extensive introgression in a malaria vector species complex revealed by phylogenomics. Science. 347 :1258524.25431491
Ford J. 1978. Geographical isolation and morphological and habitat differentiation between birds of the Kimberley and the Northern Territory. Emu - Austral Ornithol. 78 :25–35.
Fraïsse C. , PopovicI., MazoyerC., SpataroB., DelmotteS., RomiguierJ., LoireE., SimonA., GaltierN., DuretL., BierneN., VekemansX., RouxC. 2021. DILS: demographic inferences with linked selection by using ABC. Mol. Ecol. Resour. 21 :2629–2644.33448666
Good J.M. , HirdS., ReidN., DemboskiJ.R., SteppanS.J., Martin-NimsT.R., SullivanJ. 2008. Ancient hybridization and mitochondrial capture between two species of chipmunks. Mol. Ecol. 17 :1313–1327.18302691
Guerrero R.F. , KirkpatrickM. 2014. Local adaptation and the evolution of chromosome fusions. Evolution. 68 :2747–2756.24964074
Guo W. , SunD., CaoY., XiaoL., HuangX., RenW., XuS., YangG. 2022. Extensive interspecific gene flow shaped complex evolutionary history and underestimated species diversity in rapidly radiated dolphins. J. Mamm. Evol. 29 :353–367.
Guo Y. , LongJ., HeJ., LiC.-I., CaiQ., ShuX.-O., ZhengW., LiC. 2012. Exome sequencing generates high quality data in non-target regions. BMC Genomics. 13 :194.22607156
Hahn C. , BachmannL., ChevreuxB. 2013. Reconstructing mitochondrial genomes directly from genomic next-generation sequencing reads – a baiting and iterative mapping approach. Nucleic Acids Res. 41 :e129.23661685
Harrison R.G. , LarsonE.L. 2014. Hybridization, introgression, and the nature of species boundaries. J. Hered. 105 :795–809.25149255
Harrison R.G. , LarsonE.L. 2016. Heterogeneous genome divergence, differential introgression, and the origin and structure of hybrid zones. Mol. Ecol. 25 :2454–2466.26857437
Hey J. , ChungY., SethuramanA., LachanceJ., TishkoS., SousaV.C., WangY. 2018. Phylogeny estimation by integration over isolation with migration models. Mol. Biol. Evol. 35 :2805–2818.30137463
Hibbins M.S. , HahnM.W. 2022. Phylogenomics approaches to detecting and characterizing introgression. Genetics. 220 :iyab173.34788444
Hudson R. 2002. Generating samples under a Wright–Fisher neutral model of genetic variation. Bioinformatics. 18 :337–338.11847089
Huson D.H. , BryantD. 2006. Application of phylogenetic networks in evolutionary studies. Mol. Biol. Evol. 23 :254–267.16221896
Janke A. , XuX., ArnasonU. 1997. The complete mitochondrial genome of the wallaroo (Macropus robustus) and the phylogenetic relationship among Monotremata, Marsupialia, and Eutheria. Proc. Natl. Acad. Sci. U.S.A. 94 :1276–1281.9037043
Jaya F.R. , TannerJ.C., WhiteheadM.R., DoughtyP., KeoghJ.S., MoritzC.C., CatulloR.A. 2022. Population genomics and sexual signals support reproductive character displacement in Uperoleia (Anura: Myobatrachidae) in a contact zone. Mol. Ecol. 31 :4527–4543.35780470
Ji J. , DonavanJ.J., LeachéA.D., YangZ. 2023. Power of Bayesian and heuristic tests to detect cross-species introgression with reference to gene flow in the Tamias quadrivittatus group of North American chipmunks. Syst. Biol. 72 :446–465.36504374
Johnson R.N. , O’MeallyD., ChenZ., EtheringtonG.J., HoS.Y.W., NashW.J., GrueberC.E., ChengY., WhittingtonC.M., DennisonS., PeelE., HaertyW., O’NeillR.J., ColganD., RussellT.L., Alquezar-PlanasD.E., AttenbrowV., BraggJ.G., BrandiesP.A., ChongA.Y., DeakinJ.E., Di PalmaF., DudaZ., EldridgeM.D.B., EwartK.M., HoggC.J., FrankhamG.J., GeorgesA., GillettA.K., GovendirM., GreenwoodA.D., HayakawaT., HelgenK.M., HobbsM., HolleleyC.E., HeiderT.N., JonesE.A., KingA., MaddenD., Marshall GravesJ.A., MorrisK.M., NeavesL.E., PatelH.R., PolkinghorneA., RenfreeM.B., RobinC., SalinasR., TsangarasK., WatersP.D., WatersS.A., WrightB., WilkinsM.R., TimmsP., BelovK. 2018. Adaptation and conservation insights from the koala genome. Nat. Genet. 50 :1102–1111.29967444
Katoh K. , MisawaK., KumaK., MiyataT. 2002. MAFFT: a novel method for rapid multiple sequence alignment based on fast Fourier transform. Nucleic Acids Res. 30 :3059–3066.12136088
Kearns A.M. , JosephL., ToonA., CookL.G. 2014. Australia’s arid-adapted butcherbirds experienced range expansions during Pleistocene glacial maxima. Nat. Comm. 5 :1–11.
Kitchener D.J. 1989. Taxonomic appraisal of Zyzomys (Rodentia, Muridae) with descriptions of two new species from Northern Territory, Australia. Rec. West. Aust. Mus. 14 :331–373.
Köhler F. , CriscioneF. 2015. A molecular phylogeny of camaenid land snails from north-western Australia unravels widespread homoplasy in morphological characters (Gastropoda, Helicoidea). Mol. Phylogenet. Evol. 83 :44–55.25463754
Krosby M. , RohwerS. 2009. A 2000 km genetic wake yields evidence for northern glacial refugia and hybrid zone movement in a pair of songbirds. Proc. Biol. Sci. 276 :615–621.18986973
Laver R.J. , DoughtyP., OliverP.M. 2018. Origins and patterns of endemic diversity in two specialized lizard lineages from the Australian Monsoonal Tropics (Oedura spp.). J. Biogeog. 45 :142–153.
Layton K.K.S. , CarvajalJ.I., WilsonN.G. 2020. Mimicry and cytonuclear discordance in nudibranchs: new insights from exon capture phylogenomics. Ecol. Evol. 10 :11966–11982.33209263
Lin H.-Y. , HaoY.-J., LiJ.-H., FuC.-X., SoltisP.S., SoltisD.E., ZhaoY.-P. 2019. Phylogenomic conflict resulting from ancient introgression following species diversification in Stewartia s.l. (Theaceae). Mol. Phylogenet. Evol. 135 :1–11.30802596
Linck E. , EpperlyK., Van ElsP., SpellmanG.M., BrysonR.W., McCormackJ.E., Canales-Del-CastilloR., KlickaJ. 2019. Dense geographic and genomic sampling reveals paraphyly and a cryptic lineage in a classic sibling species complex. Syst. Biol. 68 :956–966.31135028
Linnen C.R. , FarrellB.D. 2007. Cytonuclear discordance is caused by rampant mitochondrial introgression in Neodiprion (Hymenoptera: Diprionidae) sawflies. Evolution 61 :1417–1438.17542850
Mallet J. , BesanskyN., HahnM.W. 2016. How reticulated are species? Bioessays. 38 :140–149.26709836
Marchi N. , ExcoffierL. 2020. Gene flow as a simple cause for an excess of high-frequency-derived alleles. Evol. Appl. 13 :2254–2263.33005222
Martin S.H. , DaveyJ.W., SalazarC., JigginsC.D. 2019. Recombination rate variation shapes barriers to introgression across butterfly genomes. PLoS Biol. 17 :e2006288.30730876
Melo-Ferreira J. , BoursotP., SuchentrunkF., FerrandN., AlvesP.C. 2005. Invasion from the cold past: extensive introgression of mountain hare (Lepus timidus) mitochondrial DNA into three other hare species in northern Iberia. Mol. Ecol. 14 :2459–2464.15969727
Melo-Ferreira J. , VilelaJ., FonsecaM.M., da FonsecaR.R., BoursotP., AlvesP.C. 2014. The elusive nature of adaptive mitochondrial evolution of an arctic lineage prone to frequent introgression. Genome Biol. Evol. 6 :886–896.24696399
Meyer M. , KircherM. 2010. Illumina sequencing library preparation for highly multiplexed target capture and sequencing. Cold Spring Harbor Protocols. 6 :pdb.prot5448. doi:10.1101/pdb.prot5448
Minh B.Q. , SchmidtH.A., ChernomorO., SchrempfD., WoodhamsM.D., von HaeselerA., LanfearA. 2020. IQ-TREE 2: new models and efficient methods for phylogenetic inference in the genomic era. Mol. Biol. Evol. 37 :1530–1534.32011700
Morales H. , PavlovaA., AmosN., MajorR., KilianA., GreeningC., SunnucksP. 2018. Concordant divergence of mitogenomes and a cytonuclear gene cluster in bird lineages inhabiting different climates. Nat. Ecol. Evol. 2 :1258–1267.29988164
Moritz C. , FujitaM.K., RosauerD., AgudoR., BourkeG., DoughtyP., PalmerR., PepperM., PotterS., PrattR., ScottM., TonioneM., DonnellanS. 2016. Multilocus phylogeography reveals nested endemism in a gecko across the monsoonal tropics of Australia. Mol. Ecol. 25 :1354–1366.26671627
Moritz C.C. , PrattR.C., BankS., BourkeG., BraggJ.G., DoughtyP., KeoghJ.S., LaverR.J., PotterS., TeasdaleL.C., TedeschiL.G., OliverP.M. 2018. Cryptic lineage diversity, body size divergence, and sympatry in a species complex of Australian lizards (Gehyra). Evolution. 72 :54–66.29067680
Oliver P.M. , JollyC.J., SkipwithP.L., TedeschiL.G., GillespieG.R. 2020. A new velvet gecko (Oedura: Diplodactylindae) from Groote Eylandt, Northern Territory. Zootaxa. 4779 :438–450.
Oliver P.M. , LaverR.J., MartinsF.D.M., PrattR.C., HunjanS., MoritzC.C. 2016. A novel hotspot of vertebrate endemism and an evolutionary refugium in tropical Australia. Divers. Distrib. 23 :53–66.
Park S. , ParkS.J. 2020. Large-scale phylogenomics reveals ancient introgression in Asian Hepatica and new insights into the origin of the insular endemic Hepatica maxima. Sci. Rep. 10 :16288.33004955
Payseur B.A. , RiesebergL.H. 2016. A genomic perspective on hybridization and speciation. Mol. Ecol. 25 :2337–2360.26836441
Peñalba J.V. , JosephL., MoritzC. 2019. Current geography masks dynamic history of gene flow during speciation in northern Australian birds. Mol. Ecol. 28 :630–643.30561150
Pepper M. , KeoghJ.S. 2014. Biogeography of the Kimberley, Western Australia: a review of landscape evolution and biotic response in an ancient refugium. J. Biogeog. 41 :1443–1455.
Peter B.M. , SlatkinM. 2013. Detecting range expansions from genetic data. Evolution 67 :3274–3289.24152007
Peter B.M. , SlatkinM. 2015. The effective founder effect in a spatially expanding population. Evolution. 69 :721–734.25656983
Pfeifer B. , WittelsbürgerU., Ramos-OnsinsS.E., LercherM.J. 2014. PopGenome: an efficient Swiss army knife for population genomic analyses in R. Mol. Biol. Evol. 31 :1929–1936.24739305
Phuong M.A. , BiK., MoritzC. 2017. Range instability leads to cytonuclear discordance in a morphologically cryptic ground squirrel species complex. Mol. Ecol. 26 :4743–4755.28734067
Portik D.M. , SmithL.L., BiK. 2016. An evaluation of transcriptome-based exon capture for frog phylogenomics across multiple scales of divergence (Class: Amphibia, Order: Anura). Mol. Ecol. Resour. 16 :1069–1083.27241806
Potter S. , BraggJ.B., BlomM.P.K., DeakinJ., KirkpatrickM., EldridgeM.D.B., MoritzC. 2017. Chromosomal speciation in the genomics era: disentangling phylogenetic evolution of rock-wallabies. Front. Genet. 8 :10.28265284
Potter S. , BraggJ.B., PeterB.M., BiK., MoritzC. 2016. Phylogenomics at the tips: inferring lineages and their demographic history in a tropical lizards, Carlia amax. Mol. Ecol. 25 :1367–1380.26818481
Potter S. , BraggJ.G., TurakulovR., EldridgeM.D.B., DeakinJ., KirkpatrickM., EdwardsR.J., MoritzC. 2022. Limited introgression between rock-wallabies with extensive chromosomal rearrangements. Mol. Biol. Evol. 39 :msab333.34865126
Potter S. , CloseR.L., TaggartD.A., CooperS.J.B., EldridgeM.D.B. 2014. Taxonomy of rock-wallabies, Petrogale (Marsupialia: Macropodidae). IV. Multifaceted study of the brachyotis group identifies additional taxa. Aust. J. Zool. 62 :401–414.
Potter S. , CooperS.J., MetcalfeC.J., TaggartD.A., EldridgeM.D.B. 2012a. Phylogenetic relationships of rock-wallabies Petrogale (Marsupialia: Macropodidae) and their biogeographic history within Australia. Mol. Phylogenet. Evol. 62 :640–652.22122943
Potter S. , EldridgeM.D.B., TaggartD.A., CooperS.J.B. 2012b. Multiple biogeographical barriers identified across the monsoon tropics of northern Australia: phylogeographic analysis of the brachyotis group of rock-wallabies. Mol. Ecol. 21 :2254–2269.22417115
Potter S. , XueA.T., BraggJ.G., RosauerD.F., RoycroftE.J., MoritzC. 2018. Pleistocene climatic changes drive diversification across a tropical savanna. Mol. Ecol. 27 :520–532.29178445
Powney G.D. , GrenyerR., OrmeC.D.L., OwensI.P.F., MeiriS. 2010. Hot, dry and different: Australian lizard richness is unlike that of mammals, amphibians and birds. Glob. Ecol. Biogeogr. 19 :386–396.
Reeves J.M. , BostockH.C., AyliffeL.K., BarrowsT.T., De DeckkerP., DevriendtL.S., DunbarG.C., DrysdaleR.N., FitzsimmonsK.E., GaganM.K., GriffithsM.L., HaberleS.G., JansenJ.D., KrauseC., LewisS., McGregorH.V., MooneyS.D., MossP., NansonG.C., PurcellA., van der KaarsS. 2013. Palaeoenvironmental change in tropical Australasia over the last 30,000 years – a synthesis by the OZ-INTIMATE group. Quat. Sci. Rev. 74 :97–114.
Rheindt F.E. , EdwardsS.V. 2011. Genetic introgression: an integral but neglected component of speciation in birds. Auk 128 :620–632.
Rieseberg L.H. 2001. Chromosomal rearrangements and speciation. Trends Ecol. Evol. 16 :351–358.11403867
Rivera D. , PratesI., FirnenoT.J.Jr, RodriguesM.T., CaldwellJ.P., FujitaM.K. 2022. Phylogenomics, introgression, and demographic history of South American true toads (Rhinella). Mol. Ecol. 31 :978–992.34784086
Rosauer D.F. , BlomM.P.K., BourkeG., CatalanoS., DonnellanS., GillespieG., MulderE., OliverP.M., PotterS., PrattR.C., RaboskyD.L., SkipwithP.L., MoritzC. 2016. Phylogeography, hotspots and conservation priorities: an example from the Top End of Australia. Biol. Conserv. 204 :83–93.
Rosauer D.F. , ByrneM., BlomM.P.K., CoatesD.J., DonnellanS., DoughtyP., KeoghJ.S., KinlochJ., LaverR.J., MyersC., OliverP.M., PotterS., RaboskyD.L., Afonso SilvaA.C., SmithJ., MoritzC. 2018. Real-world conservation planning for evolutionary diversity in the Kimberley, Australia, sidesteps uncertain taxonomy. Conserv. Lett. 11 :e12438.
Roux C. , FraïsseC., RomiguierJ., AnciauxY., GaltierN., BierneN. 2016. Shedding light on the grey zone of speciation along a continuum of genomic divergence. PLoS Biol. 14 :e2000234.28027292
Sarver B.A.J. , HerreraN.D., SneddonD., HunterS.S., SettlesM.L., KronenbergZ., DemboskiJ.R., GoodJ.M., SullivanJ. 2021. Diversification, introgression, and rampant cytonuclear discordance in Rocky Mountains Chipmunks (Sciuridae: Tamias). Syst. Biol. 70 :908–921.33410870
Seixas F.A. , BoursotP., Melo-FerreiraJ. 2018. The genomic impact of historical hybridization with massive mitochondrial DNA introgression. Genome Biol. 19 :91.30056805
Sharman G.B. , CloseR.L., MaynesG.M. 1990. Chromosomal evolution, phylogeny and speciation of rock wallabies (Petrogale: Macropodidae). Aust. J. Zool. 37 :351–363.
Shelley J.J. , SwearerS.E., AdamsM., DempsterT., Le FeuvreM.C., HammerM.P., UnmackP.J. 2018. Cryptic biodiversity in the freshwater fishes of the Kimberley endemism hotspot, northwestern Australia. Mol. Phylogenet. Evol. 127 :843–858.29953937
Shelley J.J. , SwearerS.E., DempsterT., AdamsM., Le FeuvreM.C., HammerM.P., UnmackP.J. 2020. Plio-Pleistocene sea-level changes drive speciation of freshwater fishes in north-western Australia. J. Biogeog. 47 :1727–1738.
Singhal S. 2013. De novo transcriptomic analyses for non-model organisms: an evaluation of methods across a multi-species data set. Mol. Ecol. Resour. 13 :403–416.23414390
Singhal S. , DerryberryG.E., BravoG.A., DerryberryE.P., BrumfieldR.T., HarveyM.G. 2021. The dynamics of introgression across an avian radiation. Evol. Lett. 5 :568–581.34917397
Singhal S. , MoritzC. 2012. Testing hypotheses for genealogical discordance in a rainforest lizard. Mol. Ecol. 21 :5059–5072.22989358
Slatkin M. 2005. Seeing ghosts: the effect of unsampled populations on migration rates estimated for sampled populations. Mol. Ecol. 14 :67–73.15643951
Smith K.L. , HarmonL.J., ShooL.P., MelvilleJ. 2011. Evidence of constrained phenotypic evolution in a cryptic species complex of agamid lizards. Evolut. Int. J. Org Evolution 65 :976–992.
Stamatakis A. 2014. RAxML version 8: a tool for phylogenetic analysis and post-analysis of large phylogenies. Bioinformatics. 30 :1312–1313.24451623
Sunnucks P. , HalesD.F. 1996. Numerous transposed sequences of mitochondrial cytochrome oxidase I-II in aphids of the genus Sitobion (Hemiptera: Aphididae). Mol. Biol. Evol. 13 :510–524.8742640
Swofford D.L. 2003. PAUP*. Phylogenetic analysis using parsimony (* and other methods). Version 4. Sunderland: Sinauer Associates.
Taylor R.S. , BramwellA.C., Clemente-CarvalhoR., CairnsN.A., BonierF., DaresK., LougheedS.C. 2021. Cytonuclear discordance in the crowned-sparrows, Zonotrichia atricapilla and Zonotrichia leucophrys. Mol. Phylogenet. Evol. 162 :107216.34082131
Thomas O. 1904. On a new rock-wallaby from north-west Australia. Novitates Zoologicae. 11 :365–366.
Toews D.P.L. , BrelsfordA. 2012. The biogeography of mitochondrial and nuclear discordance in animals. Mol. Ecol. 21 :3907–3930.22738314
Torkian B. , HannS., PreisnerE., NormanR.S. 2020. BLAST-QC: automated analysis of BLAST results. Environ. Microbiome 15 :15.33902722
Wang L. , LiuS., YangY., MengZ., ZhuangZ. 2022. Linked selection, differential introgression and recombination rate variation promote heterogeneous divergence in a pair of yellow croakers. Mol. Ecol. 31 :5729–5744. doi: 10.1111/mec.16693 36111361
Wascher M. , KubatkoL. 2021. Consistency of SVDQuartets and maximum likelihood for coalescent-based species tree estimation. Syst. Biol. 70 :33–48.32415974
Wilson C.C. , BernatchezL. 1998. The ghost of hybrids past: fixation of arctic charr (Salvelinus alpinus) mitochondrial DNA in an introgressed population of lake trout (S. namaycush). Mol. Ecol. 7 :127–132.
Wolf J.B.W. , EllegrenH. 2017. Making sense of genomics islands of differentiation in light of speciation. Nat. Rev. Genet. 18 :87–100.27840429
Wu C.-I. 2001. The genic view of the process of speciation. J. Evol. Biol. 14 :851–865.
