Differences in how the genome is packaged – chromosome variation – have long been known to influence how genetic variation is transmitted and redistributed within and among populations (Darlington, 1958; White, 1973). Today, with increasing availability of high quality genome assemblies, the capacity for genome-scale resequencing and the tools of molecular cytogenetics and population- and phylo-genomic analysis, we are returning to a whole-genome perspective on evolution. At the same time, evidence from genome comparisons is revealing that reticulate evolution, including introgression across distantly related species, is far more common in animals than previously thought (Mallet et al., 2016) with implications for gene-tree – species-tree discordance (Edwards et al., 2016). This directs attention to the potential for differing extents of introgression within and outside rearranged regions of the genome (Noor and Bennett, 2009; Crawford et al., 2015).

Here, we revisit the early history of thinking about how chromosomal rearrangements (CRs) affect population and speciation processes. We then highlight a case study that emphasizes how a combined knowledge of genome architecture and genomic sequence divergence is important for understanding the history of CRs in association with speciation and being able to assess whether gene flow is reduced in rearranged regions. Finally, we return to broader themes, considering what combinations of evidence and theory are necessary to gain a holistic understanding of how chromosome change can promote incipient divergence and ultimately translate into species diversification.

The rock-wallaby (Petrogale) system has been considered a classical model for chromosomal speciation due to extensive chromosome repatterning, combined with their habitat specialization (rocky environments) which causes populations to be isolated and small (King, 1993). Rock-wallabies are medium sized marsupials (1–12 kg) that inhabit complex rocky areas distributed across continental Australia and some offshore islands, (Eldridge, 2008). A strong propensity for isolation among disjunct rocky habitats (e.g., Pope et al., 1996; Hazlitt et al., 2006) is thought to increase their rate of speciation and contribute to the fixation of novel CRs (Eldridge and Close, 1993). Petrogale includes 17 recognized species corresponding to 23 chromosomal taxa (Supplementary Table 1; Figure 1). This is the most chromosomally diverse genus of marsupials, which in general have a conserved karyotype across all five Australasian and American super-families (2n = 14; Rofe and Hayman, 1985; Hayman, 1990; see O’Neill et al., 1999; Graves and Renfree, 2013). Macropodids (kangaroos and wallabies) show variable karyotypes (see O’Neill et al., 1999), but the ancestral macropodid 2n = 22 karyotype is only found in Petrogale (P. lateralis, P. persephone, P. rothschildi, and P. xanthopus; Eldridge et al., 1992a) and Thylogale (Pademelons). The 2n = 22 macropodid ancestral karyotype is itself derived from the widespread 2n = 14 marsupial karyotype by a series of fissions (Rofe, 1979; Hayman, 1990). CRs are extensive across Petrogale, and range from simple to complex. A majority of the rearrangements are Robertsonian fusions, but there is also cytogenetic evidence for inversions and centric shifts (centromeric transpositions) (Eldridge et al., 1989, 1990, 1991, 1992b; Eldridge and Close, 1992, 1993). In addition, the X chromosome is frequently variable in morphology among taxa and also sometimes within taxa (Eldridge and Close, 1997). The highest chromosomal diversity occurs in groups that are sympatric or parapatric (brachyotis and penicillata groups; Eldridge et al., 1992b; Eldridge and Close, 1993; Figure 1), and these are also the most speciose. This pattern matches the predictions of chromosomal speciation models (see Faria and Navarro, 2010) and the recombination suppression model. This could reflect yet another scenario where fixation rate of chromosomal rearrangements correlates with parapatry and sympatry – suggesting adaptation and divergent selection could be a dominant process driving fixation (e.g., Hooper and Price, 2015).

The brachyotis group of rock-wallabies includes four species and five sub-species distributed across northwestern Australia which have the most complex rearrangements found in Petrogale (Supplementary Table 2) and large amounts of centromeric constitutive heterochromatin not present in other Petrogale (Maynes, 1989; Sharman et al., 1990; Eldridge et al., 1992b; Eldridge and Close, 1993). Within the lateralis group, there are two chromosomal races and three sub-species (Eldridge et al., 1991; Figure 1). These races/sub-species are recently diverged and are distinguished by single autosomal rearrangements or fusions (Eldridge and Close, 1993, 1997; Figure 1). However, the most interesting group are the recently diverged (∼0.5–2.7 mya; Potter et al., 2012a) Queensland penicillata group taxa. Six parapatric species display extensive variation in karyotypes ranging from simple to complex – including fusions, inversions and centric shifts (see Eldridge and Close, 1993). Early research was driven by cytogenetic analyses (reviewed in Eldridge and Close, 1997) and captive breeding experiments that showed evidence of reproductive isolation including infertile male hybrids and reduced fertility of female hybrids (Eldridge and Close, 1992). This resulted in the description of three new species (Eldridge and Close, 1992) and a focus on the role of chromosomal variation in speciation. Meiotic irregularities, including problems with more extensive rearrangements and X-autosome associations have been reported (Close et al., 1996) – patterns also seen in model systems (e.g., Mus). Recent genetic analysis of the penicillata group using microsatellites and mitochondrial DNA (mtDNA) found extensive sharing of alleles between some of the most chromosomally divergent species (Potter et al., 2015). This could be a consequence of introgression or ILS. Further analysis of nuclear markers across the genome is required to assess the genomic divergence between these species and assess if speciation with gene flow is occurring between these taxa, or if more complex interactions between genomic architecture and genic divergence is at play. The characteristics of this genus, specifically their rapid radiation and extensive chromosome variation, make them a valuable model for understanding chromosome evolution and speciation.

Phylogenetic analyses of the rock-wallabies have not previously included representatives of all 23 chromosomal taxa, nor have they been able to resolve phylogenetic relationships, particularly among the more recently evolved species within the penicillata group (Campeau-Péloquin et al., 2001; Potter et al., 2012a). This, in addition to evident homoplasy of rearrangements (see Eldridge and Close, 1993), has precluded tracing the evolution of chromosomal changes. The phylogenetic relationship of P. xanthopus and P. purpureicollis has also been difficult to resolve (see Eldridge et al., 1991; Eldridge and Close, 1993; Potter et al., 2012a), which has hindered interpreting chromosome evolution as these taxa retain the ancestral chromosome number.

Here, we report results from targeted capture for ∼2000 exons from two individuals per taxon to resolve the relationships across the genus (all supplementary material and methods are outlined in Supplementary Datasheet 1). These data allow us to understand the evolution of chromosomes in this group. In particular, they provide insight into phylogenetic and sequence divergence signals of discordance across the X, rearranged and non-rearranged chromosome arms that could reflect effects of CRs on gene flow (see Supplementary Table 3). We focus on sets of concatenated loci, rather than individual gene trees as individual exons have low phylogenetic resolution at this scale. While it is desirable to use multispecies coalescent approaches (e.g., ^∗BEAST and ASTRAL), such programs are confounded by introgression across non-sister data and are therefore unsuitable for this system (see Solis-Lemus et al., 2016). Hence, we explore multispecies coalescent network approaches that allow for introgression (see below; reviewed in Nakhleh, 2013; Edwards et al., 2016). We expect to find discordant phylogenies and divergence levels between these categories of loci, particularly for the recent radiation of Queensland taxa.

The phylogenetic relationships amongst taxa using the entire dataset of 1961 exons and ∼1 million bp firmly resolves, for the first time, relationships within Petrogale (Figure 2). The brachyotis group with the most extensive chromosomal rearrangements is also phylogenetically basal, which is consistent with previous genetic data (Campeau-Péloquin et al., 2001; Potter et al., 2012a). The xanthopus chromosomal group is paraphyletic. P. rothschildi forms the sister taxon to the lateralis and penicillata groups, where as P. persephone and P. xanthopus form a well supported monophyletic group. There is deep divergence between all three taxa that retain the ancestral karyoptype. The cytogenetically conservative lateralis group has similar branch lengths among taxa to those among species within the more chromosomally diverse penicillata group. Despite recent speciation, each taxon within the lateralis and penicillata groups is monophyletic, albeit with lower support for P. lateralis lateralis, P. l. West Kimberley race and P. assimilis. Relationships amongst some of the most closely related taxa are also not strongly resolved, but we note that monophyly of P. mareeba and P. sharmani is consistent with their chromosomal evolutionary history, since both share a derived fusion between chromosomes 5 and 10. Although the phylogenetic position of P. purpureicollis has been previously unresolved (Sharman et al., 1990; Eldridge et al., 1991; Campeau-Péloquin et al., 2001; Potter et al., 2012a), these new data strongly resolve it as sister to the penicillata group. The sub-species of both P. brachyotis and P. xanthopus (with no known chromosomal differences) did not form monophyletic lineages, suggesting recent divergence or some nuclear gene flow.

Given a well-resolved phylogeny, we can now investigate the history of chromosome change in the genus. Using parsimony mapping of CRs, as identified by G-banding (Eldridge and Close, 1993) we were able to resolve ancestral nodes where rearrangements occurred, in particular a single origin of CR 7a (a = acrocentric). However, for some chromosomes, we could not distinguish between different hypotheses (see Figure 3; Supplementary Table 2). Apparent multiple independent origins of the 3a, 4a, 4sm (sm = submetacentric) and 5i (i = inversion) rearrangements suggest there could be regions of the genome susceptible to rearrangement processes (“hotspots”), which have also been implicated in chromosome change in other macropodids (Bulazel et al., 2007). This highlights the potential for convergent evolution of rearrangements, including chromosomal fusions (e.g., 6 and 10 fusion), inversions and centromere shifts. The alternate hypotheses are multiple reversals to an ancestral chromosome morphology, or the introgression of chromosomes between taxa, but further analysis (e.g., sequencing of breakpoints) is required to distinguish between them. In addition to parallel evolution, the same chromosomes are involved in fusion events in different taxa (e.g., 5, 6, 9, and 10). Cell culture experiments using mitomycin C to induce centric fusions showed that chromosome 10 fused most frequently (Eldridge and Johnston, 1993). Despite all chromosomes being involved in fusions in this experiment, the higher frequency of chromosome 10 fusions in vitro matches the larger proportion of chromosome 10 being involved in fusions in the wild in Petrogale (five out of eight fusions). This together with higher frequencies of breakpoints from gamma radiation in chromosomes 5, 6 and 10 (Eldridge and Johnston, 1993) further support the notion of a nonrandom process of CR.

Next, we partitioned the sequenced exons into autosomal-non-rearranged (N = 140 exons; 75,296 bp), autosomal-rearranged (N = 160 exons; 36,168 bp), and the X chromosome (N = 21 exons; 8,951 bp). This approach was motivated by the expectation that autosomes will have less phylogenetic signal than the X because of their higher gene flow rates and larger N_e (see Supplementary Table 3 for mapped loci). Mean divergence between the four chromosomal groups varies across the X chromosome, rearranged and non-rearranged autosomes (Figure 4). When accounting for differences in sequence length, we find that the X has reduced diversity compared with the autosomes, although only slightly compared to the rearranged autosomes. In other mammalian systems (e.g., apes – Nam et al., 2014) several selective sweeps have created large regions of low diversity on the X chromosome. Assessment of more loci along the X is necessary to explore if the lower diversity on the X within Petrogale is associated with selection. The effects of small N_e on the X will need to be assessed to distinguish between selective sweeps and neutral models of evolution. However, given that sex chromosomes contribute disproportionately to post-zygotic isolation in many taxa (e.g., Drosophila, Presgraves, 2008; Flycatchers, Saether et al., 2007), and evidence of the Petrogale system conforming to Haldane’s rule, we would expect more loci on the X to conform to the true species phylogeny than autosomal loci (as argued by Fontaine et al., 2015).

Rearranged and non-rearranged chromosomes show similar levels of divergence in the brachyotis and penicillata groups, which are among the most chromosomally diverse Petrogale (Figure 5). Conversely, rearranged chromosomes showed the greatest increase in divergence relative to non-rearranged chromosomes in the lateralis group, which has relatively few rearrangements. These observations do not tend to support the hypothesis that CRs are a primary cause of speciation. Instead, the CRs may have fixed when species were already isolated. Alternatively, the existing data may be inadequate to distinguish between these hypotheses. Further work is needed.

We expect clearer phylogenetic resolution in rearranged than in non-rearranged regions of the genome, particularly within the penicillata group (Figure 5). Neither autosomal phylogenetic reconstructions are able to resolve the relationships of the lateralis and penicillata groups. The phylogeny based on rearranged regions does however, separate P. herberti and P. penicillata from the remaining taxa, as well as P. coenensis and P. godmani, compared to the non-rearranged chromosomes. Both autosomal phylogenies resolve the brachyotis and xanthopus groups, but the rearranged phylogeny places the xanthopus group as basal instead of the brachyotis group. By contrast, the X loci resolve nodes deeper in the tree but appears to lack enough information to resolve all of the internal relationships of the lateralis and penicillata groups. The X however, does generally group individuals of a taxon together, unlike the autosomes (Figure 2).

We then asked how often the individuals sampled from each species formed a monophyletic group for the different subsets of exons. If rearrangements result in reduced introgression, we expect to see higher concordance in concatenated loci from rearranged than non-rearranged chromosome arms. Further, as the X chromosome is frequently found to be resistant to gene flow, we also expected higher congruence across the X-linked loci. On average, the rearranged chromosomes have greater monophyly of taxa than the non-rearranged for the penicillata group, supporting our hypothesis of higher concordance. We do, however, find the opposite pattern for the lateralis group (Supplementary Table 4). This may be because the lateralis group has fewer rearranged loci. Overall it had lower concordance of monophyletic individuals compared to the penicillata group. The X chromosome had the greatest average monophyly. This may result from a smaller N_e of the X, faster divergence of the X (e.g., Charlesworth et al., 1987), or greater divergent selection on the X. The evolution of the sex chromosomes needs further investigation.

Analysis of all mitochondrial coding genes (12 genes; 11,373 bp), albeit still a single linkage group, reveals some strong conflicts between mitochondrial and nuclear evolutionary history (Figure 2). Previously, it has been highlighted that introgression, retained ancestral polymorphism (or ILS) has resulted in paraphyletic species complexes within both the brachyotis group (see Potter et al., 2012a,b) and the penicillata group (Briscoe et al., 1982; Bee and Close, 1993; Potter et al., 2015). This is the first analysis using all coding genes across the mitochondrial genome, providing the most phylogenetic information and highlight discrepancies to the nuclear phylogeny, in particular – the placement of P. persephone and P. xanthopus as basal branches; the paraphyly of brachyotis group taxa; and lack of monophyly for P. assimilis, P. coenensis and P. godmani within the penicillata group. These inconsistencies further highlight areas of potential reticulation in the history of Petrogale; in particular, the potential for repeated episodes of introgression in the history of the genus.

We then explored patterns of reticulation among chromosomal groups. We based these analyses on concatenated nuclear loci, effectively ignoring coalescent variance in gene trees. We first used an exploratory statistical approach (Neighbor-Net in Splits tree – Bryant and Moulton, 2002; Huson and Bryant, 2006, based on average distances). The results suggested reticulation amongst P. coenensis and P. godmani, between P. assimilis, P. mareeba and P. sharmani, as well as between P. lateralis lateralis and P. lateralis West Kimberley race (Figure 6). The brachyotis group does not indicate any strong evidence of reticulation between species. The penicillata group results match previous results of microsatellites and mtDNA, which inferred gene flow between these taxa, even those with complex CRs (Potter et al., 2015).

We next used an approach based on the multispecies network coalescent (PhyloNet – Than et al., 2008; Yu et al., 2013, 2014). We find evidence of reticulation for all three chromosomal groups (Figure 6; Supplementary Table 5). Based on our analysis of up to three reticulation events, the results support 2–3 reticulation nodes in the brachyotis group, which suggests historical introgression may explain the discordance between mtDNA and nuclear loci. This included reticulation at nodes of ancestral branches in the brachyotis group. For both the lateralis and penicillata comparisons there was support for a single reticulation node. Within the lateralis group, the reticulation node involved P. lateralis West Kimberley race. Reticulation could reflect ILS between P. l. lateralis and P. l. hacketti or introgression with P. l. lateralis. In the penicillata group, the analyses were less concordant, and the reticulation node included P. assimilis and P. sharmani for the independent analysis. Both network analyses reflect greater reticulation across all three species and include introgression with P. mareeba and P. assimilis, as well as ILS amongst the three species. In all three chromosomal group comparisons, reticulation is evident between species with CRs (fusions and centric shifts). Further work is necessary to disentangle the effects of ILS from introgression, which will require further sampling in rearranged vs. non-rearranged regions of chromosomes.

Understanding the genome, the physical position of loci and how they interact is crucial in interpreting the evolutionary history of organisms. Our results highlight that if certain loci taken alone without any context of chromosome structure can yield completely different results and a misunderstanding of the mechanisms involved in reproductive isolation. We are still in the early stages of understanding the physical location of loci in this non-model system and as further work allows mapping of loci to chromosomes and regions of rearrangements, we will be better able to test for recombination suppression and reduced gene flow in these regions compared to the non-rearranged chromosomes. With the addition of X chromosome loci, we may be able to further understand the genetic variation on the X which itself has a high rate of intrachromosomal rearrangement in this system. The combined effects of faster X divergence and reduced recombination may affect speciation in this genus. Below, we outline the necessary steps to progress this non-model system – a framework we believe useful for any non-model speciation research.

The case study above, for a fascinating yet “non-model” system, indicates both the promise and the challenges of arriving at a holistic understanding of how CRs affect divergence and adaptation. Hybrid infertility and recombination suppression (adaptive) models have largely been treated as exclusive, but in systems with complex rearrangements, both could well be in play (Faria and Navarro, 2010; Garagna et al., 2014). Further, such hypotheses should not be treated as exclusive of accumulation of genic incompatibilities, such as large X-effects (Presgraves, 2008), effects of recessive X-linked incompatibility alleles on sterility of the heterogametic sex (Haldane’s rule; Haldane, 1922; Turelli, 1998; Presgraves, 2008, 2010) or cytonuclear incompatibilities (Turelli and Moyle, 2007). In addition to new evidence across diverse systems, we need further development of theory that incorporates these co-occurring processes (Feder and Nosil, 2009; Faria and Navarro, 2010). We should also consider extended models with meiotic drive, genomic conflict and disruption of epigenetic programming (Brown and O’Neill, 2010).

Considering all these interacting processes, we suggest a hierarchy of questions and associated requirements that progress from simple to challenging (Figure 7). The first question is the history of speciation and CRs, as we have now begun to resolve for the rock-wallabies. We move on to testing the effects of CRs on the suppression of gene flow, and then to testing mechanisms that may have suppressed gene flow. The scale and diversity of data needed expands as we move from pattern to process. Of particular note is the need to develop new theory and inference methods that can distinguish the different processes by integrating genomic data, recombination landscapes, meiotic irregularities, segregation patterns, and epigenetics. Such methods involve complicated parameter estimation and may need to rely on approaches such as approximate Bayesian computation (ABC) to compare competing models (see Beaumont, 2010). Those models, in turn, will likely employ multispecies coalescence that includes speciation and demographic history, as well as the evolution of CRs (e.g., coalescent models for CRs – Guerrero et al., 2012a,b; Ayala et al., 2013; Peischl et al., 2013). The quality of genome assemblies needed increases as we move from pattern to process. This is crucial to remove errors that can influence inference of rearrangements as well as recombination rate changes.

In the following, we explore three key issues that have emerged: (i) what do we know of differential origins of new (or recurrent) CRs; (ii) how do CRs interact with the epigenetic programming of genes and chromosome segments; and (iii) what new theory and inference methods do we need to exploit comparative genomic data in the context of individual and combined effects of CRs.

Integration of cytogenetic, genomic, and epigenetic data using a holistic approach will be crucial to improving our understanding of how genomic architecture influences and potentially drives reproductive isolation amongst organisms. The non-model system highlighted here has value for many reasons aside from its recent origin and chromosomal diversity. These include: the effect of reversed recombination rate between the sexes in marsupials relative to most other eutherians (males > females), X and Y chromosomes do not pair during meiosis (they lack the pseudoautosomal region), and epigenetic mechanisms vary significantly from eutherians. This one non-model system illustrates how biological variation can provide valuable contrasts to model systems. The growth of genomic and computational technology is opening new vistas on fundamental questions about how genomic architecture influences evolution.

Raw sequencing reads associated with this study are available at NCBI Sequence Read Archive (http://www.ncbi.nlm.nih.gov/sra; BioProject PRJNA360868; BioSamples SAMN06233288-SAMN06233338). A Dryad Repository (http://dx.doi.org/10.5061/dryad.mm856) contains code that was used in data analysis and Sequence Read Archive accession details.

This study did not use live animals and therefore does not require animal ethics approval. Samples used were tissues from the Australian Museum collection. Animal research at the Australian Museum and in Australia is governed by the “Australian code for the care and use of animals for scientific purposes.” The Code only requires a research project to be approved by an ethics committee if the research involves live animals.

SP had substantial contribution to the conception, acquisition, design, analysis and interpretation of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work. JB had substantial contribution to the conception, design, analysis and interpretation of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work. MB had substantial contribution to the analysis of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work. JD had substantial contribution to the conception, acquisition and design of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work. MK had substantial contribution to the conception of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work. ME and CM had substantial contribution to the conception, acquisition, design and interpretation of the work, drafting the work, final approval of the version to be published and agreement to be accountable for all aspects of the work.

The handling Editor declared a shared affiliation, though no other collaboration, with one of the authors JD and states that the process nevertheless met the standards of a fair and objective review. The other authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.