Download Report

Why abundant tropical tree species are
phylogenetically old
Shaopeng Wanga,1, Anping Chenb,1, Jingyun Fanga, and Stephen W. Pacalab,2
a
Department of Ecology, College of Urban and Environmental Sciences, and Key Laboratory for Earth Surface Processes of the Ministry of Education, Peking
University, Beijing 100871, China; and bDepartment of Ecology and Evolutionary Biology, Princeton University, Princeton, NJ 08544
Neutral models of species diversity predict patterns of abundance
for communities in which all individuals are ecologically equivalent. These models were originally developed for Panamanian
trees and successfully reproduce observed distributions of abundance. Neutral models also make macroevolutionary predictions
that have rarely been evaluated or tested. Here we show that
neutral models predict a humped or ﬂat relationship between
species age and population size. In contrast, ages and abundances
of tree species in the Panamanian Canal watershed are found to be
positively correlated, which falsiﬁes the models. Speciation rates
vary among phylogenetic lineages and are partially heritable from
mother to daughter species. Variable speciation rates in an otherwise neutral model lead to a demographic advantage for species
with low speciation rate. This demographic advantage results in
a positive correlation between species age and abundance, as
found in the Panamanian tropical forest community.
|
phylogeny phylogenetic age
niche hypothesis
| Barro Colorado Island (BCI) |
T
he neutral theory of biodiversity (NTB) introduced the idea
that geologic time scales may be directly relevant to ecological population dynamics (1). In NTB models, individuals produce offspring, disperse, and die at random, and species are
ecologically equivalent because they share the same per capita
birth and death probabilities. In most NTB models, new species
arise by a process analogous to mutation; every new offspring has
a (low) probability of mutating into the ﬁrst member of a new
species (1–4). In others, species randomly split into two new
species with a probability proportional to population size (1, 5,
6). NTB models are remarkably successful in predicting distributions of species abundances, particularly with the assumption of mutation speciation (4, 7, 8). Although the assumption of
ecological neutrality has been continuously challenged (9–13),
the predictions of NTB have been shown to be robust to the
existence of niches if species diversity is sufﬁciently high (2, 14,
15), because random drift can still occur between the relative
abundances of species within the same niche or between species
that share very similar niches (16, 17).
For realistically large number of individuals in a region, NTB
models also predict that the abundances of species change slowly
over geologic time before eventually drifting (randomly walking)
to extinction (1, 18, 19). The taxon cycle that has occurred over
intervals on the order of millions of years in several independent
lineages of Lesser Antillean birds (20) provides some indirect
empirical evidence for slow population dynamics over geologic
time scales. However, recent studies also suggested that taxonomic turnover in very abundant clades, like birds (21, 22) and
planktonic foraminifera (23), is sometimes much faster than that
predicted by purely ecological drift.
Here, we focus on geologic time scales and derive the predictions of NTB models for the relationship between a species’
current abundance and its phylogenetic age: the time since the
divergence of a species and its closest extant relative. When a
lineage splits into two species, the phylogenetic ages of both are
set to zero at the time of speciation, because practical phylogenetic
www.pnas.org/cgi/doi/10.1073/pnas.1314992110
trees are constructed based on trait or molecular distance among
lineages, which tells about the divergence time between relatives
but not who is the mother or the daughter. Despite their potential to bridge between ecological and evolutionary theories,
age–abundance relationships have rarely been investigated. NTB
models quantitatively couple population dynamics and speciation and thus allow age–abundance relationships to be quantitatively derived or simulated.
We also investigate the empirical age–abundance relationship
for tree species in the Panama Canal watershed. The abundance
data come from a network of 48 census plots (Fig. S1) containing
593 angiosperm species, among which 530 species are identiﬁed
to species level and thus used in the analysis (Materials and
Methods). Age estimates are derived from a phylogeny for the
1,177 angiosperm species found in this region (Materials
and Methods).
Results and Discussion
NTB models with mutation speciation predict a humped relationship between a species’ abundance and its average age (Fig.
1A). In almost all NTB models, the average rate at which
a species produces a new daughter species is proportional to its
population size (exceptions in refs. 5 and 7). Therefore, abundant species are more likely than rare species to have recently
produced an extant daughter and thus to be phylogenetically
young (Fig. 1A) (6). In contrast, rare species may range from very
young (just been born) to old (unlikely to have a recent daughter
because they are rare) (Fig. S2). The hump in average age at
intermediate abundances occurs because recently produced
daughter species are invariably rare. Also, species with intermediate abundance are less likely to have produced a recent
Signiﬁcance
Neutral population models assume that competing species are
ecologically equivalent and predict several properties of extant
ecosystems. Ecological aspects of neutral models have been
tested repeatedly, with a special focus on tropical trees, but
their macroevolutionary predictions have received little attention. We show that neutral models predict ﬂat or humped
relationships between a species’ abundance and phylogenetic
age. In contrast, abundances and ages of Panamanian tree
species are positively correlated. Similar correlations have been
reported for other groups and explained by deterministic niche
theories. However, we show that neutral models predict positive age–abundance relationships with the addition of speciation rates that vary among lineages.
Author contributions: S.W.P. designed research; S.W., A.C., and J.F. performed research;
S.W. and A.C. analyzed data; and S.W., A.C., and S.W.P. wrote the paper.
The authors declare no conﬂict of interest.
1
S.W. and A.C. contributed equally to the work.
2
To whom correspondence should be addressed. E-mail: [email protected].
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.
1073/pnas.1314992110/-/DCSupplemental.
PNAS Early Edition | 1 of 5
ECOLOGY
Contributed by Stephen W. Pacala, August 12, 2013 (sent for review March 25, 2013)
5
6
7
1
2
3
4
5
6
7
Abundance bin
Fig. 1. Simulated phylogenetic age (generations) and species abundance
under Hubbell’s NTB model with (A) point mutation speciation (per capita
speciation rate: v = 0.002) and (B) random ﬁssion speciation (v = 0.0001).
Solid and dashed lines showed the average age ± SD. The abundance bins in
A are 1, 2–100, 101–300, 301–500, 501–700, 701–1,500, and >1,500 and the
bins in B are 1–50, 51–100, 101–150, 151–200, 201–300, 301–500, and >500.
Other bin widths were also examined and produced the same qualitative
patterns. See SI Text for details of calculations. Scatter plots showing all
species can be found in Figs. S2–S5.
daughter than species with high abundance and thus tend to be
phylogenetically old.
NTB models with the ﬁssion mode of speciation produce ﬂat
age–abundance relationships (Fig. 1B; Fig. S3), because they
create new species with uniform distribution of abundance. Fission speciation divides a population of size N in two, with sizes
N − k and k, where k is uniformly distributed from 1 to N − 1,
and the smaller species is labeled as the daughter species (1, 19).
In addition, humped or ﬂat age–abundance relationships are also
obtained from NTB models in which only a subset of species
have known ages or in which a new daughter species and its
mother would be distinct enough to be classiﬁed as separate
species only after a waiting period (protracted speciation; Figs.
S4 and S5) (3). These results are especially critical for model
testing, because real phylogenies are always incomplete and may
be systematically biased to ignore cryptic young species.
In contrast to the humped or ﬂat age–abundance relationship
predicted by NTB models, the estimated age–abundance relationship for Panamanian trees shows a positive correlation
(Fig. 2). This relationship is statistically signiﬁcant (P < 0.05) by
a variety of tests, using a collection of 1,200 data sets on regional
abundance produced from different permutations of plot data
(Materials and Methods). For example, an ordinary least squares
(OLS) regression of log10(age) vs. abundance has a signiﬁcantly
positive slope for all 1,200 estimates of regional abundance
Mean age
Fraction
6
8
10
100
1.10
20
2
5
5
GLM
OLS
0
200
400
600
Species abundance
800
0.25
4
0.20
2
20
40
60
80 100
Abundance (11 ~ 100)
0.15
10 20
2
5
50
20
100
B
Abundance (1 ~ 10)
2
Phylogenetic age
A
0.30
4
1.30
3
Abundance bin
1.20
2
log10(age)
1
1-10
11-100
Fraction: age < 10My
2.0
1.5
log10(age)
1.0
1.0
0.0
log10(age)
(average P = 0.022). Because the data are highly skewed, we also
performed a Gamma family generalized linear model (GLM)
analysis of the relationship between age and abundance, which
again returned a signiﬁcantly positive slope for all possible datasets
(average P = 0.029).
One complication with the census data is that it covers a relatively wide range of annual precipitation and extends from
forests dominated by drought deciduous species to forests
dominated by tropical evergreens (24). The fact that tropical tree
diversity increases with rainfall has been reported previously
(25–27) and is conﬁrmed here. Within-plot diversity signiﬁcantly
increases with annual rainfall (P < 0.0001), which causes the
average within-plot abundance to decrease with rainfall (P <
0.01). To ensure that effects of rainfall are not solely responsible
for the age–abundance relationship explored in this paper, we
performed a partial correlation analysis of log10(age) and abundance that controlled for both species-speciﬁc mean rainfall (the
cross-plot mean of precipitation times abundance for each species) and species-speciﬁc rainfall range (maximum minus minimum rainfall for plots in which the species is present). The
partial correlation of log10(age) and abundance is signiﬁcant (for
all 1,200 datasets, average P = 0.022) after controlling for both
the mean and range of precipitation. The positive relationship
between age and abundance in the data emerges even though
estimates of both age and abundance are likely to contain large
errors, whose distributions are not possible to be estimated with
conﬁdence and even though the Panama Canal Watershed is not
a random subset of the Neotropics, which could bias the age
estimation of some clades.
The difference between the age–abundance relationships observed in the Panama Canal tree community and those predicted
by the NTB models suggests that it is necessary to modify the
NTB models by introducing some new mechanisms. Here, we
show that a simple extension of the NTB models predicts and
explains the observed positive age–abundance relationship and,
like previous NTB models, is also able to predict observed species–
abundance relationships for Panamanian tree species.
Many previous studies have shown that speciation rates vary
among evolutionary lineages and are at least partly heritable
from mother to daughter species (28–33). We therefore modiﬁed
Hubbell’s NTB model (1) by allowing the per capita speciation
probability ν to vary among species. We also developed this
model to investigate adaptive radiations (34). In each of the
three cases considered, speciation rates are conﬁned to a range
[νmin, νmax]. In the ﬁrst case (Mr), speciation rates vary among
species but are not heritable. When a new daughter species is
2.5
B
2.0
A
>100
Abundance bin
Fig. 2. Relationship between phylogenetic age (My) and relative abundance for Panamanian tree species. (A) Scatter plots of age vs. species relative regional
abundance. Red and blue lines showed the regression lines from ordinary least-square regression for log10(age) vs. abundance (OLS) and generalized linear
model (GLM) for age vs. abundance, respectively. Note the most abundant species Gustaviasuperba (average species abundance = 1,183 and phylogenetic
age = 32.5 My) is not shown in the ﬁgure. (B) The average age and the fraction of young species (phylogenetic age ≤ 10 My) in relative abundance bins: 1–10
(n = 254), 11–100 (n = 193), >100 (n = 48). Other set of bins were also examined and produced the same qualitative patterns.
2 of 5 | www.pnas.org/cgi/doi/10.1073/pnas.1314992110
Wang et al.
Speciation rate
1e-0
1e-05
5000
10000
0
15000 20000
200
1000
2.5
Mr
Ma
Mt
0.5
0.0
1.0
0.5
1.0
log10(age)
2.0
Mr
Ma
Mt
1.5
2.0
800
Abundance
D
1.5
600
400
Abundance
log10(age)
· ðn − 1Þ;
Mr
Ma
Mt
1e-03
1e-0
1e-03
1e-04
0
C
1 − vf
and thus the average speciation rate declines with abundance
and with a slope that is proportional to the variance of speciation
rates among newly born species.
If sufﬁciently strong, the negative slope of vðnÞ should cause
mean phylogenetic age to increase with n as it is observed to do
(Fig. 2A), in contrast to the humped or ﬂat age–abundance
relationships predicted by previously published NTB models
(Fig. 1 A and B). The reason is that the relatively low speciation
rates of abundant species result in relatively long waiting times
between daughters. We do not have analytical formulae for the
age–abundance relationships, but simulations conﬁrm that mean
age is indeed predicted to increase with abundance under ﬁssion
speciation (Fig. 3). Simulations of models Mr, Mt, and Ma with
ﬁssion speciation also produce the kind of upper triangular
scatter exhibited by the age–abundance data on Panamanian
trees (cf. Fig. 2 and Fig. S5). Mutation speciation reduces the
B
Mr
Ma
Mt
σ 2f
1
2
3
4
5
Abundance bin
6
7
1
2
3
4
5
6
7
Abundance bin
Fig. 3. Relationships between average speciation rate and abundance (A and B) and between phylogenetic age (generations) and abundance (C and D),
predicted by the variable speciation-rate models with mutation speciation (A and C) and random ﬁssion speciation (B and D). The trend lines in A were derived
analytically for models Mr and Mt and numerically for model Ma; all trends in B–D are moving averages of the mean values. The bins in C are set to 1, 2–100,
101–300, 301–500, 501–700, 701–1,500, and >1,500 and the bins in D are 1–50, 51–100, 101–150, 151–200, 201–300, 301–500, and >500. Other set of bins were
also examined and produced the same qualitative patterns. Scatter plots showing all species can be found in Figs. S2–S5. See SI Text for details of calculations.
Wang et al.
PNAS Early Edition | 3 of 5
ECOLOGY
vðnÞ ≈ vf − 1e-05
Speciation rate
A
abundance n, and show that it declines as n increases in all our
models (Fig. 3 A and B; Figs. S8 and S9). Abundant species
are predicted to have low speciation rates, simply because
speciation reduces a parent species’ population size by removing the individuals that are transferred to the daughter
species. This mechanism is easy to understand for ﬁssion
speciation but is also remarkably strong for mutation speciation. For example, let f(ν) be the probability density of v for
newly produced daughter species with mean vf and variance
σ 2f . If σ 2f is small, then vðnÞ is approximately (SI Text)
1e-04
produced, its speciation rate is randomly drawn from a constant
uniform probability density over the interval [vmin, vmax]. In the
remaining two cases, the speciation rate v of a new species (vnew)
is uniformly distributed around that of its parental species
(vparent), i.e., vnew ∼ U [vparent − L/2, vparent + L/2], where 1/L is
a measure of heritability. Two alternative boundary conditions
ensure that vnew stays within the interval of [vmin, vmax]. In model
Mt, the uniform distribution is truncated once reaching the
outside of [vmin, vmax], such that vnew ∼ U [vmin, vparent + L/2] if
vparent − L/2 <vmin, and vnew ∼ U [vparent − L/2, vmax] if vparent +
L/2 > vmax. In model Ma, the boundaries are absorbing: any vnew
less than vmin is set to vmin and any larger than vmax is set to
vmax. In addition, Hubbell’s original NTB model is a special form
of our model, in which νmin= νmax (hereafter referred as Mh).
Despite the modiﬁcations of the speciation rate assumptions,
our model is consistent with the assumptions of NTB in that
individuals of all species share the same per capita birth and
death probabilities.
We solved the extended NTB models for S(n, ν), the stationary
ensemble distribution of species with abundance n and speciation rate ν (closed form solutions for models Mr and Mt with
mutation speciation, numerical solution for model Ma with mutation speciation, simulations for ﬁssion speciation; SI Text).
The equilibrium distribution of species abundance (SAD) is
the marginal of S(n, ν) with respect to ν: SðnÞ =
Rsimply
vmax
vmin Sðn; vÞdv. Like Hubbell’s original NTB model, when diversity is large, our models predict log-series SADs under
point mutation and multinomial distributions under random
ﬁssion speciation (Figs. S6 and S7). Using S(n, ν), we can also
calculate vðnÞ, the average speciation rate for species with
amount under mutation speciation. This demographic cost of speciation in our model also makes the community average speciation
rate decline through geologic time (34).
Within a lineage, high speciation rates increase species diversity
and decrease average phylogenetic age. By splitting the abundance of an ancestral species when speciation occurs, high speciation rates also tend to decrease the average population size
within a lineage. Differences in the speciation rate among species thus have population dynamic implications, which are not
ecological niches in the traditional sense but which are strong
enough to alter the relationship between phylogenetic age and
population size. These ﬁndings highlight the importance of incorporating population dynamics in macroevolutionary models
of speciation and extinction.
strongly negative slope of the age–abundance relationship predicted by Mh and increases the chances that a stochastic simulation
will produce a positive slope, although it is not strong enough to
make the expected slope consistently positive (Fig. 3; Fig. S10).
Using an approximation, we also examined the age–abundance
patterns in very large communities over geologic time, which
conﬁrmed our ﬁndings of a ﬂat relationship in model Mh and
positive relationships in models Mr, Ma, and Mt under ﬁssion
speciation (SI Text; Fig. S11).
Like the original NTB models, the modiﬁed NTB models predict
log-series distributions of SAD under point mutation and multinomial distributions under random ﬁssion speciation (Figs. S6 and
S7) when diversity is large. However, compared with the original
NTB models, the variable speciation-rate models predict more
abundant species and fewer species with intermediate abundance
(Figs. S6 and S7). Thus, the new extensions of NTB models explain
how the observed age–abundance relationship could be caused by
an interaction of macroevolutionary processes and random population drift, without compromising the success of the original NTB
models in SAD prediction of species-rich communities.
Previous studies of other groups of species have shown that
broad geographic range is sometimes positively correlated with
phylogenetic age (35–38). The pattern has been explained outside of the context of NTB, with hypotheses about niches. One
hypothesis is that species with unusually broad geographical
niches (i.e., broad climate tolerance) will tend to have broad
ranges and avoid extinction during climatic ﬂuctuations (39, 40).
Another is that species able to attain large ranges might tend to
have both broad niches and the ability to disperse over obstacles,
which might reduce the rate of allopatric speciation and increase
phylogenetic age (41, 42). These hypotheses are also applicable
to the age–abundance relationship because total abundance and
range are frequently correlated (43).
Collectively, our results provide an alternative hypothesis for
the positive correlation between age and abundance, which is
more parsimonious than the niche hypotheses. Old phylogenetic
ages themselves imply rare speciation, and this is enough to
cause the pattern without invoking niche differences. One way to
falsify the speciation rate hypothesis would be to produce positive evidence for the mechanisms of the niche hypotheses. For
example, one might show that abundant species are indeed ecological generalists, perhaps in common garden experiments, and
that the ecological differences among rare and common species
are large enough to overwhelm drift.
Our results show how variation among lineages in speciation
rates can produce relatively large effects on population sizes, if
population dynamics are otherwise neutral. Like many other
traits, the rate of speciation is partially conserved within a lineage, but may change over time (29, 30, 32). These two contrary
tendencies create variation in speciation rate among different
lineages and produces phylogenies in which rapidly speciating
lineages have higher species diversity than slowly speciating ones
(30). Population dynamics in our model are consistent with the
assumptions of NTB in that individuals of all species share the
same per capita birth and death probabilities. However, our
model differs from the original NTB in the sense that high
speciation rates effectively increase realized death rates within
a species under ﬁssion speciation by a tiny amount, by removing
individuals from a species when speciation occurs. Similarly, a high
speciation rate reduces a species’ realized birth rate by a tiny
Estimating Regional Abundance. Regional abundance of each species was
estimated as a sum of its local abundance from the 45 1-ha plots and 1
random 1-ha subplot from each of the 3 larger plots. For completeness, we
constructed 1,200 datasets, using all possible combinations of single hectares
from the three plots >1 ha (1,200 = 50 × 4 × 6). The graph in Fig. 2 used the
subplots located at the southwest corner of each of the three larger plots,
but the results are robust to the choice of the subplots.
1. Hubbell SP (2001) The Uniﬁed Neutral Theory of Biodiversity and Biogeography
(Princeton Univ Press, Princeton, NJ).
2. Chisholm RA, Pacala SW (2010) Niche and neutral models predict asymptotically
4. Volkov I, Banavar JR, Hubbell SP, Maritan A (2003) Neutral theory and relative species
abundance in ecology. Nature 424(6952):1035–1037.
5. Etienne R, Haegeman B (2011) The neutral theory of biodiversity with random ﬁssion
equivalent species abundance distributions in high-diversity ecological communities.
Proc Natl Acad Sci USA 107(36):15821–15825.
3. Rosindell J, Cornell SJ, Hubbell SP, Etienne RS (2010) Protracted speciation revitalizes
speciation. Theoret Ecol 4(1):87–109.
6. Hubbell SP (2003) Modes of speciation and the lifespans of species under
neutrality: A response to the comment of Robert E. Ricklefs. Oikos 100(1):
the neutral theory of biodiversity. Ecol Lett 13(6):716–727.
4 of 5 | www.pnas.org/cgi/doi/10.1073/pnas.1314992110
Materials and Methods
Plot Data. A dataset containing 48 plots of tropical tree communities was
used to test the relationship between phylogenetic age and species abundance (24, 44, 45). The plots are located in the Panama Canal watershed,
which covered a wide range of annual precipitation from 1,760 to 3,750 mm
(Fig. S1). Most (45/48) of these plots are 1 ha in area, and three plots are
larger than 1 ha, including Barro Colorado Island (BCI; 50 ha), Cocoli (4 ha),
and Sherman (5.96 ha). In total, 610 tree size [≥10 cm in breast height diameter (DBH)] species were found in the dataset, including 593 angiosperm
species (530 of which were identiﬁed to species level and therefore used
in the analyses), 1 gymnosperm species (Podocarpus guatemalensis), and 16
unidentiﬁed species (24, 44, 45).
Phylogeny. A tree species list covering the Panama Canal watershed was
obtained from the BCI forest dynamics research project (24, 44, 45), and 1,177
angiosperm species from this list had been phylogenetically classiﬁed as part
of the Angiosperm Phylogeny Group III (APG III) (46) consensus tree. The
topology of the phylogeny for these species was obtained with the online
program Phylomatic (47), which used a based megatree derived from APG III
(46). To estimate the phylogenetic age for each species, we used the module
bladj in the software Phylocom version 4.1 to scale branch length using
known node ages (48). Here we used the age information from Wikström
et al. (49), which estimated divergence times for most angiosperm classes
and orders, as well as some families. This two-step approach for constructing
a dated phylogeny had been frequently used in recent papers (50–52). Finally, the phylogenetic age is deﬁned by the branch length of the terminal
nodes (extant species), which represent the time since the divergence of
a species from its closest extant relative.
ACKNOWLEDGMENTS. We thank R. Chisholm, D. Tilman, and R. Dybzinski
for helpful discussion and W. L. Deng and L. Lin for technical assistance.
Funding was provided by the Carbon Mitigation Initiative (CMI) of the
Princeton Environmental Institute and the National Natural Science Foundation of China (Grant 31021001). The visit of S.W. to Princeton University
was supported by the China Scholarship Council and CMI. The forest plot
data were obtained from the BCI forest dynamics research project, which
was made possible by National Science Foundation grants to Stephen P.
Hubbell, support from the Center for Tropical Forest Science (CTFS), the
Smithsonian Tropical Research Institute, the John D. and Catherine T.
MacArthur Foundation, the Mellon Foundation, the Small World Institute
Fund, and numerous private individuals, and through the hard work of more
than 100 people from 10 countries over the last two decades. The plot project
is part of CTFS, a global network of large-scale demographic tree plots.
193–199.
Wang et al.
Wang et al.
31. Magallón S, Sanderson MJ (2001) Absolute diversiﬁcation rates in angiosperm clades.
Evolution 55(9):1762–1780.
32. Savolainen V, Heard SB, Powell MP, Davies TJ, Mooers AØ (2002) Is cladogenesis
heritable? Syst Biol 51(6):835–843.
33. Sepkoski JJ, Jr. (1998) Rates of speciation in the fossil record. Philos Trans R Soc Lond B
Biol Sci 353(1366):315–326.
34. Wang S, Chen A, Fang J, Pacala SW (2013) Speciation rates decline through time in
individual-based models of speciation and extinction. Am Nat 182(3):E83–E93.
35. Böhning-Gaese K, Caprano T, van Ewijk K, Veith M (2006) Range size: Disentangling
current traits and phylogenetic and biogeographic factors. Am Nat 167(4):555–567.
36. Jablonski D (1986) Larval ecology and macroevolution of marine invertebrates. Bull
Mar Sci 39(2):565–587.
37. Paul JR, Tonsor SJ (2008) Explaining geographic range size by species age: A test
using Neotropical Piper species. Tropical Forest Community Ecology, eds Carson WP,
Schnitzer S (Wiley-Blackwell, Oxford), pp 47–62.
38. Taylor CM, Gotelli NJ (1994) The macroecology of Cyprinella: Correlates of phylogeny,
body size, and geographical range. Am Nat 144(4):549–569.
39. McKinney ML (1997) Extinction vulnerability and selectivity: Combining ecological
and paleontological views. Annu Rev Ecol Syst 28(1):495–516.
40. Purvis A, Jones KE, Mace GM (2000) Extinction. Bioessays 22(12):1123–1133.
41. Gaston KJ, Chown SL (1999) Geographic range size and speciation. Evolution of Biological Diversity, eds Magurran AE, May RM (Oxford Univ Press, Oxford, UK), pp
237–259.
42. Jablonski D, Roy K (2003) Geographical range and speciation in fossil and living
molluscs. Proc Biol Sci 270(1513):401–406.
43. Brown JH (1984) On the relationship between abundance and distribution of species.
Am Nat 124(2):255–279.
44. Hubbell SP, Condit R, Foster R (2005) Barro Colorado Forest Census Plot Data.
Available at http://ctfs.arnarb.harvard.edu/webatlas/datasets/bci. Accessed February
2, 2012.
45. Condit R (1998) Tropical Forest Census Plots (Springer-Verlag and R. G. Landes Company, Berlin).
46. Webb CO, Donoghue MJ (2005) Phylomatic: Tree assembly for applied phylogenetics.
Mol Ecol Notes 5(1):181–183.
47. The Angiosperm Phylogeny Group (2009) An update of the Angiosperm Phylogeny
Group classiﬁcation for the orders and families of ﬂowering plants: APG III. Bot J Linn
Soc 161(2):105–121.
48. Webb CO, Ackerly DD, Kembel SW (2008) Phylocom: Software for the analysis
of phylogenetic community structure and trait evolution. Bioinformatics 24(18):
2098–2100.
49. Wikström N, Savolainen V, Chase MW (2001) Evolution of the angiosperms: Calibrating the family tree. Proc Biol Sci 268(1482):2211–2220.
50. Cadotte MW, Cardinale BJ, Oakley TH (2008) Evolutionary history and the effect of
biodiversity on plant productivity. Proc Natl Acad Sci USA 105(44):17012–17017.
51. Kembel SW, Hubbell SP (2006) The phylogenetic structure of a neotropical forest tree
community. Ecology 87(7, Suppl):S86–S99.
52. Liu X, et al. (2012) Experimental evidence for a phylogenetic Janzen-Connell effect in
a subtropical forest. Ecol Lett 15(2):111–118.
PNAS Early Edition | 5 of 5
ECOLOGY
7. Etienne R, Apol M, Olff H, Weissing FJ (2007) Modes of speciation and the neutral
theory of biodiversity. Oikos 116(2):241–258.
8. Volkov I, Banavar JR, Hubbell SP, Maritan A (2007) Patterns of relative species
abundance in rainforests and coral reefs. Nature 450(7166):45–49.
9. Fargione J, Brown CS, Tilman D (2003) Community assembly and invasion: An experimental test of neutral versus niche processes. Proc Natl Acad Sci USA 100(15):
8916–8920.
10. Levine JM, HilleRisLambers J (2009) The importance of niches for the maintenance of
species diversity. Nature 461(7261):254–257.
11. Ricklefs RE, Renner SS (2012) Global correlations in tropical tree species richness and
abundance reject neutrality. Science 335(6067):464–467.
12. Turnbull LA, Manley L, Rees M (2005) Niches, rather than neutrality, structure
a grassland pioneer guild. Proc Biol Sci 272(1570):1357–1364.
13. Wills C, et al. (2006) Nonrandom processes maintain diversity in tropical forests. Science 311(5760):527–531.
14. Chave J, Muller-Landau HC, Levin SA (2002) Comparing classical community models:
theoretical consequences for patterns of diversity. Am Nat 159(1):1–23.
15. Purves DW, Pacala SW (2005) Ecological drift in niche-structured communities: Neutral pattern does not imply neutral process. Biotic Interactions in the Tropics, eds
Burslem D, Pinard M, Hartley S (Cambridge Univ Press, Cambridge, UK), pp 107–138.
16. Tilman D (2004) Niche tradeoffs, neutrality, and community structure: A stochastic
theory of resource competition, invasion, and community assembly. Proc Natl Acad Sci
USA 101(30):10854–10861.
17. Vellend M (2010) Conceptual synthesis in community ecology. Q Rev Biol 85(2):
183–206.
18. Leigh EG, Jr. (1981) The average lifetime of a population in a varying environment.
J Theor Biol 90(2):213–239.
19. Ricklefs RE (2003) A comment on Hubbell’s zero-sum ecological drift model. Oikos
100(1):185–192.
20. Ricklefs RE, Bermingham E (1999) Taxon cycles in the Lesser Antillean avifauna. Ostrich 70(1):49–59.
21. Halley JM, Iwasa Y (2011) Neutral theory as a predictor of avifaunal extinctions after
habitat loss. Proc Natl Acad Sci USA 108(6):2316–2321.
22. Ricklefs RE (2006) The uniﬁed neutral theory of biodiversity: Do the numbers add up?
Ecology 87(6):1424–1431.
23. Allen AP, Savage VM (2007) Setting the absolute tempo of biodiversity dynamics. Ecol
Lett 10(7):637–646.
24. Hubbell SP, et al. (1999) Light-Gap disturbances, recruitment limitation, and tree diversity in a neotropical forest. Science 283(5401):554–557.
25. Cayuela L, Benayas JMR, Justel A, Salas-Rey J (2006) Modelling tree diversity in
a highly fragmented tropical montane landscape. Glob Ecol Biogeogr 15(6):602–613.
26. Gentry AH (1982) Patterns of neotropical plant species diversity. Evol Biol 15(1):1–84.
27. Givnish TJ (1999) On the causes of gradients in tropical tree diversity. J Ecol 87(2):
193–210.
28. Coyne JA, Orr HA (2004) Speciation (Sinauer Associates, Sunderland, MA).
29. Davies TJ, et al. (2004) Darwin’s abominable mystery: Insights from a supertree of the
angiosperms. Proc Natl Acad Sci USA 101(7):1904–1909.
30. Heard SB (1996) Patterns in tree balance with variable and evolving speciation rates.
Evolution 50(6):2141–2148.
Supporting Information
Wang et al. 10.1073/pnas.1314992110
SI Text
Many of the elements in this SI Text were also included in
the supplemental material for Wang et al. (1), which examined
the same model to study adaptive radiations. These elements
are included here so that all information necessary to understand the age-abundance results is accessible to readers in
one place.
Model and Analytic Solutions. We modiﬁed Hubbell’s neutral
theory of biodiversity (NTB) model (2–4) by allowing per capita
speciation rate (v) to vary among species within the restricted
range [vmin, vmax], where 0 < vmin < vmax << 1. Speciation rate
hereafter always refers to the per capita probability of speciation unless otherwise speciﬁed. Under the mode of point
mutation speciation, we are able to derive the equilibrium
distribution of species abundance. We ﬁrst derive the frequency distribution of the per capita speciation rate and then
use it to obtain the species abundance distribution (see below). For random ﬁssion speciation, however, we cannot solve
the equilibrium distribution of species abundance analytically,
and thus the results for random ﬁssion speciation are based on
simulations.
Frequency distribution of per capita speciation rate. First, let quv be the
transition probability that a daughter species has a per capita
speciation rate v given that its parent had
a per capita speciaR vmax
quv dv = 1. In pretion rate u. Therefore, for any ﬁxed u, vmin
vious NTB models, speciation rates are not variable, so these
transition probabilities are irrelevant. If speciation rate is variable but not heritable (model Mr in the text), then quv is independent of u.
We now consider the probability that the population with
per capita speciation rate v in a closed and fully occupied
community of J individuals increases or decreases by one
individual given a single death-birth event in the community.
We denote gv as the probability density of individuals with
speciation rate v, and N v Δv = g v ΔvJ as the discrete population size of the individuals whose speciation rate is within v ±
Δv/2. Following a death, the probability of increasing population size by one individual with speciation rate within v ±
Δv/2 is
Probability ðincreasing population size by one individual with
speciation rate within v ± Δv=2Þ = ð1 − Nv Δv=JÞ
Zvmax Nu
p ðNv Δv=JÞ p ð1 − vÞ + ð1 − Nv Δv=JÞ p
· u · quv Δv du:
J
vmin
On the right hand, the ﬁrst term is the joint probability that an
individual with speciation rate not equal to v dies and is replaced with a new individual whose mother had speciation
rate v and who did not mutate into a new species at birth.
The second terms describe the probability that an individual with speciation rate not equal to v dies and is replaced
with a new individual whose speciation rate is v through
speciation.
Similarly, the probability that the community decreases by
one individual with speciation rate within v ± Δv/2 following
a death is
Wang et al. www.pnas.org/cgi/content/short/1314992110
Probability ðdecreasing population size by one individual
Nv Δv
Nv Δv
with speciation rate within v ± Δv=2Þ =
· 1−
J
J
"
Z
· 1−
+
vmax
vmin
Nu
Nv Δv
· v · qvv Δv
u · ðquv ΔvÞdu −
J
J
Nv Δv
1−
J
#
Nv Δv 2
· v ·ð1 − qvv ΔvÞ:
J
At equilibrium, probability (increasing population size by one individual with speciation rate within v ± Δv/2) = probability (decreasing population size by one individual with speciation rate
within v ± Δv/2), and replacing Nv/J with gv, we have
Zvmax
ð1 − gv ΔvÞ · gv Δv · ð1 − vÞ + ð1 − gv ΔvÞ ·
Z
= gv Δv · ð1 − gv ΔvÞ ·
1−
½gu · u · quv Δvdu
!
vmin
vmax
vmin
gu u ·ðquv ΔvÞdu − gv Δv · v · qvv Δv
1 − gv Δv
+ ðgv ΔvÞ2 · v · ð1 − qvv ΔvÞ:
This equation can be simpliﬁed as
Zvmax
gv v =
gu u · quv du:
[S1]
vmin
In addition, Eq. S1 could also be more easily derived if we
consider the problem from the perspective of speciation
equilibrium. When a community reaches equilibrium, speciation events that increase in the fraction of individuals with
a given speciation rate will be balanced by events that decrease it:
Zvmax
ðgu u · quv ΔvÞdu − gv Δv · v · qvv Δv = gv Δv · v · ð1 − qvv ΔvÞ:
vmin
Using Eq. S1, we can solve the equilibrium of gv after substituting in the functional form of quv (given below). For the
nonheritable case (Mr) and heritable case with truncated boundary (Mt), analytic solutions were obtained; for the heritable case
with absorbing boundary (Ma), we get numerical results.
Mr: Nonheritable random speciation rate. In the nonheritable case,
the speciation rate of a newly produced species is uniformly
sampled from a probability density h(v) on the interval [vmin,
vmax]. Hence, the transition probability (quv) is simply h(v) and
therefore equation may be written
gv v = v · hðvÞ;
where
1 of 13
difﬁcult to solve for gv analytically because of the δ function.
However, we can obtain a numerical solution from the discrete
form of Eq. S1:
Zvmax
v=
gu udu:
vmin
From
R vmax
vmin
vmin
v
dv
distribution quv = hðvÞ = vmax −1 vmin , we have
1
:
gv =
v · ðln vmax − ln vmin Þ
The average speciation rate of the community is v =
vmax − vmin
ln vmax − ln vmin .
[S2]
R vmax
vmin
gv vdv =
Mt: Heritable speciation rate with truncated boundary. In the heritable case, the speciation rate of new species (vnew) is uniformly
distributed around that of its parental species (vparent), i.e.,
vnew ∼ U[vparent − L/2, vparent + L/2]. To ensure that vnew stays
within the interval of [vmin, vmax], in model Mt, the uniform
distribution is truncated once reaching the outside of [vmin, vmax],
such that vnew ∼ U[vmin, vparent + L/2] if vparent − L/2<vmin, and
vnew ∼ U[vparent − L/2, vmax] if vparent + L/2>vmax. Therefore, for
any v satisfying ju − vj ≤ L=2, the transition probability is
8
>
>
>
<
quv =
1
:::::::::::::::::u < vmin + L=2
L=2 + u − vmin
1
L
>
>
>
:
:::::::::::::::::u ∈ ½vmin + L=2; vmax − L=2
1
::::::::::::::::::u > vmin − L=2:
L=2 + vmax − u
Substituting it into Eq. S1, we obtain
8 C 1 v−v
>
>
< v · 2+ L
min
gv =
gvi vi =
gv dv = 1, we have v = R vmax1hðvÞ . In the case of a uniform
>
>
:
::::::::::::::::v < vmin + L=2
C
:::::::::::::::::::::v ∈ ½vmin + L=2; vmax − L=2
v
C 1 vmax − v
· +
::::::::::::::::v > vmax − L=2:
v 2
L
[S3]
R vmax
gv dv =
Here, C is the normalization parameter to ensure that vmin1
ðvmax − L=2Þ
vmin
vmax
1
+
ln
+
v
ln
1; hence, C1 = 12 ln vvmax
v
min
max
vmin R+ L=2
vmax − L=2 .
L
min ðvmin + L=2Þ
vmax
gv vdv = C · ðvmax −
The average per capita speciation rate v = vmin
vmin − L=4Þ.
Ma: Heritable speciation rate with absorbing boundary. Model Ma is
also a heritable case, which differs from Mt only in the treatment
of boundary constraints. Here the boundaries are absorbing: any
vnew less than vmin is set to vmin, and any vnew larger than vmax is
set to vmax. Therefore, for any u satisfying ju − vj ≤ L=2, the
transition probability is
1 u − vmin
−
:::::::::::::::::v = vmin
δv−vmin
2
L
8
>
>
<
quv =
>
>
:
1
::::::::::::::::::::::::::::::::::::::::v ≠ vmin ; vmax
L
1 vmax − u
−
:::::::::::::::::v = vmax
δv−vmax
2
L
Here quv follows a mixed distribution. The probabilities at v =
vmin or vmax are denoted by a delta function δxR that returns to
zero if x ≠ 0 and inﬁnity at x = 0, which satisﬁes f0g δx dx = 1. It is
Wang et al. www.pnas.org/cgi/content/short/1314992110
X
gvj vj · qvj vi :
[S1′]
j
Here, we divide [vmin, vmax] into K equal segments, and gvi and
qvj vi are the mean of gv and quv for each segment. We replace the
left side of Eq. S1′ with a vector w = ðw1 ; w2 ; :::; wK Þ, where
wi = gvi vi . Hence, Eq. S1′ can be written in matrix form as
w = wQ, where Q is a K*K transition matrix. Because the transition probability is absorbed at the boundary, but equally distributed within [vparent − L/2, vparent + L/2] when (vparent − L/2) and
(vparent + L/2) do not reach beyond the boundaries, the transition
matrix Q is deﬁned as
0
W
BW −1
B
B ⋮
B
B 1
B
B
B
B
B
B
1 B
B
Q¼
2W − 1 B
B
B
B
B
B
B
B
B
B
B
@
1
1
⋮
1
1
⋯
⋯
⋯
⋯
1
1
⋮
1
1
1
1
⋮
1
1
⋱
⋱
⋯
⋯
1
1
⋱
1
⋱
⋯
1 ⋯
⋱
1
1
⋮
1
1
1
⋮
1
1
⋯
⋯
⋯
⋯
1
1
⋮
1
1
C
C
C
C
C
C
C
C
C
C
C
C
C;
C
C
C
C
C
C
C
1 C
C
⋮ C
C
W −1A
W
where 2W − 1 = K· vmax L− vmin is the length (number) of possible
discrete segments that daughter species may fall in. Hence, Q
has an eigenvalue of 1, and w is the relevant left eigenvector. Because of the scarcity of Q (W << K), w can be solved very efﬁciently in MatLab, even when K is as large as 105. After solving
XK
g =
for w, we then obtain g = ðgv1 ; gv2 ; :::; gvK Þ from (1 =
i = 1 vi
XK wi
XK
) and the average speciation rate v =
g v.
i=1 v
i = 1 vi i
i
Species abundance distribution. To derive the equilibrium distribution of species abundance in a neutral community, a commonly
used formulation is to consider the one-step stochastic process
governed by a master equation (4–6). Under the mode of point
mutation speciation, the master equation governing the probability that a run of the model contains a species with population size n and speciation rate v in the community at time t,
S(n, v, t)Δv, is
d½Sðn; v; tÞΔv
= In−1;v Sðn − 1; v; tÞΔv + Dn+1;v Sðn + 1; v; tÞΔv
dt
− In;v + Dn;v Sðn; v; tÞΔv:
[S4]
In,v and Dn,v are one-step transition probabilities for a species
with abundance n and speciation rate v to increase or decrease
population size by one individual given a death because of ecological drift. When n > 0, we have In;v = J −J n · J −n 1 · ð1 − vÞ and
Dn;v = nJ JJ −− n1 + nJ −− 11 p v (6). When n = 1, in the ﬁrst item on the
right side of Eq. S4, I0,v represents the probability of speciation
that produces
R vmax a new species with speciation rate v. Therefore,
I0;v Δv = vmin
gu u · quv Δvdu = gv vΔv according to Eq. S1. Finally,
with the boundary condition I-1,v = 0, D0,v = 0, IJ,v = 0, and
= 0 for all n > 0, we can
DJ+1,v = 0, at equilibrium when d½Sðn;v;tÞΔv
dt
easily obtain Dn;v Sðn; vÞΔv − In−1;v Sðn − 1; vÞΔv = 0. Consequently,
the equilibrium of S(n, v,t) is
2 of 13
R vmax
n−1
n−1 I
Ii;v
1
i;v
Sðn; vÞΔv = ∏
Δv = I0;v · ∏
·
Δv
i = 0 Di+1;v
i = 1 Di;v Dn;v
n−1
= gv v · ∏
i=1
vmin
and σ 2f , respectively. From Eq. S1, when J is large, S(n, v) can be
approximated by (7)
1−v
1
Δv
· 1 + iJ−−1i v nJ· JJ −− n1 + nJ −− 11 v
gv v · ð1 − vÞn−1 n−1
= nJ · J − n n − 1 · ∏
i=1
J −1 + J −1 v
1
1 + iJ−−1i v
f ðvÞdv = 1. The mean and variance of f(v) are denoted as vf
[S5]
J
Sðn; vÞΔv = · gv v · ð1 − vÞn−1 Δv:
n
Δv:
Thus, we have
Z vmax
It is also easily shown that the equilibrium S(n, v)Δv is also the
probability that a randomly chosen species in a randomly chosen
run of the model has abundance n. Therefore, we obtained the
equilibrium distribution of abundance that is simply the marginal
of S(n, v) with respect to v:
vðnÞ =
Z
Sðn; vÞ · vdv
vmin
Z
Zvmax
SðnÞ =
vmin
vmin
Substituting the analytic or numerical solution of gv into
Eqs. S5 and S6, we obtained stationary distribution of S(n, v)
and S(n).
The species abundance distributions derived above for mutation speciation and simulated for ﬁssion species are similar
to those from the original NTB models, which predict a log
series distribution for mutation speciation and a multinomial
distribution for ﬁssion speciation (Figs. S6 and S7) when diversity is large. However, compared with the original NTB
models, the variable speciation rate models predict more abundant species and fewer species with intermediate abundance
(Figs. S6 and S7).
Speciation rate–abundance distribution. Based on Eq. S6, we can also
calculate the average per capita speciation rate for the species
with abundance n:
Z vmax
Z vmax
Sðn; vÞ · n · vdv
Sðn; vÞ · vdv
= vmin
:
[S7]
vðnÞ = vmin
SðnÞ · n
SðnÞ
2
2
2
The variance of vðnÞ
R can also be derived: σ v ðnÞ = E½v ðnÞ − vðnÞ ,
vmax
vmin
vmax
v
= Z minvmax
[S6]
Sðn; vÞdv:
vðnÞ ≈ vf −
σ 2f
2
σ2
= vf − 2f
·
1 − vf
n−1
−2
· 2n
=v
ð1 − vf Þ f
Z
vðnÞ = 1 − Z
Wang et al. www.pnas.org/cgi/content/short/1314992110
vmax
vmin
vmax
−
ð1 − vf Þ
J
·gv v · ð1 − vÞn−1 dv
n
gv v2 · ð1 − vÞn−1 dv
gv v · ð1 − vÞn−1 dv
:
[S9]
f ðvÞ · ð1 − vÞn dv
f ðvÞ · ð1 − vÞn−1 dv
:
[S10]
vmin
Under the assumption that σ 2f is small, we expanded (1 − v)n to
the second order around v = vf
Zvmax
Zvmax
n
f ðvÞ · ð1 − vÞ dv ≈
vmin
h
n n−1
f ðvÞ · 1 − vf − v − vf · n · 1 − vf
vmin
+
ðv − vf Þ
2
2
n−2 i
dv
nðn − 1Þ 1 − vf
n σ 2
n−2
= 1 − vf + 2f nðn − 1Þ 1 − vf
:
Therefore,
Z
vðnÞ = 1 − Z
vmax
vmin
vmax
f ðvÞ · ð1 − vÞn dv
f ðvÞ · ð1 − vÞn−1 dv
vmin
n σ 2f
n−2
1 − vf + nðn − 1Þ 1 − vf
2
=1−
:
n−1 σ 2f
n−3
+ ðn − 1Þðn − 2Þ 1 − vf
1 − vf
2
We further expanded this expression to the ﬁrst order around
σ 2f
2
=0
n−2 n
n−3
nðn − 1Þ 1 − vf
− 1 − vf ðn − 1Þðn − 2Þ 1 − vf
2n−2
1 − vf
σ 2f
· gv v · ð1 − vÞn−1 · vdv
Substituting gv*v = f(v)/C0 into Eq. S9, we have
Sðn;vÞ · v2 dv
where E½v2 ðnÞ = vmin SðnÞ
:
Numerical and simulation results showed that the average
speciation rate generally declined with abundance under all cases
of the variable speciation rate model (Fig. 3; Figs. S6 and S7).
Here, to investigate how average speciation rate changes with
abundance analytically, we study Eq. S7 under some restrictions.
We ﬁrst derive a general analytical approximation for vðnÞ assuming small variance of the speciation rate. We then derive
a general expression for vðnÞ under model Mr in which Vnew is
uniformly distributed.
If f(v) is the probability density of v for newly produce
daughterR species, then, using the derivation of gv above,
vmax
gu uquv du = C0 · gu u, where C0 is determined by
f ðvÞ = C0 vmin
J
vmin n
Z vmax
=
SðnÞ
vmax
[S8]
[S11]
· ðn − 1Þ:
3 of 13
Eq. S11 clearly demonstrates that vðnÞ will decrease with n
when speciation rate v is variable.
Next, for model Mr with a uniform distribution for f(v), Eq. S9
becomes
Z vmax
v · ð1 − vÞn−1 dv
vmin
vðnÞ = Z vmax
ð1 − vÞn−1 dv
vmin
n ð1 − vmin Þn vmin − ð1 − vmax Þn vmax
1
=
+
·
:
ð1 − vmin Þn − ð1 − vmax Þn
n+1
n+1
[S12]
By studying the ﬁrst derivative, we can prove that Eq. S12 is
a monotonically decreasing function of n.
Model Simulations. To investigate the relationship between speciation rate, phylogenetic age, and species abundance, we performed simulations for Mh, Mr, Ma, and Mt. In all simulations,
we set the community size at J = 20,000. We also set the speciation rate heritability parameter L at L = 0.001 for Ma and Mt.
With point mutation speciation, we simulated three levels of speciation rates for Mh, v = 0.0005, 0.002, and 0.01, and three speciation rate ranges [vmin, vmax] for Mr, Ma, and Mt: small (0.001,
0.0035), intermediate [0.0001, 0.01], and large [0.00001, 0.015].
With random ﬁssion speciation, we did parallel simulations with
parameters v = 0.00001, 0.0001, and 0.001 for Mh and [vmin, vmax] =
[0.001, 0.005], [0.0005, 0.005], and [0.00001, 0.01] for Mr, Ma,
and Mt. For Mh, simulations were run for 50,000 generations (1
generation = J random death and replacement steps), which is
sufﬁcient for convergence. For other models, simulations were
run until the average speciation rate stopped decreasing, which
required up to 1,000,000 generations. Each of the 24 cases was
repeated 200 times.
To examine the relationships between phylogenetic age and
abundance, we plotted the average age and its SD in different
abundance bins. For simulations with mutation speciation (Figs. 1
and 3; Figs. S2 and S4), the bins are 1, 2–100, 101–300, 301–500,
501–700, 701–1,500, and >1,500; for ﬁssion speciation (Figs. 1
and 3; Figs. S3 and S5), bins are 1–50, 51–100, 101–150, 151–200,
201–300, 301–500, and >500. Other sets of bins were examined
and produced the same qualitative patterns.
In addition, we simulated age–abundance relationships in randomly chosen subsets of the species present (20%) to mimic incomplete phylogenies and reconstructed the phylogeny in the
model based only on the subsample (Figs. S4 and S5). We also
1. Wang S, Chen A, Fang J, Pacala SW (2013) Speciation rates decline through
time in individual-based models of speciation and extinction. Am Nat 182(3):
E83–E93.
2. Hubbell SP (2001) The Uniﬁed Neutral Theory of Biodiversity and Biogeography
(Princeton Univ Press, Princeton, NJ).
3. Hubbell SP (2003) Modes of speciation and the lifespans of species under neutrality: A
response to the comment of Robert E. Ricklefs. Oikos 100(1):193–199.
4. Volkov I, Banavar JR, Hubbell SP, Maritan A (2003) Neutral theory and relative species
abundance in ecology. Nature 424(6952):1035–1037.
Wang et al. www.pnas.org/cgi/content/short/1314992110
simulated protracted speciation by assuming that new daughter
species are recognizable as separate species only after they have
diverged from their parents for at least 10 generations (200,000
random death and replacement steps). New species younger than
10 generations were merged with their parental species when
calculating the phylogeny and abundances.
We also used an approximate algorithm (1, 8) to study the age–
abundance relationship at realistically large temporal and spatial
scales under the random ﬁssion speciation. Speciﬁcally, the
change of population size from time T to time T + 1 for each
species can be approximated as a Brownian motion process, with
the mean change equal to zero and variance proportional to the
duration (1,000 generations in our simulations) and initial population size (NT). Meanwhile, there is a small probability (i.e.,
the product of total births 1,000NT and per capita speciation rate
v) that the species undergoes ﬁssion and produces a new species
(population size: n0). The approximate algorithm is
NT+1 =
Rnormð NT ; 1;000NT · s Þ; with probability : 1 − 1;000NT · v
:
Rnormð NT − n0 ; 1;000NT · s Þ; with probability : 1;000NT · v
where Rnorm(u, w) represents the normal distribution with
mean u and variance w, and s denotes the per capita (demographic and environmental) stochasticity. Without considering
environmental stochasticity, s is equal to 2 because each individual undergoes on average two events of birth and/or death in one
generation.
With this approximate algorithm, we simulated a neutral community with size J = 1010, a range of per capita speciation rates
[vmin, vmax] = [10−18, 10−11], and heritability parameter L = vmax/5.
For each simulation, the community starts from one species with
per capita speciation rate v0 = 5 × 10−12 and is run for 5,000
time steps (i.e., ﬁve million generations). Without environmental
stochasticity, some species abundances in simulations sometimes
reached ten millions (107), and some phylogenetic ages reached
millions of generations. Like the exact results from simulations
of small neutral communities, the age–abundance relationships
from the approximations of large communities were ﬂat in model
Mh and increasing in the extended models (Mr, Ma, and Mt; Fig.
S11). When environmental stochasticity was included (8), average species abundance increased and average phylogenetic age
decreased. However, the predictions for the shape of the age–
abundance relationship (ﬂat or increasing) were qualitatively
unchanged (Fig. S11).
5. McKane AJ, Alonso D, Solé RV (2000) Mean-ﬁeld stochastic theory for species-rich
assembled communities. Phys Rev E Stat Phys Plasmas Fluids Relat Interdiscip Topics
62(6 Pt B):8466–8484.
6. Vallade M, Houchmandzadeh B (2003) Analytical solution of a neutral model of
biodiversity. Phys Rev E Stat Nonlin Soft Matter Phys 68(6 Pt 1):061902.
7. Alonso D, McKane AJ (2004) Sampling Hubbell’s neutral theory of biodiversity. Ecol
Lett 7(10):901–910.
8. Allen AP, Savage VM (2007) Setting the absolute tempo of biodiversity dynamics. Ecol
Lett 10(7):637–646.
4 of 13
Fig. S1. Spatial distribution of the 48 plots in the Panama Canal watershed. The map is adapted from the Center for Tropical Forest Study (CTFS) website
(https://ctfs.arnarb.harvard.edu/Public/Datasets/PanamaPlotSummary/PanamaPlots.html).
Wang et al. www.pnas.org/cgi/content/short/1314992110
5 of 13
5
1
3
5
1
3
5
Abundance bin
log10(age)
1.0
6
2.0
3
5
7
0
3
5
7
1
3
5
7
log10(age)
2.0
2.0
Abundance bin
1.0
log10(age)
7
5
1.0
log10(age)
log10(age)
2.0
1
0.0
2.0
Mt
1.0
log10(age)
Large
1
Abundance bin
0.0
Intermediate
4
Abundance bin
1.0
log10(age)
7
Abundance bin
Small
7
0.0
2.0
Ma
1.0
Large
log10(age)
Intermediate
5
Abundance bin
0.0
Small
3
3
0.0
2.0
1.0
log10(age)
1
7
Abundance bin
2
Abundance bin
3
3
1
2
1
7
0.0
2.0
Mr
1.0
Large
log10(age)
Intermediate
5
Abundance bin
0.0
Small
3
0.0
2.0
1.0
1
7
1
5
1
3
5
Abundance bin
7
1.0
3
Abundance bin
0.0
1
log10(age)
Mh
0.0
3.0
v=0.01
1.5
v=0.002
0.0
v=0.0005
log10(age)
B
A
1
3
5
7
Abundance bin
Fig. S2. (A) Scatter plots of simulated phylogenetic age (generations) and abundance with mutation speciation for 200 runs of each model. For model Mh,
three levels of speciation rate were simulated: v = 0.0005, 0.002, and 0.01; for models with variable speciation rates (Mr, Ma, and Mt), three speciation rate
ranges (vmin, vmax), small (0.001, 0.0035), intermediate (0.0001, 0.01), and large (0.00001, 0.015), were simulated. Other parameters are community size J =
20,000 and heritability parameter L = 0.001. All simulations were repeated 200 times. (B) Mean age (solid lines) along with abundance (x axis) binned from
rarest to the most abundant by 1, 2–100, 101–300, 301–500, 501–700, 701–1,500, and >1,500. The dashed lines mark ±SD.
Wang et al. www.pnas.org/cgi/content/short/1314992110
6 of 13
2
3
4
1
2
3
4
5
1
2
3
4
Abundance bin
5
6
2.0
1
3
5
7
log10(age)
2.5
1
2
3
4
5
6
1
log10(age)
1
2
3
4
5
Abundance bin
3
5
7
Abundance bin
1.5
log10(age)
5
4
Abundance bin
0.5
Mt
1.5
log10(age)
Large
6
Abundance bin
0.5
Intermediate
5
1.5
log10(age)
6
Abundance bin
Small
4
0.5
Ma
1.5
log10(age)
Large
3
Abundance bin
0.5
Intermediate
2
3
0.5
log10(age)
1
2
Abundance bin
1.5
log10(age)
5
Abundance bin
Small
1
2.0
1
7
0.5
Mr
1.5
Large
log10(age)
Intermediate
5
Abundance bin
0.5
Small
3
1.5
log10(age)
1
0.5
2.0
7
0.5
5
6
2.0
3
Abundance bin
0.5
1
log10(age)
Mh
1.0
v=0.001
2.5
v=0.0001
log10(age)
B
v=0.00001
1.5
A
1
3
5
7
Abundance bin
Fig. S3. (A) Scatter plots of simulated phylogenetic age (generations) and abundance with ﬁssion speciation for 200 runs of each model. For model Mh, three
levels of speciation rate were simulated: v = 0.00001, 0.0001, and 0.001; for models with variable speciation rates (Mr, Ma, and Mt), three speciation rate ranges
(vmin, vmax), small (0.001, 0.005), intermediate (0.0005, 0.005), and large (0.00001, 0.01), were simulated. Other parameters are community size J = 20,000 and
heritability parameter L = 0.001. All simulations were repeated 200 times. (B) Mean age (solid lines) along with abundance (x axis) binned from rarest to the
most abundant by 1–50, 51–100, 101–150, 151–200, 201–300, 301–500, and >500. The dashed lines mark ±SD.
Wang et al. www.pnas.org/cgi/content/short/1314992110
7 of 13
7
1
7
1
3
5
1
1
3
5
Abundance bin
7
2
3
log10(age)
Sampling: 20%
1
3
5
Abundance bin
3
5
7
Abundance bin
1
1.0
log10(age)
Mt
1
2.0
7
4
2.0
5
Abundance bin
0.0
log10(age)
3
7
2.5
log10(age)
3
7
5
Age > 10
3.0
5
3.0
4.0
5
Sampling: 20%
Abundance bin
3
Abundance bin
1
log10(age)
3
2
log10(age)
1
3
2.0
log10(age)
1
7
0
1
7
Age > 10
Abundance bin
Ma
5
1.0
3.5
Sampling: 20%
Abundance bin
3
Abundance bin
1.0
5
1
7
2.0
log10(age)
1.0
0.0
3
5
0.5
2.0
Mr
1
3
Abundance bin
7
Age > 10
2.0
5
1.0
3
Abundance bin
log10(age)
log10(age)
0.0
1
Age > 10
1.0
3.0
log10(age)
Sampling: 20%
1.5
2.0
Mh
1.0
log10(age)
B
0.0
A
1
3
5
7
Abundance bin
Fig. S4. (A) Scatter plots of simulated age (generations) and abundance under mutation speciation (ﬁrst column), mutation speciation with only a subset of
20% of species sampled (second column), and with protracted mutation speciation, where daughter species were considered to be indistinguishable form their
parents for 10 generations (third column). For Mh, we used the speciation rate v = 0.002; for Mr, Ma, and Mt, we used the scenario with largest variation in
speciation rate in Fig. S2. (B) Mean age (solid lines) along with abundance (x axis) binned from rarest to the most abundant by 1, 2–100, 101–300, 301–500, 501–
700, 701–1,500, and >1,500. The dashed lines mark ±SD.
Wang et al. www.pnas.org/cgi/content/short/1314992110
8 of 13
7
1
7
1
3
5
1
7
3.0
2.0
log10(age)
1
7
3.0
log10(age)
Sampling: 20%
1
3
5
Abundance bin
3
5
7
Abundance bin
2.5
log10(age)
5
7
1.0
2.0
log10(age)
0.5
3
Abundance bin
5
7
Age > 10
Abundance bin
Mt
1
3
5
1.0
Sampling: 20%
Abundance bin
3
Abundance bin
2.5
log10(age)
5
3.0
1
4.0
3
2.0
log10(age)
7
1.0
2.0
log10(age)
0.5
1
7
Age > 10
Abundance bin
Ma
5
1.0
3.0
Sampling: 20%
Abundance bin
3
Abundance bin
2.0
log10(age)
5
5
1.0
2.0
log10(age)
0.5
3
3
Abundance bin
Mr
1
log10(age)
1.0
1
7
2.0
3.5
2.5
log10(age)
5
Abundance bin
7
Age > 10
2.0
3
Age > 10
1.0
1
Sampling: 20%
1.5
Mh
2.0
log10(age)
B
1.0
A
1
3
5
7
Abundance bin
Fig. S5. (A) Scatter plots of simulated age (generations) and abundance under ﬁssion speciation (ﬁrst column), ﬁssions speciation with only a subset of 20% of
species sampled (second column), and with protracted ﬁssion speciation, where daughter species were considered to be indistinguishable form their parents
for 10 generations (third column). For Mh, we used the speciation rate v = 0.0001; for Mr, Ma, and Mt, we used the scenario with largest variation in speciation
rate in Fig. S3. (B) Mean age (solid lines) along with abundance (x axis) binned from rarest to the most abundant by 1–50, 51–100, 101–150, 151–200, 201–300,
301–500, and >500. The dashed lines mark ±SD.
Wang et al. www.pnas.org/cgi/content/short/1314992110
9 of 13
4096
1
4
16 64
512
4 16
128
40
1
2048
4 16
128
2048
abundance
I
Species richness
Mt
Mh
0
0
64 256
30
2048
10 20 30 40 50
Species richness
40
30
20
10
16
abundance
8 10
0
1
Species richness
4
Ma
Mh
abundance
0
1
2048
6
Species richness
30
4096
H
128
F
abundance
G
4 16
abundance
20
Species richness
30
20
10
0
4
1
0 5 10
40
E
1
20
512 4096
abundance
D
Species richness
16 64
4
512
2
16 64
abundance
10 20 30 40 50
4
0
0
1
Mr
Mh
10
20
30
40
Species richness
C
10
10
20
30
40
Species richness
B
0
Species richness
A
1
4
16 64
512 4096
abundance
1
4 16
128
2048
abundance
Fig. S6. Comparing the species abundance distributions for variable speciation rate models and the original neutral theory of biodiversity models (Mh), under
mutation speciation. Thick bars are averages from simulations for each variable speciation rate model capped by thin bars showing 1 SD. Green, orange, and
blue lines shows our analytic solutions, whereas black lines show the log-series curve of model Mh with same average speciation rate (4). Results for models Mr,
Ma, and Mt are shown in A–C, D–F, and G–I, respectively. Results under scenarios A, D, and G correspond to small; B, E, and H correspond to intermediate; and C,
F, and I correspond to large variation in speciation rate, which are deﬁned in the legend for Fig. S2.
Wang et al. www.pnas.org/cgi/content/short/1314992110
10 of 13
1 2 4 8
128
32
250
1
4
64
32
abundance
64 256
2048
128
350
250
Mh
150
150
250
Species richness
I
1 2 4 8
16
abundance
0
8 16
abundance
250
128
50
Species richness
250
150
0 50
4
256
150
Species richness
1 2 4 8
H
2
64
Mh
abundance
G
16
0 50
250
128
abundance
1
4
F
150
Species richness
250
150
0 50
32
1
abundance
E
1 2 4 8
150
128
0 50
D
Species richness
32
abundance
0 50
32
abundance
Species richness
Species richness
150
0
1 2 4 8
Mh
0 50
250
C
50
150
250
Species richness
B
0 50
Species richness
A
1
4
16
64
256
abundance
Fig. S7. Comparing the species abundance distributions for variable speciation rate models and the original neutral theory of biodiversity models (Mh), under random
ﬁssion speciation. Thick bars are averages from simulations for each variable speciation rate model capped by thin bars showing 1 SD. Black lines show multinomial
curve of model Mh with same average speciation rate (5). Results for models Mr, Ma, and Mt are shown in A–C, D–F, and G–I, respectively. Results under scenarios A, D,
and G correspond to small; B, E, and H correspond to intermediate; and C, F, and I correspond to large variation in speciation rate, which are deﬁned as in Fig. S3.
Fig. S8. Relationship between speciation rate and abundance under mutation speciation. Points are simulated data (200 repeats of each model), and red lines showed
analytic solutions of the corresponding model (SI Text). Results for models Mr, Ma, and Mt are shown in A–C, D–F, and G–I, respectively. Results under scenarios A, D, and
G correspond to small; B, E, and H correspond to intermediate; and C, F, and I correspond to large variation in speciation rate, which are deﬁned as in Fig. S2.
Wang et al. www.pnas.org/cgi/content/short/1314992110
11 of 13
Slope
Ma
Mt
Mh
Mr
Ma
Mt
No
Small Intermediate Large
-0.002
0.002
Mh
Mr
0.006
B
-0.0020
Slope
A
-0.0005 0.0005
Fig. S9. Relationship between speciation rate and abundance under random ﬁssion speciation. Points are simulated data (200 repeats of each model), and red
lines show the moving average. Results for models Mr, Ma, and Mt are shown in A–C, D–F, and G–I, respectively. Results under scenarios A, D, and G correspond
to small; B, E, and H correspond to intermediate; and C, F, and I correspond to large variation in speciation rate, which are deﬁned as in Fig. S3.
No
Small Intermediate Large
Variation in speciation rate
Variation in speciation rate
Fig. S10. Distribution of regression slopes from ordinary least squares (OLS) regressions [log10(age) vs. abundance] of the simulated data from each of the 200
runs summarized in Figs. S2 and S3 under mutation (A) and ﬁssion speciation (B). See the main text for model deﬁnition.
Wang et al. www.pnas.org/cgi/content/short/1314992110
12 of 13
Fig. S11. Age–abundance relationship approximately simulated for large communities without (A: s = 2) and with environmental stochasticity (B: s = 100).
Parameters: J = 1010, [vmin, vmax] = [10−18, 10−11], L = vmax/5. Abundance bins are 1–2 × 106, 2 × 106–4 × 106, 4 × 106–6 × 106, 6 × 106–8 × 106, and >8 × 106 for A
and 1–2 × 107, 2 × 107–4 × 107, 4 × 107–6 × 107, 6 × 107–8 × 107, 8 × 107–108, and >108 for B. Solid lines and dashed lines showed the average and SD of (logtransformed) age in each abundance bin.
Wang et al. www.pnas.org/cgi/content/short/1314992110
13 of 13