Species‐range‐size distributions: Integrating the effects of speciation, transformation, and extinction

Abstract The species‐range size distribution is a product of speciation, transformation of range‐sizes, and extinction. Previous empirical studies showed that it has a left‐skewed lognormal‐like distribution. We developed a new mathematical framework to study species‐range‐size distributions, one in which allopatric speciation, transformation of range size, and the extinction process are explicitly integrated. The approach, which we call the gain‐loss‐allopatric speciation model, allows us to explore the effects of various speciation scenarios. Our model captures key dynamics thought to lead to known range‐size distributions. We also fitted the model to empirical range‐size distributions of birds, mammals, and beetles. Since geographic range dynamics are linked to speciation and extinction, our model provides predictions for the dynamics of species richness. When a species‐range‐size distribution initially evolves away from the range sizes at which the likelihood of speciation is low, it tends to cause diversification slowdown even in the absence of (bio)diversity dependence in speciation rate. Using the mathematical model developed here, we give a potential explanation for how observed range‐size distributions emerge from range‐size dynamics. Although the framework presented is minimalistic, it provides a starting point for examining hypotheses based on more complex mechanisms.


| INTRODUC TI ON
The distribution of species-range sizes is a key pattern of interest in macroecology, providing a basic spatial profile of the earth's biodiversity (Gaston, 1998). Despite the broad range of biological dynamics that potentially influence such a macroscopic pattern, a relatively consistent form is often observed across taxa and regions: Speciesrange sizes have a lognormal distribution with left skew (Brown et al., 1996;Gaston, 1998;Gaston & Blackburn, 1996Gaston & He, 2002;Noonan, 1990;Orme et al., 2006;Ruggiero, 1994). Gaston (1998) provides skewness values from a wide range of taxonomic groups where 19 out of 22 range-size distributions have a skewness between −1.26 and −0.03 on a log-axis (note right-skewed distributions were also reported but these were minor).
This generality in the form of the pattern raises the possibility of some similarly general mechanism. While there is a long list of processes that potentially influence range-size distributions, the shapes of these distributions are ultimately a product of speciation, extinction, and transformation of the range (Gaston, 1998). Therefore, exploring these patterns leads to a central question in ecology and biogeography: How do the processes of speciation, extinction, and transformation mutually interact to shape range-size distributions? (Brown et al., 1996;Gaston, 1998).
Given the large spatial and temporal scales involved, quantitative approaches have been central to the study of species ranges. These include studies with a particular focus on inferring speciation modes (Cardillo, 2015;Phillimore & Price, 2008;Skeels & Cardillo, 2019), phylogenesis (Albert et al., 2017;Pigot et al., 2010), niche evolution (Rangel et al., 2007(Rangel et al., , 2018, size-age relationship of geographic ranges (Pigot et al., 2012), and the heritability of range size (Borregaard et al., 2012). In addition to the focus on individual species-range sizes, previous studies also shed light on the shape and the underlying mechanisms of the emergent distributions of range size across many species (e.g., Alzate et al., 2019;Anderson, 1985;Gaston & He, 2002;Pigot et al., 2010;Rangel et al., 2007). In early work, Anderson (1985) focused on range-size distributions using an algorithm with nine different scenarios for range dynamics (speciation is associated with an extinction event in the algorithms) to generate range-size distributions resembling faunal data for North American vertebrates. Later, Gaston & He (2002) developed a stochastic model that describes range-size dynamics driven by population dynamics of a single species, and they showed that the generalized equilibrium range-size distributions produced by the model fit well to range-size distributions of several taxonomic assemblages. Pigot et al. (2010) built a spatially explicit numerical algorithm to simulate geographic range evolution where new species arise via vicariance or peripatry, and they reproduced empirically observed range-size distributions across bird genera that lead to diversification slowdown in the reconstructed phylogenetic tree.
While these efforts have been illuminating, and detailed spatially explicit algorithms are available (e.g., Rangel et al., 2018), we still lack a comprehensive understanding of the assumptions necessary and sufficient to realize observed patterns of range-size distributions.
For example, we do not have concrete insight into the mutual importance of key processes, range transformation, extinction, and speciation (Gaston, 1998). Also, there are contrasting views regarding the effect of range size on speciation rate (Gaston, 1998;Jablonski & Roy, 2003;Rosenzweig, 1995;Tokeshi, 1996), and the effect of speciation rates combined with other key processes on emergent range-size distributions has not yet been well explored.
We develop a new class of macroecological model describing key processes of range-size dynamics that we call the gain-lossallopatric speciation (GLAS) model. The GLAS model accounts for the key processes shaping ranges: the gain in area through dispersal, the loss of area due to local extinction, and the splitting of a single range into multiple ranges through allopatric speciation. Our mathematical framework is simple and flexible to accommodate various assumptions, and we examine several scenarios of allopatric speciation. While a priori we could expect any number of range-size distributions to be realized, we ask whether the fact that ranges are subject to gain, loss, and speciation dynamics places any constraints on their distributions. While in principle such a modeling approach could become very complex when incorporating a broad range of pertinent processes occurring in a region (e.g., Rangel et al., 2018), we take a minimalist approach to address general cases. We do not deal explicitly with species interactions (as, e.g., Gaston & He, 2002;Pigot et al., 2010), the genetic mechanisms of speciation, environmental heterogeneity, or change, or any number of other potential elaborations. Rather, we model range-size distributions in terms of phenomenological rates of range growth, contraction, and speciation, which are implicitly affected by all those mechanisms.
The dynamics of species-range-size distributions can also affect the diversification process and structure of the associated phylogeny (Albert et al., 2017;Pigot et al., 2010). If speciation and extinction are linked to range size, then lineage diversification rates can be affected by dynamical changes in the range-size distribution itself. This can occur by concentrating species at sizes with lower or higher speciation and extinction rates. Hence, investigating this theoretically provides an opportunity to link range dynamics to macroevolutionary processes. The dynamics of species diversity are also a central focus in macroecological studies (Nee, 2006;Nee et al., 1992;Rabosky, 2009). Diversification slowdowns, often discussed via lineage-through-time plots (Harvey et al., 1994;Nee et al., 1995), are a common feature of inferred lineage dynamics (Cusimano & Renner, 2010;Moen & Morlon, 2014;Morlon et al., 2010;Nee, 2006;Yoder et al., 2010), yet biological explanations for the underlying mechanisms (e.g., diversity-dependent speciation) are still under active investigation (Condamine et al., 2019;Moen & Morlon, 2014;Pannetier et al., 2021). We demonstrate that the GLAS model provides an opportunity to examine the effect of range-size dynamics on the diversification rate.
We find that left-skewed distributions on a log-axis are predicted under various speciation scenarios and parameter sets. The leftskewed distributions are realized by the similar rates of the expansion and contraction of geographic ranges with the moderate supply of new species. We demonstrate that these three parameters regulate the shape of emergent range-size distributions. Interestingly, our model can generate diversification slowdown in the dynamics of the number of species under all speciation scenarios investigated, even though the model has neither a diversity dependence effect nor a species interactions.

| Gain-loss-allopatric speciation model
Here, we give an overview of the gain-loss-allopatric speciation (GLAS) model. See Appendix 1 for complete discussion of the model. The GLAS model contains minimal but essential factors to describe range-size dynamics and to recover empirical range-size distribution: allopatric speciation, extinction, and transformation of the range size (Gaston, 1998).
Let us assume that a species has a geographic range size r at time t, and experiences gain in its range size by one unit (r → r + 1) or loss by one unit (r → r − 1) at rates g and l, respectively. In addition, allopatric speciation can occur at a size-dependent rate a r , causing the subdivision of a geographic range to produce two smaller range sizes r 1 and r 2 with one new species (Figure 1a). We assume for simplicity that the sum of the two subdivided range sizes is equivalent to the parent range size r 1 + r 2 = r, and all combinations of split size are equally likely. We define an extinction as the event that a species reaches range size 0, and it only occurs from r = 1 (Figure 1b). We assume that each species-range size changes independently of other species. These assumptions are translated into the following stochastic process that describes the dynamics of P(r, t), the expected number of species with range size r at time t: Equation (1) describes the dynamics of geographic range-size distribution over time. It is a class of continuous-time random walk (Gardiner, 2009;Karlin & Taylor, 2012) with an additional term for the effect of allopatric speciation.
It is worth noting that we assume that individual processes within each geographic range induce change of geographic range by "one unit," and hence, the size of change is a constant regardless of the geographic range size. However, it may also be sensible to assume that environmental factors cause the size shift, and lead to a size-dependent rate of change. We will discuss this effect in the Discussion and provide an example.
Because the GLAS model describes speciation and extinction events, it also generates dynamics of the total number of species, an exponential increase, decrease, or stable species number depending on the value of the dominant eigenvalue (see Appendix 1 for details and an example of eigenvalues provided in Figure A1).

| Continuum limit of the GLAS model
It is reasonable to consider that the geographic range size changes continuously rather than in a discrete manner. When the step size of the range size is small (r ∼ r + 1) in Equation (1), we introduce a new variable u(x, t) where the number of species with range size between x and x + Δx at time t is described by u(x, t)Δx. Then, we obtain the following integro-differential equation that describes the F I G U R E 1 Schematic diagram of the model. Each species has a geographic range size at given time t that is labeled by a number. (a) Within a small time period Δt, a species range can increase in size, decrease in size, or undergo allopatric speciation leading to two smaller range sizes with one new species. (b) Extinction can occur if a species has the smallest allowable range size (labeled 1 in this example) and then decreases in size. (c) Four scenarios for the dependence of speciation rate on normalized range size (0 ≤ x ≤ 1). The probability distribution is described by a beta distribution and the parameter values for each scenario are Linear increase: = 2, = 1 ; Linear decrease: = 1, = 2; Parabola: = 2, = 2; and Constant: = 1, = 1. See the main text for more details where V = (g − l)Δx and D = (g + l)Δx 2 ∕2. Note the assumption that the ancestral range is divided into two ranges with an equal probability of any split size, and so yields the average split proportion 25:75. Possibility of split asymmetry in marine mollusks was suggested in Pigot et al. (2012). Also, Budding speciation, by which a reproductively isolated smaller-ranged species is originated from a larger-ranged species by a highly asymmetric fashion, may be common in plant species (Anacker & Strauss, 2014;Grossenbacher et al., 2014), and our model can capture this speciation mode as this is a potentially (spatially overlapped) subdivision of a geographic range.
Here, without loss of generality, we set the range size to lie in the interval 0 ≤ x ≤ 1, with 1 being the size of the entire domain (the maximum possible range, e.g., the size of a continent). Assuming that an extinction event occurs only when a range size reaches x = 0 and no species exits the domain from the upper bound x = 1 (i.e., no-flux condition), we need mixed boundary conditions.
Namely, a Dirichlet boundary condition for the lower boundary u(0, t) = 0 and the Robin boundary condition for the upper bound- We can derive an analytical solution of Equation (2) when speciation increases linearly with range size by imposing extra conditions. However, the analytical form is rather complicated and computationally expensive due to an infinite summation. For the reader's convenience, we provide the derivation in Appendix 2. In the following, we will perform numerical analysis of Equation (2). Also, we will focus on deterministic characteristics of the GLAS model to directly discuss mutual importance of speciation, transformation, and extinction events, while Equation (1) allows us to investigate stochastic processes.

| Scenarios of allopatric speciation
To complete the model Equations (1) and (2), we need to define the size dependency of speciation rate a x . However, as discussed in the Introduction, this is rather contentious. We assume four potential scenarios to capture some of the existing hypotheses, where the probability of speciation: (i) increases with range size (Rosenzweig, 1995;Tokeshi, 1996); (ii) decreases with range size (Jablonski & Roy, 2003); (iii) peaks at an intermediate size, which is attributed to a hypothesis that a high dispersal ability can lead to large ranges and inhibit speciation (e.g., Claramunt et al., 2012;Gaston, 1998;Mayr, 1963); or (iv) is independent of range size. We examine how these four speciation scenarios affect the form of the range-size distribution. We use the beta distribution to model these four probability distributions with a single probability distribution function (PDF).

| Range-size distributions
Having completed Equation (2)  Although we use a normalized model to investigate range-size distributions across different taxonomic groups, this approach does not lose any qualitative property of the original model. The same analysis could be applied for any specific case with empirically estimated parameter values. Within these settings, we examine multiple sets of gain rates g, including the case g > l, g = l, and g < l, and multiple orders of the underlying speciation rates including the case where there is no speciation a = 0. We also perform model fitting to data for mammals and birds in the Americas, and Harpalus carabids in North America north of Mexico. Since empirically observed range-size distributions often span multiple orders of magnitude, we employ the logarithmically transformed GLAS model (Equation A13) to compute the range-size distribution in the analysis where all the characteristics remain unchanged (see Appendix for details).

| Species number
Since the GLAS model characterizes speciation and extinction events, it produces the dynamics of total number of species N(t) = ∫ 1 0 u(x, t) dx where the number of species is also a continuous value as u(x, t) is a continuous value. Given such dynamics, it is possible to calculate diversification rate (t) over time, the difference between the speciation rate at time t and extinction rate at time t. Since the diversification rate has the following relationship (t) = dlogN(t)∕dt and the extinction rate is calculated by , it is possible to compute all these time-dependent rates along with the lineage-through-time plot. Along with these continuous representations, we also provide the statistic r (Etienne & Rosindell, 2012;Pigot et al., 2010) where r < 0 and r > 0 indicate a slowdown and speedup of the diversification rate, respectively, over the period of calculations.

| Species-range data
Species-range data for birds (n = 365) and mammals (n = 628) in  Noonan (1990) where the range size was calculated based on occurrence data.
We used 17,683,892 km 2 to normalize the range sizes (i.e., to organize all the range sizes between 0 and 1). Note the choice of the size for the normalization does not affect the shape of range-size distributions (e.g., variance, skewness, and kurtosis), but it merely causes a shift on the x-axis.

| Model fitting
Although we do not have an explicit form of the likelihood function with arbitrary parameters for speciation scenarios α and β of Equation (2), we can perform data fitting via simulated annealing, an optimization algorithm, based on a log-likelihood using a numerically calculated range-size distribution. Detailed settings of the optimization algorithm are found in Appendix 3. In this process, we set the range of speciation parameters α, β ∈ [1,50]. We sampled the range-size distribution when the changes in skewness and kurtosis between time steps become smaller than 10 −5 (i.e., convergence to equilibrium) and eliminated parameter sets which cause species number smaller than 0.9 until the convergence. Note the value to judge the numerical convergence is different from the value used in other numerical simulations to reduce the simulation time of the optimization algorithm. Since we are interested in relative significance of gain, loss, and speciation, we again set the loss rate to 1.

| Model predictions
Given a speciation scenario, and the three factors affecting rangesize dynamics, gain rate g, loss rate l, and an underlying allopatric speciation rate a, we can numerically solve the GLAS model. We observed that the range-size distribution approaches a unique equilibrium distribution regardless of initial conditions after a certain simulation period. We sampled the range-size distribution when the absolute difference in skewness (kurtosis) at two arbitrary time points (we sampled every 1000 time steps for numerical convenience) became smaller than 10 −7 . Figure 1a shows example dynamics of the range-size distribution under the scenario of linearly increasing probability of speciation with range size and underlying speciation rate a = 0.1. Starting with a single species with a range size 0.1, the range-size distribution converges to an equilibrium distribution that is comparable to some existing data: a left-skewed curve on a log-scale (Brown et al., 1996;Gaston, 1998;Gaston & Blackburn, 1996Gaston & He, 2002;Noonan, 1990;Orme et al., 2006;Ruggiero, 1994).  (Table 1). However, a smaller underlying speciation rate causes a right shift of the distribution, since it suppresses the chance of allopatric speciation, which is the mechanism that pulls the distribution toward the left. As a result, the distributions of the linear decrease and constant speciation scenarios come closer to the other scenarios. We provide heat maps of skewness and kurtosis of the species-range-size distributions across parameters in Figures A2 and A3 showing the robustness of these summary parameters provided by the left-skewed distribution.
If the underlying speciation rate is zero or sufficiently small, the mechanism to reduce the range size by supplying a new species is suppressed and the number of species declines over time (if an initial population size is larger than 1). As a result, the left skew vanishes and many species arrive at the largest possible range size when g ≥ l ( Figure A4) and this corresponds to the left-top region of each panel in Figures A2 and A3. Figures A2 and A3 also show that a similar result occurs when the gain rate is substantially larger than the loss rate and the underlying speciation rate, since this also suppresses the mechanism to reduce the range sizes.
The left-skewed range-size distribution is widely observed across the parameter space, including the case of nonlinear speciation scenarios ( Figure A5), as long as there is a mechanism that avoids many species growing to the maximum possible range size. Although the gain rate tends to be larger than the loss rate (g > l) to have a positive species growth rate, it can still show a leftskewed distribution even when the species growth rate is negative ( Figure A6) and when speciation rate is zero and with a smaller gain rate than the loss rate (g < l; Figure A4) in the course of all species extinction.
The GLAS model characterizes speciation and extinction events, and it produces the dynamics of total number of species. Starting with a single species, we found the diversification rate can decrease, without a diversity-dependent effect, as it approaches the equilibrium diversification rate ( Figure 3). Typically, diversification slowdown is associated with decline of the realized speciation rate and increase in the extinction rate as it approaches an equilibrium.
When the speciation rate increases over time, slowdown is less likely to occur but it is still possible if the extinction rate also increases fast enough to suppress the effect of speciation as in Figure   A8f. Speciation scenarios also affect these dynamics. The linear increase scenario tends to require larger initial range size to show the slowdown, while the linear decrease still can show the strong slowdown with a smaller range size. This trend can be explained by the match between the initial range size and the possibility of speciation. Namely, if an initial range has larger likelihood for speciation, it can produce a higher accumulation of species at an initial phase than the range size with lower likelihood of speciation. We provide these numerical results with two other initial range sizes (0.05 and 0.5) for each combination of the three underlying speciation rates For example, the slowdown can occur when both speciation and extinction show a decreasing trend ( Figure A11a,c) that is not observed with linear speciation scenarios. Also, nonlinear increasing and decreasing speciation scenarios generally induce more significant differences compared to these linear counterparts.  Table 2. We found that our model tends to show narrower equilibrium distributions (spanning around 3 orders of magnitude) than the empirical datasets (5 or more orders of magnitudes). The speciation scenarios selected are nonlinear-decreasing functions, but the parameter values are close to the domain boundaries of parameter space defined. Widening parameter spaces for speciation scenarios (i.e., α and β) allows more extreme, and perhaps unrealistic, speciation scenarios and may not improve prediction.

| Model fits to empirical data
This casts a limitation of the model fit to a continental scale dataset with wide taxonomic scales. However, this may not be surprising as the model is minimalistic and assuming homogeneous biological parameters between species. In fact, our model leads to a better fit F I G U R E 2 Normalized range-size distributions that typically show a lognormal feature but are left skewed. (a) The time evolution of a range-size distribution starting with a single species with range size 0.1 when the speciation scenario is linear increase. The rangesize distribution is recorded every 25 time units for visualization until t = 400. (b) Equilibrium range-size distributions with four speciation scenarios. Other parameter values used are l = 1, g = 1.2 for (a) and (b), and g = 1.1 for (c) and (d). Also, parameters for the speciation scenarios (α, β) are (2, 1) (linear increase); (1, 2) (linear increase); (2, 2) (parabola); and (1, 1) (constant) TA B L E 1 Summary of (skewness, kurtosis) from each curve in Figure 2 to the carabid dataset, which shows narrower distribution due to narrower taxonomic scales in the dataset. Also, the limitation in the  Table 2 taxonomic groups. Our model provides a better understanding at the level of such a component.

| Left-skewed range-size distributions
The gain-loss-allopatric speciation (GLAS) model is a general mathematical framework to model the key factors of transformation, speciation, and extinction processes in the dynamics of range-size distributions.
While the components of the model are ultimately phenomenological and not based on modeling the mechanisms driving those processes, it allows us to investigate how their relative rates and functional forms affect range-size distributions. We found a left-skewed and lognormallike distribution, which has been identified as one of the general patterns in previous studies (Brown et al., 1996;Gaston, 1998;Gaston & Blackburn, 2008;Gaston & He, 2002;Noonan, 1990;Orme et al., 2006), and equivalently a positive skew on a normal axis ( Figures A5 and A6), yet these parameters and speciation scenarios affect the shape statistics of the range-size distributions (Table 1). Typically, this left-skewed pattern is observed when the rates of gain and loss of range size are similar in magnitude, and allopatric speciation rates are also of comparable order to these rates. Balanced gain and loss rates avoid a disproportional chance of extinction and an arrival at the large range-size limit. Allopatric speciation also reduces the chance of a species arriving at the large limit of range size. It induces a split of a range size into two smaller range sizes, pushing the range-size distribution toward the left side, and we attribute the left-skew distribution to this mechanism.
A lower rate of allopatric speciation allows many species to occupy a larger range ( Figure A4). It suggests that the knowledge of the rangesize distribution could be used to diagnose diversification strength relative to the rate of the range-size transformation. Also, a comparable order of magnitude in the rate of allopatric speciation to the transformation rate is necessary to maintain the diversity. For example, a much larger speciation rate than the transformation rates causes increasingly smaller ranges and it can lead to a larger extinction rate than the speciation rate. On the other hand, biodiversity cannot be maintained with a too small speciation rate due to the lack of a source of new species.

| Fitting to range-size distribution datasets
Our attempts of model fitting produced narrower range-size distributions than the empirical datasets for mammals, birds, and Harpalus carabids. The datasets of Harpalus carabids show narrower distributions than the others, and consequently, it leads to a relatively better fit. We hypothesize potential explanations of this variable fitting performance.
First, we assume that all species have the same rate of changes in geographic range size and the speciation rate, and the same speciation scenario. This is a simplification, and wider range-size distributions of birds and mammals may result in larger heterogeneity in these rates. The number of species is also larger for these datasets than the datasets of Harpalus carabids, and these may include more diverse functions, phenotype, evolutionary strategies, and so on, and it leads to more variable parameter values.
Second, we mainly focused on an equilibrium range-size distribution. However, its transient range size does not yield a unique distribution pattern starting with an arbitrary initial range-size distribution. In reality, each rate may fluctuate over time and the observed pattern may not be an equilibrium distribution.
Third, environmental factors, for example, occupancy of glaciers (Noonan, 1990), drive transformation of the range sizes. While we assume that the change in range size is based on individual processes where transformation of range size is by one unit, it may not be the case for environmental change where its increments/decrements are proportional to the range size. We provide an example with an assumption of proportional change in range size that leads to wider range-size distributions in Figure A12 in the Appendix and some technical details below.

| Diversification slowdown
We show that the GLAS model also produces the signature of diversification slowdown in the phylogeny as in Pigot et al. (2010). In

the GLAS model, diversification slowdown occurs during a transient
phase as the diversification rate approaches a stable value. Several hypotheses have been discussed to explain these phenomena including density-dependent speciation (Phillimore & Price, 2008;Rabosky & Lovette, 2008;Weir, 2006), failure to adapt to a changing environment (Quental & Marshall, 2013), and protracted speciation (Etienne & Rosindell, 2012). In our simulations, this slowdown behavior is observed under a wide range of parameter sets and speciation scenarios when the initial range size is closer to a size with high speciation probability. Since speciation tends to move species ranges away from the sizes with highest speciation probability, it tends to push the distribution away as well, slowing down the rate of speciation.
For example, when the original species speciates under the linear increase speciation scenario, the two new smaller range sizes may still have a high chance of producing new species before going extinct if

TA B L E 2
Estimated parameters by the model fitting in Figure 4 the initial species had a relatively large range size. Extinction events are more likely to occur after this transition when there are more species with smaller range sizes, and this eventually suppresses the diversification rate. Although this slowdown is pervasive, opposite trends in the diversification rate can occur (e.g., Figure A7(i) and (k) in Appendix). For example, when the initial species has relatively small range size and speciation probability increases with the range size, there is an initial extinction-prone period since smaller range sizes are less likely to speciate. After the initial transition period, the speciation rate increases due to the existence of large range sizes, and this leads to an upward curve in the lineage-through-time plot. This opposite trend was reported in previous studies and was explained by, for example, climatic cycles and disruptive mountainous terrain where rapid recent speciation might have been induced (Linder et al., 2003;Weir, 2006).

| Mechanisms to extend the GLAS model
While our model is a simple approach to modeling the three key factors as discussed in (Gaston, 1998), the results suggest that the model produces widely observed range-size distributions and diversification slowdown. Thus, it allows us to obtain general insight with minimal assumptions. However, our model can further be extended to investigate a wide range of assumptions, scenarios, and speciation modes. It is possible to assume that the range size x in Equation (2) is an arbitrary parameter to characterize a geometry of the range size (e.g., length and perimeter), not necessarily area itself. If the area is a nonlinear function of x, then the realized range-size distribution may have a wider distribution. For example, consider circular ranges, assuming that Equation (2) describes the dynamics of the radius (x) of the range, rather than range size ( x 2 ) itself, Figure A12 shows a wider range-size distribution under this assumption. This simple example induces the loss of area of the original range size after speciation: If a range size with a radius x splits into two ranges with radius x 1 and x 2 , then we have x 2 ≥ (x 2 1 + x 2 2 ). While we do not have a clear idea if geographic range size is decreased after allopatric speciation, we can assume various geometries that induces different amounts of area loses after allopatric speciation.
Split asymmetry of ancestral species influences the range-size evolution (Pigot et al., 2012), and a large asymmetry may also cause a wider range-size distribution. While our assumption of the equal split probability gives the asymmetric split proportion 25:75, this can be controlled: This effect is incorporated into Equation (2) as where a 2 (x) is any symmetric probability distribution function in 0 ≤ x ≤ x 1 with original range size x 1 , regulating the degree of the split asymmetry. The average split proportion In our analysis, we used a 2 (x) = 1∕x 1 .
We can incorporate range-size-dependent gain and loss rates by making them a function of range size g(x) and l(x). It would be possible to relate to population dynamics where, as in Gaston and He (2002), these rates are characterized by colonization and equilibrium phases. Similarly, diversity-dependent speciation is realized by defining the underlying speciation rate as a function of the total number of species a(N). Other speciation modes such as the point mutation (Hubbell, 2001) are possible to take into account by setting g(0) > 0 as in (Volkov et al., 2003). Peripatry may induce larger skew in the range-size distribution (Pigot et al., 2010), and it can be incorporated by adding a term to Equation (2) that governs the rate of peripatry of a species with a range-size x. Finally, although diversification slowdown is often discussed using reconstructed phylogeny with information extracted from extant species, our discussion is based on the true species number. Developing a method that links to phylogenetic reconstructions would be beneficial for analyzing diversification slowdowns, to directly link to the existing insights.
The shape of species-range-size distributions helps us understand the structure of biodiversity across scales. Species-range-size distributions are intertwined with other macroecological and community patterns such as species-area and endemic-area relationships and species abundance distributions (Takashina et al., 2019).
Thus, further developing our understanding of the dynamics of the range-size distribution will also promote integrated understanding of many basic and applied macroecological questions.

ACK N OWLED G M ENTS
We would like to thank A. Chao, A. Carnaval, and S. J. Franks for their thoughtful comments. NT was supported by Grant-in-Aid for the Japan Society for the Promotion of Science (JSPS) Fellows. Financial support was provided by JSPS (No. 17K15180 to EPE;No. 21K17913 to NT). EPE and NT were additionally supported by subsidy funding to OIST.

CO N FLI C T O F I NTE R E S T
The authors declare that they have no conflict of interest.  Noonan (1990).

A PPE N D I X 1 G A I N -LOSS -A LLO PATR I C S PECI ATI O N PRO CE SS
Let us assume that each species has a geographic range size with state r at time t that experiences gain in its range size by one unit (r → r + 1) or loss (r → r − 1) at rates g and l, respectively. In addition, we assume that allopatric speciation can occur at a size-dependent rate a r , and it produces two smaller range sizes r 1 and r 2 with one new species so as to satisfy r 1 + r 2 = r. All combinations of split size are equally likely. These events and probabilities are summarized as follows: We define an extinction as the event that a species achieves range size 0, which can only occur from r = 1. These events are translated into the following stochastic process that describes the dynamics of the expected number of species with range size r at time t: where g and l are the gain and loss rate of geographic range, and a r are the size-dependent allopatric speciation rate. Here, we perform a continuum limit of Equation (1).

S O M E DY N A M I C A L PRO PE RTI E S O F TH E G L A S PRO CE SS
Since Equation (A2) is linear, we can rewrite the equation with a matrix form where P = (P(1, t), P(2, t), ⋯ , P(r max , t)) T and M is the coefficient matrix of Equation (A2). Therefore, the real part of the dominant has a stable number of species (Re(λ max ) = 0). Figure A1 shows such values. The curve in each panel corresponds to the Re( max ) = 0, and the upper and lower regions are Re(λ max ) > 0 and Re( max ) < 0, respectively.

A PPE N D I X 2 CO NTI N U U M LI M IT O F TH E G L A S PRO CE SS
It is reasonable to consider that the geographic range size changes continuously rather than a discrete manner. The continuum limit of the GLAS process is expected to inherit dynamical property of the GLAS process mentioned above. In the following, we will use a standard procedure to obtain such a model (see, e.g., Gardiner, 2009).
First, let Δx be a small step size, and let x = rΔx. Then, we define the new variable u(x, t) where the number of species between range size x and x + Δx at time t is u(x, t)Δx. By expanding u(x, t) to the second order of Δx, we have the following integro-differential equation: where V = (g − l)Δx and D = g + l 2 Δx 2 , and we replace Σ r � >r with ∫ x max x and a x is the size-dependent split rate.
As in the main text, let us normalize the range size x ∈ [0, 1], and if we assume that a loss of species only occurs when the range size of a species achieves x = 0. To guarantee this condition, we require a Robin boundary condition at the upper boundary and we have the following mixed boundary conditions: where u ′ indicates the derivative with respect to x.

DY N A M I C S O F TOTA L S PECI E S
The dynamics of the total number of species N(t) = ∫ 1 0 u(x, t)dx is derived by integrating both sides of Equation (A4) with respect to x to obtain (A4) Heat map of skewness of the normalized range-size distribution under four speciation scenarios. The left-top region in each panel tends to show the species accumulation at the largest possible range size. The loss rate is l = 1. Parameters for the speciation scenarios (α, β) are (2, 1) (linear increase); (1, 2) (linear increase); (2, 2) (parabola); and (1, 1) (constant) where we have used the condition − Vu(1, t) + Du � (1, t) = 0.
Interchanging the order of the double integral, this becomes The first term of the right-hand side of Equation (A6) corresponds to species loss due to extinction, and the second term represents species increment due to speciation events. Per-species rates of extinction and speciation events are obtained by dividing by the total number of species at given time N(t).
Using Equation (A6), the diversification rate at time t can be computed. Diversification rate at given time (t) is the rate of increase in the log-number of species dlogN(t)∕dt and has the following relationship: (t) = speciation rate per species at time t − extinction rate per species at time t .
It is expected that the diversification rate asymptotically converges to the dominant eigenvalue max = lim t→∞ logN(t)∕t of the coefficient matrix in Equation (A3) as t becomes sufficiently large.

S O LUTI O N O F EQ UATI O N ( A 4) I N A S PECI A L C A S E
It is possible to derive an analytical solution when the speciation scenario is linear increase in a special situation. Although the assumptions allow us to analyze restricted situations, this derivation could help further analytical development of the GLAS process with more general assumptions.
Overall, the following assumptions are made for make the analysis of the Laplace-transformed integration term possible. First, we assume the boundary condition at the upper limit is either Dirichlet boundary condition u(x max , t) = 0 or bounded at infinity lim x→∞ u(x, t) = 0. As an example, we will use the Diriclet boundary condition in the following. Next, we transform the dynamics of species number into fraction of species. To do this, let ũ(x, t)Δx be the fraction of species that range size is between x and x + Δx,

F I G U R E A 3
Heat map of kurtosis of the normalized range-size distribution under four speciation scenarios. The left-top region in each panel tends to show the species accumulation at the largest possible range size. The loss rate is l = 1. Parameters for the speciation scenarios (α, β) are (2, 1) (linear increase); (1, 2) (linear increase); (2, 2) (parabola); and (1, 1) (constant). where the Laplace transformation of integration term is obtained by the integration by parts and used the above-mentioned assumptions ∫ x max 0ũ (x, t)dx = 1 and ũ(x max , t) = 0. Equation (A9) is actually a firstorder linear ordinary differential equation, and hence, it is easy to solve.
With an initial condition that range size of all species is concentrated at a single value ũ(x, 0) = 0 (x 0 ) and Ũ (s, 0) = 1, we obtain the following solution To perform inverse Laplace transform, we further need to arrange Equation (A10): where we used multi-index notations |k| = k 1 + k 2 + k 3 + k 4 , To derive the second line, we applied the multinomial theorem to A. Equation (A10) is now the form to be able to inverse Laplace transform. By applying the inverse Laplace operator to both side, and calculate it separately for each term contains s, we get the final result where Γ(n) is the gamma function.

LO G A R ITH M I C TR A N S FO R M ATI O N O F EQ UATI O N ( A 4)
Since datasets of species-range-size distribution often show a diverse variations of range sizes on the logarithmic axis, it is convenient to employ a log-transformed model in model fitting where range sizes are equally spaced on a logarithmic axis. The logarithmic transformation of the GLAS process is achieved by introducing a new variable y = lnx∕x 0 where x 0 is the reference point. Under this transformation, we describe the number of species between y and y + Δy at time t by v(y, t)Δy, and replacing u (x, t) in Equation (A4) with v(y, t) provides the log-transformed GLAS process: where the boundary condition as Equation (A5) is

F I G U R E A 4
Normalized range-size distributions with (a) no speciation and (b) a = 0.0001. All scenarios cause a negative growth rate of the number of species. We sampled the range-size distribution when either species number becomes 0.01 or rangesize distribution converged to equilibrium regardless of the negative species growth. Other parameter values used are l = 1 and (α, β) are (2, 1) (linear increase); (1, 2) (linear increase); (2, 2) (parabola); and (1, 1) (constant).

(a) (b)
A PPE N D I X 3

TECH N I C A L D E TA I L S O F TH E O P TI M IZ ATI O N A LG O R ITH M
We discuss here the settings of simulated annealing (Kirkpatrick et al., 1983) used in data fitting. This is a heuristic searching process of local optima of an objective function f(S) until the temperature, T, becomes lower than predetermined minimum temperature T min or an objective function (a log-likelihood) meets the predetermined target.
We provide pseudocode as follows: In our case, there are four parameters to optimize (g, s, , ) and these are chosen randomly at each step for updating the value. We heuristically define parameters in the algorithm that induce convergence trends in the optimization (T, T min , Target, c) = (10 2 , 10 −5 , − 1, 0.99) .

F I G U R E A 5
Equilibrium normalized range-size distributions with nonlinear speciation scenarios: Likelihood of speciation monotonically increases with range size (left), monotonically decreases with range size (center), or peaks at intermediate range size (right). The form of each speciation scenarios is shown in the top panels. Skewness and kurtosis of each equilibrium distribution is provided in

F I G U R E A 6
Normalized range-size distributions in case g = l = 1 (left) and g < l = 1 (right). All scenarios cause negative growth rate of the number of species. We sampled the range-size distribution when either the species number becomes 0.01 or range-size distribution converged to equilibrium regardless of the negative species growth rate. Parameters for the speciation scenarios (α, β) are (2, 1) (linear increase); (1, 2) (linear increase); (2, 2) (parabola); and (1, 1) (constant) TA B L E A 1 Summary of (skewness, kurtosis) from the each curve in Figure A5 (α, β) a = 10 −1 a = 10 −2 a = 10 −3