Mito-nuclear selection induces a trade-off between species ecological dominance and evolutionary lifespan

Mitochondrial and nuclear genomes must be co-adapted to ensure proper cellular respiration and energy production. Mito-nuclear incompatibility reduces individual fitness and induces hybrid infertility, suggesting a possible role in reproductive barriers and speciation. Here we develop a birth-death model for evolution in spatially extended populations under selection for mito-nuclear co-adaptation. Mating is constrained by physical and genetic proximity, and offspring inherit nuclear genomes from both parents, with recombination. The model predicts macroscopic patterns including a community's long-term species diversity, its species abundance distribution, speciation and extinction rates, as well as intra- and inter-specific genetic variation. We explore how these long-term outcomes depend upon the microscopic parameters of reproduction: individual fitness governed by mito-nuclear compatibility, constraints on mating compatibility, and ecological carrying capacity. We find that strong selection for mito-nuclear compatibility reduces the equilibrium number of species after a radiation, increases the species' abundances, while simultaneously increasing both speciation and extinction rates. The negative correlation between species diversity and diversification rates in our model agrees with the broad empirical pattern of lower species diversity and higher speciation/extinction rates in temperate regions, compared to the tropics. We therefore suggest that these empirical patterns may be caused in part by latitudinal variation in metabolic demands, and corresponding variation in selection on mito-nuclear function.


Introduction
The origin, maintenance, and extinction of species involve diverse ecological and evolutionary mechanisms, so that they are rarely explained by a single hypothesis [1,2].
Mathematical models have still been helpful in understanding empirical patterns of biodiversity, such as species abundance distributions [3][4][5] and the species-area curve [6,7], which are well described by both neutral and niche-based theories of community assembly [8,9]. However, models that explain these ecological features in a community seldom make predictions for long-term evolutionary patterns, such as speciation rates [10] and species lifetimes [11][12][13], nor do they provide a mechanistic description for the formation of species [14,15]. Studies describing mechanisms of speciation, on the other hand, using both genetic [16,17] and ecological [18][19][20] drivers, rarely connect these processes to macro-evolutionary patterns [21][22][23][24], whereas macro-ecological models succeed in reproducing such patterns [25,26] but may lack a detailed connection with organism and population properties to explain variation in diversification rates [27,28]. In this context, linking micro-evolutionary processes to macro-evolutionary patterns is key to understanding what drives diversification [29][30][31][32][33], especially when similar evolutionary signatures can result from different population processes [22,23,34].
Here we study spatially distributed populations that evolve under selection for mating compatibility between parents as well as mito-nuclear compatibility within individuals [35].
Genetic interactions between the mitochondrial (mtDNA) and nuclear (nDNA) genomes is necessary for cellular respiration, which depends on biochemically compatible proteins [36,37]. We develop an individual-based model to describe selection on mito-nuclear compatibility and we follow the process of radiation and eventual equilibrium towards a stationary diversity of species. We incorporate the mitochondrial genetic material and its interaction with a portion of the nuclear genome in a phenomenological approach [35] based on the lock-and-key principle that guides the interaction between proteins coded by both genomes [38,50]. Simulations of our model record all forms of diversification events -speciation, extinction, and hybridization -and evaluate the impact of selection on the emergence and persistence of species. In addition, we can use the model to reconstruct the history of an extant community, analyzing the dominance of species throughout their evolutionary lifetime. The model reveals how mito-nuclear selection can jointly influence species abundances and diversity in a community as well as long-term patterns of species lifetimes.
We use the results of our model to help explain longitudinal patterns in species diversity and speciation rates, related to a cline of environmental harshness. Selection on mito-nuclear compatibility is expected to be stronger in harsh, temperate environments that demand optimal metabolic function. Whereas tropical environments, where resources are abundant and energy efficiency is less important, present weaker selection on mito-nuclear compatibility and function.
There are well-known patterns in diversity from tropical to temperate regions, with the former displaying high biodiversity, old species, and low speciation and extinction rates, and the latter with fewer species and elevated recent speciation rates [52][53][54][55][56][57]. We hypothesize that mito-nuclear selection might play a role in forming these patterns, as resources availability should correlate to selection on metabolic efficiency which, in turn, governs patterns of diversity and speciation rates in our model.

Model for mito-nuclear interaction
We employed an agent-based model (Fig. 1, Methods) based on [6], in which neutral processes driven by recombination and mutation change the allele frequencies in a population.
The environment is homogeneous and the population size is kept constant at the carrying capacity. Mating is restricted by both genetic similarity and spatial proximity, so that species emerge through isolation by distance [22]. Here we start from the same neutral assumptions and introduce the mito-nuclear genetic compatibility as a measure of fitness of individuals.
Simulating populations that evolve under these micro-evolutionary dynamics, we catalogue the resulting macro-evolutionary patterns under weak and strong mito-nuclear selection.
Following [35], each individual is described by a haploid, bi-allelic sequence of B sites that represent its nuclear genome (nDNA). To describe the mito-nuclear interaction we also consider a shorter string of B M sites representing the individual's mitochondrial genome (mtDNA). Sexes are separated, and sexual reproduction with non-overlapping generations occurs as follows: during the reproduction phase, all individuals have a probability to randomly choose a mating partner within a spatial range of radius S, called the mating neighborhood. Genetic compatibility between individuals is determined by the Hamming distance between their nDNA's, i.e., the number of loci bearing different alleles, and must be below the threshold of G sites for successful mating. The nuclear genome of the offspring is built from recombination of the nuclear material of the parents, with a probability of mutation µ per site. The mtDNA is directly inherited from the mother, with mutation probability per site µ M . Starting from identical individuals, species emerge as clusters of genetically compatible individuals, i.e., as groups with ongoing gene flow among individuals, determined by the genetic mating restriction G, but no gene flow between different groups (see Methods Section 4.1).
We model mito-nuclear compatibility in an individual as a locus-by-locus interaction between the mtDNA and the first B M sites of the nuclear genome (the remaining (B − B M ) nuclear genes do not participate in the interaction); the interacting pair of loci is compatible if they have the same allele value (Fig. 1b). The mito-nuclear distance d is then the fraction of incompatible sites. Such a scheme is a phenomenological description of the biochemical and structural compatibility between proteins coded by both genomes necessary for the respiration process [37,50]. In the absence of selection, all individuals have an equal probability of reproduction, and nDNA and mtDNA evolve independently. When mito-nuclear compatibility is under selection, fitness w is assigned to individuals accordingly to their mito-nuclear distance, with strength expressed by a parameter σ w , such that w = e −d 2 /2σ 2 w , i.e., small σ w implies strong selection. Individuals with higher fitness, normalized within the mating range, have a higher probability of reproduction. We tracked all events during the dynamics, registering speciation, hybridization, and extinctions, and recorded the ancestry of all species.   Figure 1: Model of selection on mito-nuclear compatibility. a, Individuals are randomly distributed in a square domain and carry nuclear (nDNA) and mitochondrial (mtDNA) genomes described by binary sequences of B and BM sites, respectively. Sexual reproduction occurs with non-overlapping generations. Individuals look for a mating partner within a spatial range of radius S, and they must be genetically compatible (within nuclear distance G) to mate. Offspring are generated with recombination of the nDNA and mutation probability µ per locus, with direct inheritance of maternal mtDNA with mutation µM . b, Mito-nuclear compatibility of an individual is measured by the distance between its mtDNA and BM corresponding within the nDNA. In the absence of mito-nuclear selection, all individuals have equal fitness. Under selection for mito-nuclear compatibility fitness is assigned to an individual according to its mito-nuclear distance, with strength of selection σw, with higher fitness implying higher probability to reproduce. c, Species are classified as groups of individuals with ongoing gene flow, and their abundances are recorded at each generation. We register all speciation, hybridization, and extinctions events, and record the lifespans of all species.

Diversification rates and species richness
We simulated the evolution of populations under various strengths of selection, σ w . In the non-interacting scenario (NI), the mito-nuclear distance is not considered for reproduction (equivalent to σ w → ∞). Selection was simulated with σ w ranging from 0.175 (weak selection) to 0.025 (strong selection). Results are shown for 50 independent realizations for each set of parameters. Starting from a population of identical individuals, all simulations are characterized by a period of rapid increase in species numbers (a radiation period), followed by relaxation to a steady-state species richness. During the dynamics, all diversification events were registered (Fig. 2a). The populations were simulated for T = 2, 000 generations, a time longer than necessary for the completion of the radiation. The equilibrium, defined by the moment when the number of species reaches a plateau, occurred after about T = 500 generations in all cases (Fig.   2b). Because the population size is nearly constant, equilibrium results from a balance between speciation, extinction and hybridization events (compatible with the hypothesis of ecological limits [58]). Selection had the clear effect of reducing the number of species in equilibrium and introducing a delay in the radiation process.  We define a speciation event, or branching, as the moment when a species splits into two or more reproductively isolated groups. The lifetime of the mother species is considered over at this point and the species disappears as two (or more) sister species are born, with a positive net balance in the number of species. Species also disappear by extinction, when the population becomes very small and eventually does not produce descendants due to random ecological drift and by the accumulation of mito-nuclear incompatibilities that lead to low fitness and low fecundity; or by hybridization when two or more species reestablish gene flow. When species hybridize, we consider that the more abundant species persists while the other one disappears as it merges. Linear regression of the number of events in equilibrium provided the rates of speciation, extinction, and hybridization (Fig. 2c). In the equilibrium, when the average number of species becomes constant, the rate at which species are born (speciation) is equal to the rate at which species die (extinction plus hybridization). In this regime, speciation, extinction and hybridization rates per species increase with the strength of selection (Fig. S1). However, since the number of extant species decreases with selection, the nominal rates have a maximum for intermediate values of σ w , as shown in Figure 2c. The decrease in species richness, on the other hand, results from the radiation process, because selection decreases the net diversification (speciation minus extinction and hybridization) during the transient (Fig. S2).

Species abundance distributions and lifespans
Selection on mito-nuclear compatibility decreased the number of species in equilibrium, and because the population size is fixed, this led to an increase in the average species abundance and variance of the species abundance distributions (SADs) (Fig. 3a). Since our simulations resulted in an average number of species no greater than 25, we accumulated data over several independent runs to obtain a clear signal for the expected distribution of abundances, for each value of σ w . This procedure is justified because the lattice is much larger than the mating range and species' range size are smaller than the available space (see Supplementary Information, Section S2).
We recorded species abundances in generation T = 1, 000 and found a roughly lognormal shaped distribution of species abundances, even in the presence of strong selection on mito-nuclear compatibility (Fig. 3a). Selection increased both the mean species abundance and the variance: the histograms became wider and right-shifted under strong selection (Fig. S5).
The SADs are asymmetric, and, although they present a tail on the left side, low-abundant species were infrequent because speciation occurs in our model by fission of abundant species, producing relatively few rare species (similarly to a stick breaking speciation model in a community with finite resources [60]). SADs were stable during equilibrium, such that distributions taken at any time step after T > 500 generations were very similar (Fig. S4).
The time of birth and disappearance of all species were also recorded, then species age and remaining lifespan were known for all species present at time T = 1, 000 (Fig. S6). The distribution of species lifespans with species binned by same abundance class (Fig. 3b) shows that selection tended to reduce species ages and lifetimes in all abundances classes. In all cases, species born with intermediary relative abundance had the highest longevity within their communities. In the absence of selection (NI), species lifetimes were short for low-abundant species and were equally distributed between species of intermediate or large abundance. Some of those species persisted until the end of the simulation (T = 2, 000). Surprisingly, large species evolving under strong selection showed drastically reduced longevity. In the non-interacting case (no selection), the most-right filled bin (log2 abundance class 7) had a slight reduction of lifespan, indicating that there was also a small effect due to species abundance, although selection seems to be the greater cause of lifespan reduction. With regards to species ages at the  diversification events (Fig. 4), i.e., the mode of disappearance of each species (branching, extinction, or hybridization) or whether it persisted through generation T = 2, 000. Regardless of selection strength, most species had two possible endings: by branching, if it was dominant within the community, or by extinction, if it had low abundance. Remarkably, in the absence of mito-nuclear selection, a significant portion of species persisted through to generation T = 2, 000. Selection suppressed the persistence of species, even when weak (σ w = 0.175).
Selection also increased the rate of hybridization (note that hybridization events were counted only for low-abundant species due to our convention of considering the least abundant species involved as dying during the fusion, whereas the more abundant species persists with increased abundance). Notably, species belonging to the same class of abundances had different fates, depending on whether they were under strong, weak, or no selection. For instance, a species belonging to bin 6 in the non-interacting scenario could persist for 1,000 generations ahead, while a species of the same abundance was most likely to go extinct in any scenario with selection. Whereas abundance class 7 was the largest one in simulations without selection, more likely to speciate, but it was not dominant and more likely to go extinct under strong selection.
Greater species abundances are expected when the number of species is smaller, because we assume a constant population size; but it is not clear if the macro-evolutionary patterns we observe are attributable to species abundances per se, or to mito-nuclear selection. To tease apart these mechanisms, we simulated populations in the absence of mito-nuclear selection, relaxing the genetic threshold for mating compatibility (from G = 75 to G = 220) so that they produced a similar number of species in steady-state as communities under strong mito-nuclear selection (Fig. S11). Species-abundance distributions were shifted towards higher abundances as expected, but all diversification rates were reduced, and species lifespans were significantly longer (hence turnover was not increased), even though radiation was delayed (Fig. S10). Therefore, we conclude that the relationship between high speciation rate and low species diversity results directly from the selection on mito-nuclear compatibility.

Genetic diversity of the populations
The analysis of genetic patterns can help to elucidate the mechanism by which mito-nuclear selection promotes or hinders diversification. We evaluated genetic diversity by measuring the nuclear genetic distance between pairs of individuals at time T = 2, 000. Mito-nuclear selection drastically reduced the genetic distances between individuals belonging to different species (Fig.   5a). In the absence of selection (NI), species were genetically well isolated from each other, with a prominent peak around 650 loci. Increasing the strength of selection (reducing σ w ) made species progressively more similar to each other genetically. However, the effect was opposite for conspecific individuals (Fig. 5b): under selective pressure, we observed a slightly larger nuclear divergence within a species than in the non-selective case. Notably for strong selection, several pairs of individuals had distances exceeding the genetic threshold (G = 75) required for mating.
In this case, species cohesion was kept by ongoing gene flow.
The patterns of diversity between different species were similar for the mtDNA content fitness (low mito-nuclear distance), fixation of a mutation in any of the genome is hampered, because it must be followed by a mutation in the corresponding site of the interacting pair to maintain the mito-nuclear compatibility. As a result, interacting sites had lower substitution rates in both the mitochondrial and nuclear genomes (Fig. S8). Within a species, the mean nDNA distance is positively correlated to the species' abundance ( Fig. 5c), independent of the selection strength. As selection increased species abundances, the effect on the population was to increase the nDNA distances within each species. Therefore, the mito-nuclear interaction promoted the intraspecies genetic diversity but reduced the global diversity in the community.
A chain of effects connects selection on mito-nuclear genetic compatibility within an individual to skewed species abundance distributions in a community, shorter species lifespans, and elevated macro-evolutionary turnover (chart at Fig. 6). Strong selection promotes the conservation of the nuclear genetic content in the population because mutations accumulate more slowly. Consequently, the disruption of gene flow is hindered, and the number of species  T = 2, 000 generations. We also show the mean nDNA distance within each species, plotted against the species abundance classes (c). nDNA distances between different species are much lower under strong selection, while the nDNA distances within species are slightly higher. The latter effect is largely due to species abundances: as shown in c, more abundant species allow for greater intraspecific diversity, independent of selection strength. diminished. With a reduced number of species, they are naturally more abundant and support more diversity, favoring evolutionary branching events. Therefore, the speciation rate increases, and abundant species have reduced lifetimes. Low-abundant (non-dominant) species are also short-lived for being more susceptible to fluctuations in the abundance. The reduced genetic distances between different species also promotes hybridization, with sister species merging more easily after speciation. As a global result, communities under mito-nuclear selection were characterized by a low number of species with larger abundances, short species lifespans and high rates of speciation, extinction, and hybridization, i.e., high turnover rate.

Discussion
We have investigated how micro-evolutionary processes impact ecological and macro-evolutionary patterns in a model with selection for mito-nuclear compatibility.
Speciation, hybridization, and extinction are emergent phenomena in our model, byproducts of changes in allele frequencies in the populations, as well as spatial and genetic restrictions on mating. We have shown that strong mito-nuclear selection at the individual level triggers a succession of effects in populations and species, with low species richness in equilibrium, reduced species lifespans, and high evolutionary turnover (Fig. 6). We investigated the mechanisms by which selection induces these effects, and how they are correlated across scales.
The negative correlation between species diversity and speciation rates predicted by our model has been observed empirically -in particular, across latitudinal gradients.  6). These empirical patterns are likely the result of many different ecological, evolutionary, and biogeographical processes [1,25,26], or even neutral processes [62], such that there is no consensus about the causes of the richness gradient [61,63,64] or even the homogeneity of patterns [65] (ecototherms [55] versus planktonic foraminifera [66], for instance, exhibit a positive correlation between richness and speciation rate). We suggest that selection for mito-nuclear compatibility is a possible mechanism that connects environmental harshness, population genetics and species trajectories. In this way, our study helps to link macro-evolutionary patterns of diversification and species turnover, from first principles of mating and fitness operating within a population [27].
Temperate zones are generally regarded as harsh environments with low resource availability [53, 54,67], in which metabolic function is critical [49]. We therefore assume that selection on mito-nuclear co-adaptation is stronger in such zones. The empirical correlation between environmental harshness and the differentiation within species [53,54] suggests that speciation is frequent in those regions, which is supported by the observation of young species ages [68]. Cutter and Gray [67] scrutinized how harsh environments promote the recurrent fragmentation of species in space, which may sequentially develop reproductive barriers [10,23]. Those incipient species evolve through local adaptation and ecological speciation, and are therefore more prone to hybridization. Our model produces similar patterns in a process driven by stabilizing (purifying) selection, in contrast to the mechanism of ecological speciation (based on divergent selection) -which shows that the same macro-evolutionary pattern can result from distinct population processes.
The genetic patterns in our simulated populations reflect the primary causes of variation in diversification rates. A negative correlation between speciation rates and species richness based on a spatially explicit genetic model was also obtained in [24], as a result of sexual reproduction (independent of selection strength). In our model, by contrast, strong mito-nuclear selection reduces the nuclear genetic distance between species (Fig. 5a), as the maintenance of mito-nuclear compatibility slows the fixation of mutations for both nDNA and mtDNA (Fig.   S8). Substitutions occur more slowly, delaying and prolonging radiation and promoting the conservation of genetic content in the population. Since gene flow is higher, there are fewer species that are each more abundant and contain greater intraspecific genetic diversity (Fig.   5b,c). Those abundant species branch more frequently, shortening their lifetimes and resulting in a high speciation rate despite the lower substitution rate of the nuclear alleles. Notably, the average rate of substitutions was not positively correlated with the speciation rate [69], which has previously been proposed to be a driver of biodiversity (revised by [63]). The large intra-species genetic diversity in communities with low richness in our model are similar to wider niches [70,71] and are consistent with predictions of the niche hypothesis (which proposes that ecological opportunity is greater where diversity is low, leading to faster evolution at high latitudes [69]).
Although our model reproduces the broad empirical signatures of mito-nuclear coevolution across the genomes [72,73], it does not produce strong signatures of positive selection associated with the compensatory evolution of nDNA in response to mutations in mtDNA [41,74]. This is because we have simulated an adaptive radiation starting from perfectly coordinated mtDNA and nDNA, so that the majority of selection on interacting mito-nuclear sites is purifying selection.
In contrast to rare species that are more extinction-prone, abundant species are expected to be old and long-lived, i.e. more robust to the ecological drift [4,11]. However, we found that strong mito-nuclear selection reduced species' lifespans even though they were abundant (Fig.   3). We conclude that, for abundant species, the intrinsic robustness to the ecological drift was opposed by purifying selection against mito-nuclear incompatibilities that rapidly decreased individual fitness as mutations accumulated. Frequent speciation by species' fission was responsible for the death of abundant species under strong selection, which explains the relationship between their short lifetime and the corresponding high speciation rate (also suggested by [13]).
Genetic divergence between species is generally greater in the tropics than temperate regions, which is usually explained by the long age of tropical species [55]. Likewise, intraspecific divergences are generally higher in temperate zones, associated with recent speciation rates [53,54,68,75,76]. Our results suggest that the process of accumulating intraspecific genetic diversity may be the cause, rather than the consequence, of a species being older or younger.

Model
We use an individual-based model of reproduction and eventual speciation adapted from [6,22], following the methodology of previous work [35]. We consider a spatially distributed population evolving in a homogeneous environment, in the absence of geographic barriers, and under the influence of neutral micro-evolutionary processes -mutation and genetic drift. Resources are finite but readily available, so that the population is always close to the carrying capacity M [58]. Reproduction occurs through pairwise mating of parents.
Derrida & Peliti [77] first demonstrated this simple model of sympatric speciation with infinite bi-allelic genomes in the absence of natural selection. For sexual reproduction, the only requirement for speciation to occur is genetic compatibility, a weak form of assortative mating requiring a threshold genetic distance G between mating pairs [77,78]. Implementation of finite genomes [35] requires also spatial proximity of mating individuals, which must be separated by a maximum spatial distance S, changing from a sympatric to a parapatric mode of speciation [78].
We call the area of radius S around an individual its mating neighborhood. These are the minimal conditions that guarantee the evolution of reproductive isolation [79], and they allow us to introduce selection on mito-nuclear compatibility into the model. Species emerge from the dynamics as groups of individuals connected by gene flow and separated from all others by the genetic mating restriction G. Not all individuals of a species need to be compatible since gene flow can be established through intermediary individuals. In other words, species correspond to the components of a network where individuals (nodes) are connected if their nuclear genetic distance is less than or equal to G. The model predicts patterns of biodiversity, such as species abundance distributions (SADs) and species-area relationships [6], similar to previous neutral models [3,4].
In order to include mito-nuclear selection, individuals are described by a nuclear DNA, which undergoes recombination and defines genetic compatibility for mating, and by a mitochondrial DNA, inherited from the mother. Both genomes are modeled as binary strings with sizes B and B M respectively. Sexes are separated and individuals are randomly distributed over a square area of size L × L (see Fig. 1a) with periodic boundary conditions (individuals can, by chance, occupy the same lattice site). Generations do not overlap and the population is fully replaced by the offspring. Each individual has a probability P to reproduce and Q = 1 − P of dying without leaving offspring. For reproduction, a mating partner is randomly chosen among genetically compatible individuals in its mating neighborhood. The nuclear genetic distance D n between individuals is defined as the Hamming distance between their nuclear genomes, and compatibility requires D n ≤ G. If a genetically compatible mate is found within the mating neighborhood of the 'focal' individual, an offspring is produced with locus by locus recombination of the nuclear genomes, followed by mutation with probability µ per locus. The mtDNA is inherited from the mother, with mutations with probability µ M (Fig.   1a) We compare the neutral process, when mito-nuclear selection is absent, with the scenario where mito-nuclear compatibility promotes fitness differences in the population, testing different levels of selective pressure while keeping all other parameters fixed. The mito-nuclear distance d measures the compatibility between the mtDNA and the corresponding loci of nDNA, assigned to each individual, ranging from 0 (totally unmatched sequences) to 1 (fully compatible sequences); for completely random sequences, the expected value of d is 0.5. The individual's fitness w depends on d according to a Gaussian function with width σ w , which quantifies the selective strength of the mito-nuclear interaction (Fig. 1b). The parameter σ w is related to the population average d in equilibrium (see Fig. S9), such as the absence of selection leads to d ≈ 0.5 and increasing the selection diminishes d towards 0. Selection over mito-nuclear compatibility changes the probability of an individual to reproduce, i.e., it affects its fecundity. The focal individual has its chance for reproduction modified by its fitness as follows: first the fitness of all individuals in the population are computed as w(i) = e −d(i) 2 /2σ 2 w ; second, Q for individual i is modified to Q w (i) = 2Q(w max − w(i))/(w max − w min ) where w max and w min are the maximum and minimum fitness in the population, respectively. Individuals with w(i) = w max have Q w (i) = 0 and probability of reproducing P w = 1 − Q w = 1. Those with w(i) = w min have Q w (i) = 2Q and P w = 1−2Q. Therefore, individuals with low fitness still have a small chance of reproducing.
If the individual dies without reproducing, another one is selected in its mating neighborhood, also according to fitness, to reproduce in its place, keeping the total population constant. Finally, a mating partner is selected from the compatible individuals in the mating neighborhood, with probability proportional to their fitness.  Branching events increase species richness, while extinction and hybridization decrease the number of extant species, reducing biodiversity. Hybridization is an emergent event in our model and only cause loss of diversity. In nature, however, hybridization may increase species richness inducing the evolution of reproductive barriers [80]. Here we emulated the loss of diversity due to hybridization from lineage fusion, speciation reversal [81], and hybrid breakdown [40]. We included in this category events of failed speciation attempts (when two or more recent sister species merge back) and events of reticulate evolution (when there are two speciation events in a row and a species fuses with the one with the most distant ancestor)

Species abundances distributions
In our simulations the number of species in equilibrium (T ≥ 500) ranged from about 20 in the non-interacting case to about 5 for strong selection (Fig. 2b). To generate abundance distribution plots with statistical significance, we have run 50 independent simulations for each set of parameters and accumulated the results. This increased the number of species in the histograms to about 1,000 for the non-interacting case and 100 for the case of strong selection.
Section S2 of Supplementary Information justifies the validity of this procedure.
To generate the species' abundances histograms, data were binned in abundances classes (log2). There are several methods for choosing the bins' boundaries, and we followed the method adopted by [3]: bins were built by using the powers of two as the center of bins. The boundaries then come as 2 n−1/2 and 2 n+1/2 . Thus, bin=0 stores species with 1 individual, bin=1 counts species with 2 individuals, bin=2 species with 3, 4, and 5 individuals, and so on.

Data availability
There are no empirical data associated with this study.

Code availability
All simulations are run in Fortran. All the codes for simulations and scripts for data treatment

S1 Diversification rates
In this section we show how the diversification rates -speciation, hybridization and extinction -depend on the selection strength at equilibrium and how they behave along the evolutionary process (radiation and equilibrium). Figure S1 shows the diversification rates per species in the equilibrium regime for different selection strengths σ w and in the absence of selection (non-interacting case, N I) (analogous to the rate of speciation per lineage, as usually shown for comparison of global patterns [1]). We note that the number of species in equilibrium varies considerably with σ w (Fig. 2b of the main text). Figure S2 show the diversification rates along the evolutionary time (instantaneous rates evaluated for time bins ∆t = 10 generations), displaying peaks during the radiation and reaching a steady state afterwards. The net diversification rate, result from speciation minus extinction and hybridization, peaks at the radiation and it is null in the equilibrium, when species richness reaches a plateau. Ns is the total number of species, Mt is the total number of individuals andms is the mean species' abundance over all simulations.

S3 SADs at intermediary times
Here we show that, after equilibration at T ≈ 500 generations, the abundances distributions do not change significantly, reaching a steady configuration. Figure S4 shows that, despite small fluctuations due to the stochastic nature of the model, the overall shape of the distributions do not change.  Figure S4: Abundance distributions computed at three different times after the radiation period.Ns is the mean number of species andms is the mean species' abundance. Although small fluctuations are observed, the overall shape does not change significantly over time.

S4 Mean, variance and skewness of the SADs
The differences between the histograms of Figure 3a (Species abundances distributions at T = 1000) can be quantified simply by the first, second, and third momentum of the distribution, i.e., the mean value, the variance, and the skewness of the histograms (Fig. S5). We considered the center points of the bar diagrams where the abundances were binned in log2 abundances classes, as described in Section 4.3. Therefore, we did not assume any best fit to the data (lognormal, Fisher distribution, etc.) because we aimed only to see the effect of selection on the distribution. Previous work on this model discussed the resulting patterns of species abundances (in the absence of mito-nuclear selection) and compared them with empiric data [2].
Selection had the role of increase the mean size of species (Fig. S5, top) and the distributions also became wider (Fig. S5, middle). However, the skewness did not preset a monotonic tendency with increasing selection (Fig. S5, bottom), all distributions were left-skewed (skewness with negative value) and the non-interacting model is the most skewed.

S5 Species' ages and time for death
Here we investigate how species abundances correlate with age at the time of sampling and remaining lifetime. Species were detected at T = 1, 000 and their abundances and ages computed. Figure S6 shows that species of intermediate sizes are older (Fig. S6a) and live longer ( Fig. S6b) when selection over the mito-nuclear interaction is considered. The non-interacting case shows no significant correlation between lifespans and abundances, except for small species, that tend to go extinct faster.  Figure S7 shows the distribution of mitochondrial distances between (Fig. S7a) and within ( Fig. S7b) species. Selection over mito-nuclear interaction drastically reduces the interspecies distances, implying in more genetically similar species. Intraspecies distances, on the other hand, change from an exponential distribution with large average distance (non-interacting case) to a unimodal distribution with much smaller average distance. Species, therefore, become more uniform in terms of their mitochondria as selection strength increases. For the mitochondrial DNA, which does not recombine and is not under mating restriction, selection utterly increased similarity within species, spontaneously forming clusters of mitochondria (Fig. S7b).

S7 Substitutions rates
In our simulations the mutation probabilities per locus were fixed at µ = 0.00025 for the nDNA and µ M = 0.00075 for the mtDNA. Mutation reversals, however, decrease the effective number of mutations and selection filters out 'deleterious' mutations that decrease fitness, also decreasing the net rate of substitutions.
The number of substitutions is computed as the number of loci that changed from the initial state 0 to state 1 at T = 2, 000. For the non-interacting case, the probability that a single mutation occurs and is not reversed after T generations is T 1 µ(1 − µ) T −1 . The probability that a mutation occurs and is reversed twice is T is an estimate for the average number of substitutions per generation. The expected substitution rates of the nuclear and mitochondrial genomes in the non-interacting case are then 1.5 × 10 −4 and 2.3 × 10 −4 per generation, respectively [3]. Figure S8a shows the substitution rate per locus in the mtDNA in the non-interacting case (red) and the strongly interacting case (blue). Figure S8b shows similar information for the nDNA, separated into coupled sites (first 500 loci) and uncoupled sites (remaining 1,000 loci). The resulting population-averaged mito-nuclear distance, computing the Hamming distance between coupled sites in the nDNA and the respective sites in the mtDNA, is shown in Figure S9. In the absence of selection, the sequences evolve independently, and d reaches the maximum value of 0.5 (expected value for completely random sequences). When selection is present, the steady-state value of d is smaller with stronger the selection, i.e., the mito-nuclear compatibility is preserved.

S8 Non-interacting case with reduced number of species
In the main text, we compared the evolution of populations under variable selection strength while fixing all other parameters, which resulted in a decrease in the number of species as the intensity of selection increases (Fig. 2b). Here we compare the case of strongest selection (σ w = 0.025) with a neutral case presenting similar number of species. The purpose of such comparison is to evaluate whether the effects on diversification rates and species' lifespans and fates are simply a consequence of the decrease in species richness. The non-interacting case (N I) presented a similar number of species of σ w = 0.025 when the mating restriction was decreased, setting G = 220 (G = 0.15B, in contrast with G = 75 = 0.05B in the main text), as shown in Figure S10b. The radiation process was delayed, and the number of diversification events drastically decreased (Fig. S10a), also reducing the rates at the steady-state (Fig. S10c). Due to the reduced number of species, they became more abundant, right-shifting the SADs ( the distribution was not much wider (Fig. S12 middle). More important, species' lifespans were not reduced: comparing N I (G = 220) and σ w = 0.025, species in the same abundance class lived longer in the absence of selection than under strong selection (Fig. S11b). The remaining lifetimes were also longer (Fig. S13b), and, even with the delay in the radiation, species at T = 1000 were older without selection (Fig. S13a). Therefore, the reported reduction of lifespan and diversification rates were due to the selection over mito-nuclear compatibility.    These evolutionary patterns connect once again to the genetic structure of the population: even though the increase in the genetic threshold G made nDNA distances within species larger ( Fig. S15b), the distance between different species did not change, compared to N I (G = 75) (Fig. S15a). We observe that: larger nDNA distances within species hinder speciation (N I (G = 220)), while shorter distances between different species promotes speciation and hybridization (σ w = 0.025). Mitochondrial distances for N I (G = 220) were similar to N I (G = 75), supporting more diversity intraspecies.