Meiotic drive shapes rates of karyotype evolution in mammals.

Chromosome number is perhaps the most basic characteristic of a genome, yet generalizations that can explain the evolution of this trait across large clades have remained elusive. Using karyotype data from over 1000 mammals, we developed and applied a phylogenetic model of chromosome evolution that links chromosome number changes with karyotype morphology. Using our model, we infer that rates of chromosome number evolution are significantly lower in species with karyotypes that consist of either all bibrachial or all monobrachial chromosomes than in species with a mix of both types of morphologies. We suggest that species with homogeneous karyotypes may represent cases where meiotic drive acts to stabilize the karyotype, favoring the chromosome morphologies already present in the genome. In contrast, rapid bouts of chromosome number evolution in taxa with mixed karyotypes may indicate that a switch in the polarity of female meiotic drive favors changes in chromosome number. We do not find any evidence that karyotype morphology affects rates of speciation or extinction. Furthermore, we document that switches in meiotic drive polarity are likely common and have occurred in most major clades of mammals, and that rapid remodeling of karyotypes may be more common than once thought. 2019 The Authors Evolution published by Wiley Periodicals, Inc. on behalf of The Society for the Study of Evolution.

Chromosome number is one of the first genomic characteristics identified (Goldsmith 1919; Painter 1921), and by the 1950s chromosome number had been investigated in several hundred species (Makino 1951). Despite this early start, the dynamics of chromosome number evolution have been recalcitrant to generalizations across large taxonomic scales. Researchers have argued that chromosome number is a direct agent driving the evolution of other traits (Blackmon et al. 2015), an indirect target of selection (Sherman 1979; Ross et al. 2015), or a mostly neutral trait whose dynamics are dominated by drift (Bickham and Baker 1979; Rockman and Rowell 2002). Early cytogeneticists suggested that changes in chromosome number might be an engine of speciation (Stebbins 1971; Árnason 1972; White 1973, 1978), but depending on the underlying assumptions, theoretical and empirical studies have shown mixed support for changes in chromosome number driving speciation (Lande 1979; Hedrick 1981; Rieseberg 2001; Guerrero and Kirkpatrick 2014; Freyman and Höhna 2018).

Meiotic drive is perhaps one of the best explanations for the rapid evolution of chromosome number. In almost all organisms female meiosis is asymmetric—the primary oocyte with four copies of each chromosome will pass only one copy to the haploid gamete, and the remaining three will be eliminated with the polar bodies during meiosis. If there are consistent differences between the spindles that lead to the polar bodies versus the spindles that lead to the ovum, then any chromosome that can preferentially attach to the spindle leading to the ovum will be overrepresented in gametes and will “drive” to high frequency in a population (Rhoades 1942; Sandler and Novitski 1957). In fact, limited data from mouse, chicken, and human indicate that oogenesis in these species favors the transmission of certain chromosome morphologies. Specifically, oogenesis in mouse preferentially passes monobrachial (acrocentric or telocentric) chromosomes to the egg while oogenesis in both human and chicken preferentially passes bibrachial (metacentric) chromosomes to the egg (Pardo‐Manuel de Villena and Sapienza 2001). In the human and chicken cases, we expect oogenesis to favor a decrease in chromosome number by preferentially passing two monobrachial chromosomes as a fused bibrachial chromosome. We use the term fusion for verbal simplicity, but simple fusion between two chromosomes is unlikely to describe the process at the chromosomal level—instead, a Robertsonian translocation likely with some loss of nonessential DNA underlies this process (Schubert and Lysak 2011). In contrast, in the case of mice, we expect oogenesis to favor an increase in chromosome number by preferentially passing fissioned chromosomes. In this case, our terminology is well aligned with the process at the chromosomal level as breaking in the centromeric region followed by de novo telomere gain is well documented (Harrington and Greider 1991; Tsujimoto et al. 1999). This tendency for meiotic drive to favor one type of chromosome morphology is termed meiotic drive polarity.

These observations lead us to hypothesize a taxon with a karyotype having homogeneous chromosome morphology has reached a point where their meiotic drive polarity is matched with the morphology of their chromosomes, leading to relative stability in chromosome number. In contrast, we hypothesize that a taxon with a karyotype having heterogeneous chromosome morphology may be in the process of transitioning from one type of karyotype morphology to the other, leading to a high rate of chromosome number evolution through either fusions or fissions. Based on this hypothesis, when all chromosomes of a taxon have a single morphology we will refer to the karyotype as “matched,” and a taxon with a mixture of chromosome morphologies will be termed “mismatched.” Under this hypothesis, we predict that the rate of chromosome number evolution in taxa with mismatched karyotype morphology should be higher than in those with matched karyotype morphology. An alternative hypothesis is that species with mismatched karyotypes may simply lack this type of meiotic drive. In this case, rates of chromosome number evolution may not be consistently higher for mismatched karyotypes than for matched karyotypes.

Despite the large amount of karyotype data available, the macroevolutionary relationship between karyotype morphology and rate of chromosome number evolution has not been explored. To address this gap in our understanding, we developed a probabilistic model of chromosome evolution that allows for a binary trait to evolve on a phylogenetic tree and for each state of that trait to have unique rates of chromosome number evolution, speciation, and extinction. In our analyses, the two states of the binary trait are mismatched karyotype morphology (both mono‐ and bibrachial chromosomes) versus matched karyotypes (all mono‐ or all bibrachial chromosomes). This coding is intended to separate species where a recent change in the meiotic drive polarity may be driving change in chromosome number (mismatched) from those where the meiotic drive polarity may act to increase karyotype stability (matched). We used this model to evaluate the evolution of chromosome number across a supertree phylogeny of 1059 mammals. To examine whether the observed pattern is consistent within the major mammalian clades, we also conducted subclade analyses of chromosome number evolution using independently inferred phylogenies for Cetartiodactyla, Marsupialia, Carnivora, Chiroptera, and Primates. We find that the rate of chromosome number evolution is dramatically increased in taxa with mismatched karyotypes. We suggest that the patterns we observe are best explained by frequent switches in the polarity of meiotic drive across mammals.

Methods

DATA COLLECTION AND SCORING

We collected chromosome number and morphology data from existing cytogenetic compilations (Pardo‐Manuel de Villena and Sapienza 2001; Tree of Sex Consortium 2014). We define chromosome number as the haploid chromosome count of females. We use the female haploid chromosome count and morphology because males contain Y chromosomes that are not exposed to the effects of female meiotic drive. Furthermore, females produce gametes with an equal number of chromosomes while males in species with multiple X or Y chromosomes (e.g., XXY sex chromosome systems) produce haploid gametes with different numbers of chromosomes. We classify chromosome morphology as either monobrachial (acrocentric or telocentric) or bibrachial (metacentric). In cases where original compilations lacked a description of chromosome morphology, the morphology reported in the Atlas of Mammalian Chromosomes was used (O'Brien et al. 2006). We were able to collect data for 1060 species (Table S1). Species names were updated to match accepted taxonomy (Wilson and Reeder 2005). We used a stringent threshold to discretize the morphology of chromosomes present in a karyotype. We scored the karyotype as “matched” if all chromosomes were of one type (monobrachial or bibrachial) and “mismatched” if both types of chromosomes were present (but see alternative coding below).

To analyze the evolution of chromosome number in a phylogenetic framework, we searched for the most complete and recently published time‐calibrated species‐level phylogenies. For our mammal‐wide analysis, we used the supertree from Fritz et al. (2009) as modified into a pseudo‐posterior distribution by random resolution of polytomies and branch length adjustments (Kuhn et al. 2011). The mammal supertree includes virtually all species, but many are placed solely based on taxonomy, and the supertree construction methods leave artifacts in the branch lengths. To test if our conclusions are robust to those concerns, and to assess their consistency within mammals, we also conducted our analyses in five large subclades with phylogenies based on molecular data. For the analyses of Primates, Cetartiodactyla, and Carnivora, we downloaded 100 trees from a posterior distribution using the 10kTrees website; these trees contained 84, 102, and 69 species, respectively (Arnold et al. 2010). The species on these trees represent 17, 30, and 24 percent of extant species in each clade, respectively. The Cetartiodactyla trees had branch lengths in units of substitutions, so we transformed them into ultrametric trees using the chronos function in the R Package APE (Paradis 2013). This function is a penalized likelihood approach, and we assumed correlated variation in substitution rate among branches and a smoothing parameter of 1. The Chiroptera tree was downloaded from Dryad and contained 155 tips, or 12% of extant species (Shi and Rabosky 2015a, b). The tree for Marsupialia was supplied by the authors of the original study and contained 40 tips, or 12% of extant species (Mitchell et al. 2014). Phylogenetic uncertainty could not be incorporated into our analysis for Chirpotera and Marsupialia because only a single maximum likelihood tree was available for each.

MODEL DESCRIPTION

Our model builds on existing dynamic phylogenetic models of chromosome number evolution (Mayrose et al. 2010; Glick and Mayrose 2014), and, like the model of Zenil‐Ferguson et al. (2017), it incorporates a binary trait that can affect rates of chromosome number change. Unlike those models, our model can additionally allow for the binary trait to affect rates of speciation and extinction, using the mathematical framework introduced by Maddison et al. (2007). Our model implementation is available in an R package called chromePlus (github.com/coleoguy/chromePlus). The chromePlus package uses the existing likelihood functions for a multistate character that may or may not affect rates of speciation and extinction (Lewis 2001; FitzJohn 2012), and it automatically constrains them to represent biologically realistic models of chromosome number evolution. It is built on the diversitree package (FitzJohn 2012), which provides these likelihood functions and additional amenities such as allowing for uncertainty in the binary character state, sampling incompleteness, and MCMC (Markov Chain Monte Carlo) or likelihood maximization functions to estimate model parameters.

A graphical schematic of our model is shown in Fig. 1, depicting the ten possible state transitions. While in binary state 1, a lineage with i chromosomes may either increase chromosome number to (ascending dysploidy, with rate γ1), reduce chromosome number to (descending dysploidy, δ1), or increase chromosome number to 1.5i or 2i (demi‐polyploidy, ν1, and polyploidy, ρ1, respectively). Likewise, a lineage in binary state 2 with i chromosomes may increase chromosome number to (ascending dysploidy, γ2), reduce chromosome number to (descending dysploidy, δ2), or increase chromosome number to 1.5i or 2i (demi‐polyploidy, ν2, and polyploidy, ρ2, respectively). Additionally, a lineage that is currently in one of the binary states can transition to the other state with no change in chromosome number (transitions from state 1 to state 2 at rate q 12, or from state 2 state 1 at rate q 21). Because chromosome number is limited to whole numbers, demi‐polyploidy as described above would only be possible for even values of i. As in previous work, we assume that demi‐polyploidy when i is odd leads to either of the closest whole numbers with equal probability. The state‐dependent diversification version of our model has four additional parameters: speciation and extinction rates in binary state 1 (λ1 and μ1, respectively), and speciation and extinction in binary state 2 (λ2 and μ2, respectively).

Figure 1

Schematic representation of the model for joint evolution of a binary character and chromosome number. The model depicts the possible transitions of a lineage that is in character state 1 or 2 and has i haploid number of chromosomes. For each binary state, there can be as many as four transition rates: descending dysploidy δ leading to chromosomes, ascending dysploidy γ leading to chromosomes, demi‐polyploidy ν leading to 1.5i chromosomes, and polyploidy ρ leading to 2i chromosomes (with subscripts on those rates indicating they can differ for binary states 1 and 2). Transition rates between the binary character states are q 12 and q 21. The chromePlus package can implement this model of state changes, and it can include state‐dependent speciation and extinction effects. In this case, the binary character states 1 and 2 have speciation rates λ1 and λ2 and extinction rates μ1 and μ2.

The above is the general description of chromePlus, but our analyses here do not include polyploidization and thus omit the four ploidy‐dependent parameters (ν1, ρ1, ν2, ρ2). Although polyploidy is an important force in the evolution of chromosome number in many groups (Glick et al. 2016; Zenil‐Ferguson et al. 2017), there is little evidence for it in mammals. The one possible exception is the red viscacha rat (Tympanoctomys barrerae) from Argentina. However, the evidence for and against polyploidy in this species is arguably inconclusive (Gallardo et al. 1999; Svartman et al. 2005; Teta et al. 2014; Evans et al. 2017). We chose to exclude this species from our dataset for three reasons. First, estimating a transition rate in a model for an event that has only occurred once is at best difficult and at worst may bias other parameter estimates. Second, it allows us to estimate the same set of rates (i.e., excluding and ) in all subclade analyses. Finally, we are attempting to understand broad patterns across all mammals rather than one exceptional event.

SIMULATION STUDY

We conducted a set of four simulation tests to evaluate the efficacy of our method. In tests one and two, we worked with four sizes of phylogenies: 50, 100, 200, and 500 extant taxa. For each size, we simulated 100 phylogenies, generated using a birth‐death model with speciation rate three times greater than extinction rate. Each tree was then scaled to a root age of one, to ease interpretation of the transition rate parameter values. In test one, transitions of the binary character were symmetric, with , and we evaluated four scenarios for chromosome number evolution. In the base model simulation, rates of chromosome gain and loss were the same, with and . In the three alternative models, we set the rates of chromosome number change in binary state 2 (γ2 and δ2) to be higher than those in state 1 by factors of 2, 5, or 10. In test two, transitions of the binary character were asymmetric, with and , and we evaluated the same four scenarios for chromosome number evolution as we did in test one. These rate parameters were chosen so that the smallest trees with 50 taxa would still exhibit variation in chromosome number (median range of eight under a model with equal rates of chromosome evolution) and at least 10% of taxa in each binary state. The remaining scenarios result in datasets with increasingly large effect sizes that should be easier to detect during testing. Because acceptance of a trait simulation was conditioned on the presence of both states of the binary character at the end of the simulation, 2–5% of simulations were discarded. This procedure provided a total of 3200 simulated datasets. In both the symmetric and asymmetric scenarios, we used the base model to measure the rate of type I error: falsely reporting different rates of chromosome number evolution associated with the states of the binary character.

Models of trait evolution can produce biased results if one of the states is associated with higher diversification (Maddison 2006). To assess this effect on inferences about chromosome number evolution, we performed several additional sets of simulations. Test three evaluates performance on simulated phylogenies under three scenarios. The first scenario is a test of the base model where one binary state is associated with a threefold increase in speciation (), but the rates of character and chromosome number evolution are the same in both binary states (, ). In the second scenario, one binary state was again associated with a threefold increase in speciation (), and this same binary state also had a threefold increase in the rate of chromosome number evolution (). In the third scenario, one binary state was again associated with a threefold increase in speciation (), but it was then the slower diversifying binary state that had a threefold increase in the rate of chromosome number evolution (). We simulated trees with 200 extant tips under each of these models, and again scaled each to a root age of one. We fit our model to each simulated dataset in a Bayesian framework using the mcmc function in the R package diversitree. For the equal diversification scenarios we estimated parameter values under a model with the six transition parameters: δ1, γ1, δ2, γ2, q 12, and q 21. For the differential diversification scenarios, we inferred parameters under two models: both contained the six transition parameters above, but one allowed for different rates of speciation and extinction in the two binary states, and the other constrained speciation and extinction rates to be equal. Preliminary analyses indicated that a small fraction of MCMC chains would mix poorly and that in these cases some parameter estimates were excessively high (and biologically unrealistic). To eliminate this behavior, we placed a uniform prior from 0 to 200 on each parameter. The MCMC chain was initialized with parameter values drawn from a uniform distribution between 0 and 5. Preliminary analysis indicated that MCMC chains reached convergence within 500 steps. The MCMC was allowed to run for an additional 500 steps, and parameter estimates are based on this postburnin portion of the chain. For each dataset, the inference was performed on the original simulated tree.

In test four, we focused on the impact of an empirical phylogeny that contains a strong signal of diversification rate heterogeneity. We used a time‐calibrated phylogeny of 73 cetaceans that has been shown to generate a high false‐positive rate using standard BiSSE analyses (Steeman et al. 2009; Rabosky and Goldberg 2015). On this phylogeny, again scaled to a root age of one, we simulated two scenarios with equal rates of chromosome evolution in both binary states. The first set was generated with low rates of chromosome number evolution and transitions in the binary state (, , and ) and the second set had high rates (, , and ). For each set of rates we performed 100 simulations. The MCMC chain was initialized with parameter values drawn from a uniform distribution between 0 and 5. Preliminary analysis indicated that MCMC chains reached convergence within 100 steps. The MCMC was allowed to run for an additional 300 steps, and parameter estimates are based on this postburnin portion of the chain. Scripts for the simulations analyses are provided as a supplement on DRYAD: https://doi.org/10.5061/dryad.rg5s170.

Because our focal question is whether one of our binary states has different rates of chromosome number evolution (fusions or fissions), we report our results in terms of a mean rate difference statistic, . For each postburnin sample we calculated as In addition to estimating the magnitude of this effect, we also reduced it to a simple test of the motivating hypothesis by comparing the 95% credible interval of (i.e., the 95% highest posterior density) with zero. If the entire 95% credible interval of is positive, we interpret it as support for a higher rate of chromosome evolution when the binary character is in state 2. If the entire 95% credible interval of is negative, we interpret it as support for a higher rate of chromosome evolution when the binary character is in state 1. Otherwise, we conclude there is not support for the rate of chromosome evolution being different in the two states of the binary character. Using the statistic we are able to combine results from analyses across a posterior distribution of trees accounting for phylogenetic uncertainty.

ANALYSIS OF MAMMALIAN DATA

Diversification analysis

To determine whether there is any difference in the net diversification rate of taxa with matched versus mismatched karyotypes, we first fit a simplified model to our full mammal dataset, considering the matched or mismatched binary trait but not chromosome number. With this model we estimated the transition rate between matched and mismatched karyotypes as well as the speciation and extinction rate in each of these states, using the BiSSE model (Maddison et al. 2007) as implemented in the R package diversitree (FitzJohn 2012). Parameter estimates were sampled using MCMC with 1000 steps and a burnin of 500 steps on each of the 100 mammalian supertrees. We calculated the net diversification rate for each binary state as speciation minus extinction. To test for a significant difference in diversification, we calculated the 95% credible interval of the difference of net diversification between taxa with matched karyotypes and taxa with mismatched karyotypes. If the credibile interval overlaps zero, this suggests no significant difference in rates of diversification. If it does not, it will be important to allow for state‐dependent diversification in our analysis of chromosome number evolution.

Chromosome rate analysis

We analyzed our dataset of 1059 mammalian karyotypes on a pseudo‐posterior distribution of 100 supertrees. Based on the results from the BiSSE analyses (reported below), we did not need to use the version of our model that allows rates of speciation and extinction to depend on the binary karyotype trait. We therefore obtained estimates of six parameters: rates of chromosome number increase (γ1 and γ2), rates of chromosome number decrease (δ1 and δ2), and rates of change in karyotype state (q 12 and q 21). An uninformative, unbounded improper prior (assuming all non‐negative values are equally likely) was applied for all parameters. The MCMC chain was initialized with parameter values drawn from a uniform distribution from 0 to 10, which is broad but biologically reasonable. Preliminary analysis indicated that MCMC chains reached convergence quickly in most cases (usually less than 500 generations), but that a small proportion of runs got “stuck” in a certain region of the parameter space. To avoid this problem, four independent MCMC runs on each tree were completed (run to 1000 steps) and compared to ensure that they had converged. The last 250 generations of one successful MCMC chain per tree were retained, and these were combined from all 100 trees to form the posterior sample. We then calculated the statistic described above to assess if the binary karyotype trait is associated with different rates of chromosome number evolution.

Subclade analyses were performed similarly. For Cetartiodactyla and Primates, parameters were estimated with an MCMC chain that ran for 5000 generations with the first 10% discarded as burnin on each of 100 trees, and postburnin samples were combined from all trees. For Carnivora, preliminary analyses indicated that the MCMC chain was sampling two different regions of the parameter space. To ensure there was not a failure of convergence or mixing, we ran longer chains (10,000 generations each with 10% treated as burnin) and inspected their movement in parameter space. After removing burnin which was typically less than 500 generations, assessed by eye, we confirmed that movement between the two peaks of was not extremely rare: there was an average of 23.5 jumps between the two peaks per tree. For Chiroptera and Marsupialia, only single trees were available. Parameters were estimated with an MCMC that ran for 10,000 generations with a burnin of 1,000 generations. All trees were rescaled to unit length for the model‐fitting, but the rates we report for empirical analyses are all transformed to absolute time units (i.e., per million years). For this time‐scaling, we used the median root age estimate for each tree from the Time Tree of Life database (Kumar et al. 2017).

Results

In this work, we developed the chromPlus model, which allows studying the possible interplay between a binary character trait and rates of chromosome‐number evolution, while additionally allowing for diversification rate differences between the two character states. Our simulations indicate that our method is not prone to false positives and has reasonable power (Fig. 2). As expected, the propensity to detect an association between the binary character and rates of chromosome number evolution increased with tree size and effect size. For example, in test one, if the rates differed by a factor of two between the binary states, the effect was detected in only 8% of simulations with 50‐tip trees, but in 85% of simulations with 500‐tip trees. Increasing the rates ratio to a factor of 5 made the effect detectable in 47% of 50‐tip trees and 100% of 500 tip trees. Comparing tests one and two, we found no difference in the power or false positive rate of our method when the rates of transitions in the binary character were either symmetric or asymmetric (Fig. 2A,B).

Figure 2

Simulation analysis of chromePlus model. In (A–C) the vertical axis is the proportion of simulations in which the states of the binary character were inferred to have different rates of chromosome number evolution. In (A) and (B), results are shown for test one (symmetric transitions in the binary trait) and test two (asymmetric transitions in the binary trait) with tree sizes of 50, 100, 200, and 500 taxa. The rates ratio describes the difference in the rates of chromosome number evolution in the two states of the binary character. A rates ratio of 1, bottom lines, illustrate the false‐positive rate. The three other lines in each plot represent power with different magnitudes of rate differences. In (C) are results for test three where speciation in one binary state was three times higher than in the other, and trees contained 200 taxa. Rates of chromosome number evolution are equal (left), greater for the state with higher diversification (center), or greater for the state with lower diversification (right). Labels on the horizontal axis indicate whether the inference model was with or without state‐dependent diversification. In (D) are results for test four where simulations with equal rates of chromosome evolution were performed on an empirical phylogeny. Labels on the horizontal axis indicate whether the simulation was for low or high rates of chromosome number and binary trait evolution. The vertical axis indicates the proportion of simulations in which the states of the binary character were inferred to have different rates of diversification or chromosome number evolution.

Test three, which allowed for state‐dependent diversification, also yielded satisfactory results (Fig. 2C). We found no evidence of an increased false positive rate when simulating under a scenario where one character state affects speciation rate but not rates of chromosome evolution. Whether or not the inference model allowed for different speciation and extinction rates associated with the binary character, the false‐positive rate was 6%, and the estimated rates of chromosome number evolution were largely the same. Failing to account for state‐dependent diversification did, however, bias the transition rate estimates, as could be expected (Maddison 2006). Specifically, when the inference model constrained equal speciation rates and equal extinction rates, parameter estimates were larger and more variable for the transition into the binary state with higher speciation rate (Fig. S1). In the next two simulation scenarios, the states of the binary character had different rates of both diversification and chromosome number evolution. When simulating a scenario where state 1 of the binary character had a threefold increase in both diversification and chromosome number evolution, our power to detect a significant difference in rates of chromosome evolution was reasonable (power of 0.55 when the inference model constrained speciation and extinction to be constant, and 0.56 when the model allowed for separate rates of speciation and extinction in the two binary states; Fig. 2C, center). Higher power to detect a significant difference in rates of chromosome evolution was obtained when state 1 of the binary character had a threefold increase in diversification rates but a threefold decrease in rates of chromosome number evolution (0.7 when the inference model constrained speciation and extinction to be constant and 0.66 when the model allowed for separate rates of speciation and extinction in the two binary states; Fig. 2C, right). In our state‐dependent speciation and extinction simulations, the median number of taxa with the lower diversifying binary state was only 49 while the median number of taxa with the alternative state was 151. Therefore we suggest that the slight increase in power when the slower diversifying binary trait is associated with higher rates of chromosome evolution is likely a result of having sufficient chromosome number changes on both binary state backgrounds to reliably infer rates.

Because many empirical phylogenies contain diversification rate heterogeneity unrelated to the traits of interest, we also explored the performance of our approach on a time‐calibrated phylogeny of cetaceans (test four). We found that it was unlikely to incorrectly infer different rates of chromosome number evolution for the two binary states (false‐positive rates of 0% for both the low and high rates of chromosome number evolution). Although this component of inference is robust, we did find an elevated rate of false positives for inference of differential diversification associated with our binary trait. Differences in net diversification rate were falsely associated with the binary trait in 40% and 70% of simulations for, respectively, the low and high rates of chromosome evolution and transitions in the binary state. These results underscore the caution that must be applied when interpreting results from state dependent speciation and extinction models like chromePlus and BiSSE which it is built upon (Rabosky and Goldberg 2015).

Distribution of karyotype morphologies

Mammal karyotype morphologies display a distribution consistent with widespread meiotic drive. Specifically, the distribution of chromosome number in mammals is unimodal and approximately normal, but the distribution of karyotype morphologies is bimodal, with most species having either all bibrachial chromosomes or all monobrachial chromosomes (Fig. 3). It is also important to note that species with homogeneous karyotype morphology exhibit a wide range of chromosome numbers, indicating that chromosome number is not primarily dictated by karyotype morphology. To varying degrees, the bimodal distribution of karyotype morphologies is replicated at lower taxonomic levels within mammals (Fig. S2). Mapping karyotype morphology onto phylogenies reveals that species with variation in chromosome morphology are often spread broadly across the tree (Fig. S3)

Figure 3

Distribution of karyotype morphology in mammals. The horizontal axis indicates the proportion of chromosomes that are monobrachial, and the vertical axis provides the total count of haploid chromosomes in a female genome. The margins show histograms for each of these characters, emphasizing the unimodal distribution of chromosome number and the bimodal distribution of karyotype morphology. Values for each species are plotted, and where they coincide, point size indicates the number of species. The apparent banding pattern results from the fewer possible values of the proportion as the haploid number decreases.

Our analysis using the BiSSE model indicated a slightly higher rate of net diversification (speciation minus extinction) for species with mismatched karyotypes than for species with matched karyotypes. However, our estimates of these two diversification rates were largely overlapping (Fig. S4A), and the 95% credible interval of the difference in diversification rates overlapped zero (−2.73 to 5.73; Fig. S4B). These results suggest that mammals with matched or mismatched karyotypes do not have detectably different net diversification rates. Based on these results, we performed the subsequent analyses of chromosome number without accounting for state‐dependent diversification.

Rates of chromosome number evolution

Our analysis of all mammalian data indicates a strong pattern of greater rates of chromosome number evolution in species with mismatched karyotypes (Fig. 4A). However, the use of a tree where polytomies have been randomly resolved is expected to increase transition rate estimates and total numbers of transitions (Rabosky 2015). We address this potential positive bias in the rate estimate for the mammal‐wide analysis by conducting a separate analysis for each of the five mammalian subclades. These subclade analyses further account for rate heterogeneity and could further examine whether our results obtained for the entire mammalian clade are consistent across the major subclades. Chiroptera, Marsupialia, Primates, and Cetartiodactyla each show clear support for higher rates of chromosome number evolution in taxa with mismatched karyotypes, with the 95% credible interval of entirely positive (Fig. 4B–D, F).

Figure 4

Differences in rates of chromosome number evolution for taxa with mismatched and matched karyotypes. The difference, , is calculated as the mean rate in taxa with mismatched karyotypes minus the mean rate in taxa with matched karyotypes. The 95% credible interval of each distribution is indicated by the matching line beneath the distribution. Panel (A) shows results for all mammals, and each subsequent panel shows the results from the independent analysis of a subclade, as indicated on top of each panel.

Carnivora present a more complex result because of bimodality in the posterior distributions, causing the 95% credible interval of to include zero (Fig. 4E). This fails the simple test of whether rates of chromosome number evolution differ in taxa with matched or mismatched karyotypes. However, the posterior distribution of drops dramatically around zero, and the 95% credible interval of the absolute value of excludes zero (Fig. S5). Examining our results reveals that 70.5% of the posterior distribution of is greater than zero and 29.5% is less than zero. We examined whether some trees were biased toward positive or negative values of but found that 100% of trees sampled both regions ( and ) during the postburnin portion of the analysis. Finally, we found that while positive estimates are associated with slightly higher posterior probabilities, the distribution of posterior probabilities for MCMC samples with negative and positive estimates are largely overlapping (Fig. S6). We interpret these results as support for a difference in rates, but that the data can be explained by having high rates for either binary state. In contrast to the mammal‐wide analysis and the other four subclades, we thus conclude that Carnivora provides only weak support for faster rates in species with mismatched karyotypes.

Taken together, the results obtained for the five mammal subclades are largely consistent with our analysis of all mammals, and they should not be affected by the possible supertree artifact described above. We also calculated the mean of the posterior distribution and credible interval for each of the individual transition rates included in the statistic, and these individual rates estimates are in concordance with the results from the analysis of (Table S2). As part of our subclade analyses, we also estimated the rates of transition between matched and mismatched karyotypes. The rate that taxa switch from a matched karyotype to a mismatched karyotype should represent the rate that taxa experience switches in meiotic drive polarity (whether meiosis favors transmission of bibrachial or monobrachial chromosomes to the gametes). We found the highest rate of polarity switching in Cetartiodactyla, where the mean waiting time for a transition was 10.8 million years (Fig. 5A). The lowest rate of polarity switching was in Primates, where the median rate was almost an order of magnitude longer at 90.9 million years. We did not find a strong correlation between the rate of karyotype switching and the proportion of taxa with mismatched karyotypes (Fig. 5B). Instead, clades with low rates exhibit the full range of variation in the proportion of taxa with mismatched karyotype morphologies.

Figure 5

Meiotic drive polarity switching. (A) The rate of transition from matched karyotype to mismatched karyotype. (B) The relationship between rate of karyotype switching and the proportion of taxa with mismatched karyotypes.

Character coding

One potential concern with the application of our model to the mammal dataset is that the binary trait coding of matched or mismatched is fundamentally linked to chromosome number evolution, which could lead to identifiability issues. For instance, in Akodon boliviensis the karyotype is made up of 38 monobrachial chromosomes, and this could evolve into a karyotype with 19 bibrachial chromosomes through 19 successive fusion events. However, the first and last of these 19 fusion events will also lead to a change in the binary trait, as the lineage evolves from matched to mismatched and back to matched. To explore whether transitions in these extreme karyotypes affect our inference of transition rates, we investigated an alternative threshold for discretization of karyotype morphology into our binary states of matched and mismatched. Exploratory analyses were performed in Primates and Chiroptera where we relaxed the karyotype morphology coding, including as matched any karyotype that was at least 90% bibrachial or 90% monobrachial. Under the stringent coding approach, 26% (22 of 84) of primates had matched karyotypes. Using the relaxed coding approach, 38% (32 of 84) of primates had matched karyotypes. A smaller proportion of taxa were impacted by the alternative coding in Chiroptera. Under the stringent coding approach, 34% (52 of 155) of bats had matched karyotypes. Using the relaxed coding approach, 41% (63 of 155) of bats had matched karyotypes. To evaluate the impact of coding alternatives, we compared the posterior distribution of the rate difference statistic, , under both schemes (Fig. S7). We found minor and opposite effects in the two groups examined. In Primates, the rate difference was somewhat larger with the relaxed coding, and in Chiroptera it was somewhat larger with the stringent coding. However, in both clades the 95% credible interval of these two distributions were largely overlapping, and in no case did it include 0. Like the primary results presented above, analysis under the relaxed coding scheme also provided support for a higher mean rate of chromosome evolution in taxa with mismatched karyotypes. This suggests that our inferences are robust to minor differences in the codings of karyotype morphology.

Discussion

Our findings suggest that meiotic drive may be a widespread force in shaping chromosome number evolution. The striking difference in rates of chromosome number evolution between taxa with matched and mismatched karyotypes suggests that any attempts to identify the role of alternative forces like population structure or selection (Bengtsson 1980; Petitpierre 1987; Martinez et al. 2015; Ross et al. 2015; Martinez et al. 2016) in large‐scale karyotype evolution should allow for the possible impact of meiotic drive.

Many theoretical models have invoked chromosomal evolution as a driver of speciation. We have identified a significantly higher rate of chromosome evolution in taxa with mismatched karyotypes but no increase in speciation associated with this binary state. Taken together, this suggests that chromosome number evolution is unlikely to be a dominant force in the generation of extant mammal biodiversity. However, the application of a model that divides chromosome number into cladogenetic and anagenetic chromosome number change (Freyman and Höhna 2018) may reveal a role for chromosomal speciation in some cases.

The high rates of chromosome number change in taxa with mismatched karyotypes are generally consistent with the rapid changes of karyotype morphology that have been identified in some mouse populations (Britton‐Davidian et al. 2000), suggesting that if a meiotic polarity change occurs at certain times, it could create an effective barrier to gene flow in some cases. For instance, if a species is characterized by a patchy distribution and one deme experiences a transition in meiotic polarity, this deme may fix a sufficiently altered karyotype that subsequent gene flow with other demes is unlikely. Indeed, the rapid karyotypic divergence documented in rodents (Burt et al. 2009; Garagna et al. 2014) suggests that this may be a promising clade for further investigation of this process, although a large species‐level phylogeny will be needed for an analysis like ours.

Recent work in mice has revealed some of the molecular processes that are responsible for meiotic drive (Chmátal et al. 2014). Briefly, it has been shown that the polarity of meiotic drive is based on the relative strength of centromeres, where strength is the ability to express kinetochore proteins and interact with spindle fibers. In this experimental system, it was shown that a fusion with the same centromere strength could be either favored or disfavored depending on the genetic background that it was segregating within. With this understanding, how do we interpret the rates of polarity switching in mammals? Polarity switching may represent a possible side effect of centromeric drive (Henikoff et al. 2001; Fishman and Saunders 2008), occasionally shifting to a point where chromosome morphologies that are not present in the current karyotype are favored if present during meiosis. Traditionally, when closely related lineages hybridize, differential introgression across the genome has been attributed to either positive selection on favored alleles or purifying selection on deleterious alleles (Nance et al. 1990; Li et al. 2010). However, our results suggest that meiotic preference may be widespread and should be considered as a possible explanation for unusual patterns of introgression. Our estimates of polarity switching rates suggest that the rate of this process is considerably higher in some clades. Interestingly, the highest rate of switching in karyotype morphology was in Cetartiodactyla (median waiting time 10.9 myr), and this is also the clade with the highest proportion of mismatched karyotypes (Fig. 5). It remains unclear what characteristic of the genome might explain the marked differences in the frequency of polarity switching in this clade.

Our analysis of carnivore data highlights some of the challenges of applying complex comparative methods. The posterior distribution of rate estimates encompasses two distant regions of parameter space that have very similar posterior probability, leading to difficulty in mixing and convergence (Fig. S6 and Fig. S8). This difficulty illustrates the importance of examining trace files to ensure that individual MCMC runs are not sampling a single local optimum. Our solution was to extend the length of the MCMC chains to allow for sufficient sampling of both regions. However, a more comprehensive solution would be to apply a metropolis coupled MCMC strategy where multiple chains would search the parameter space simultaneously (Gilks and Roberts 1996). One of the chains is the “cold” chain and behaves normally, while several other chains would be heated (the likelihood ratio is raised to a power of a temperature coefficient). This has the effect of flattening the likelihood surface, thus allowing the heated chains to explore a broader area of parameter space. Then at defined intervals, the cold chain is compared to one of the heated chains and a swap can occur so that the cold chain begins exploring a new area of parameter space. Though not yet implemented, this is certainly an area of future development. These problems are not unique to Bayesian approaches. Maximum likelihood inference must also be conducted in a fashion that accounts for the possibility of local optima. The best solution in this framework is to perform multiple fittings with randomly drawn starting values to ensure that the maximum likelihood estimates are consistent and are globally optimal.

Increasingly large phylogenies provide the opportunity to test ever more complex and nuanced models of trait evolution. However, as phylogenies grow larger they cover such large spans of time that it becomes ever more unlikely that a single model or set of parameters describes the evolution of a trait (Felsenstein 2004). For this reason, models like ours that allow for parameter values to vary across the tree are becoming increasingly important to our understanding of the dynamics of trait evolution. Because our analysis effectively combines taxa whose meiotic drive increases chromosome number with those whose meiotic drive decreases chromosome number, we might be underestimating the true rate of chromosome number change. We do not expect this effect to drive our conclusion, however, because it would affect both of the binary states. Nonetheless, several extensions of our model could provide better understanding of the evolution of genome structure at the level of karyotypes. First, based on the molecular mechanisms described above, it seems likely that meiotic drive polarity is controlled by an underlying continuous variable (e.g., average centromere strength). This would fit well with an application of Sewall Wright's (1934) threshold model as extended into a comparative framework (Felsenstein 2012; Revell 2014). Using the threshold approach, a model could be developed where lineages are able to transition into a state without meiotic drive rather than assuming that either type of meiotic drive should operate. This would be a valuable extension that might help us understand how common meiotic drive is in nature. Additionally, a more nuanced view of the karyotype could track chromosome arms (fundamental number) rather than simply chromosome number. The need for this model is based on two factors. Abundant evidence from groups as distant as flies and humans indicates that mutations changing the fundamental number are likely of a different nature than those that change the chromosome number (Matthey 1945; Sved et al. 2016). Furthermore, our data collection suggests that at least two separate processes are occurring. In fact, if we fit a simple linear model to our data (i.e., the interior of Fig. 3), where chromosome number is predicted by the proportion of monobroachial chromosomes in a genome, we find that we have an adjusted R‐squared of 27%. However, if only fusions and fissions altered the haploid number, these variables should be perfectly correlated (i.e., a genome with 100% monobroachial chromosomes should have a chromosome number twice that of a genome with 100% bibrachial chromosomes). The lack of a perfect correlation suggests that additional processes affect these characteristics. A promising future direction would thus be the application of an expanded model that could track fusions and fissions of existing arms as well as gains of new arms and of whole chromosomes through aneuploidy or other remodeling processes.

In conclusion, we have demonstrated that the rate of chromosome number evolution in mammals is associated with the morphology of the karyotype. We argue that the most parsimonious explanation is that karyotypes with both mono‐ and bibrachial chromosomes are transient in nature, and that meiotic drive will eventually push these intermediate karyotypes toward either purely monobrachial or purely bibrachial. Finally, our analysis indicates that the polarity of meiotic drive has switched on many occasions, leading to evolution of both types of chromosome morphology in most mammal subclades.

Associate Editor: M. R. Servedio

Handling Editor: C. Ané

Figure S1. Estimate of transition rate into higher diversification state under two models.

Figure S2. Histograms of karyotype morphologies across all mammal species and within subclades.

Figure S3. Chromosome number and karyotype morphology in mammals.

Figure S4. Diversification parameter estimates.

Figure S5. Posterior distribution of .

Figure S6. Comparison of the likelihood associated with different values of .

Figure S7. Impact of alternative character coding.

Figure S8. Example of MCMC trace for Carnivora analysis.

Table S1. The count for haploid number includes all chromosomes (sex chromosomes and autosomes).

Table S2. For each parameter the mean of the posterior distribution is given followed by the 95% credible interval (i.e., highest posterior density).

Click here for additional data file.