Introduction

The concept of beta-diversity was introduced by Whittaker (1960, 1972) to define the variation in community composition among sites in a region. Several other beta-diversity indices have subsequently been proposed (reviewed by Koleff et al. 2003; Anderson et al. 2011; Jost et al. 2011; Legendre and De Cáceres 2013). Recent studies have focused on temporal changes in community composition at both single sites and multiple sites surveyed repeatedly over time (Magurran 2011; Legendre and Condit 2019). Temporal changes in community composition are referred to as “temporal beta-diversity” (Hatosy et al. 2013; Legendre and Gauthier 2014; Shimadzu et al. 2015), which is simply an extension of spatial beta-diversity. Only a few studies have focused on methodological developments of temporal beta-diversity (Shimadzu et al. 2015; Legendre 2019; Nakadai 2020). The persistence and turnover of individuals over time is a key feature of temporal beta-diversity (Magurran et al. 2019; Nakadai 2020). The speed and frequency of compositional change over time is associated with the speed of individual turnover, and must be considered in community comparisons because even randomly high individual turnover can result in high temporal beta-diversity (Nakadai 2020).

To address this issue, Nakadai (2020) proposed the “individual-based beta-diversity” concept, and provided novel indices to evaluate individual turnover and compositional shift by comparing individual turnover between two periods at a given site. Specifically, I developed an individual-based temporal beta-diversity approach for the Bray–Curtis dissimilarity index (i.e., an abundance-based temporal beta-diversity index) (Nakadai 2020). Over time, some individuals replace others, and therefore, ecological communities are dynamic and vary to some degree both spatially and temporally (Mori et al. 2018; Tatsumi et al. 2019). In addition, some individuals of long-lived species (e.g., trees and vertebrates) persist over the long term, whereas individuals of short-lived species are replaced with conspecifics or individuals of other species (Nakadai 2020).

The average values of pairwise dissimilarity have been widely used to evaluate compositional heterogeneity among multiple sites but cannot accurately reflect the overall compositional heterogeneity within a pool of more than two sites (Baselga 2013a). Diserud and Ødegaard (2007) first proposed a multiple-site (or more generally, a multiple-unit) similarity measure, which has conceptual and methodological advantages over conventional approaches based on average pairwise similarities (Koleff and Gaston 2002; Gaston et al. 2007; Baselga et al. 2007). Many dissimilarity indices have been extended to multiple-unit indices (Baselga 2013b, 2017). Multiple-unit dissimilarity indices can take into account components shared among more than three units (i.e., sites or temporal units) (Baselga et al. 2007). The pairwise individual-based indices proposed by Nakadai (2020) may be limited by the typical problems associated with the application of multiple-unit datasets. In particular, even if the numbers of mortal and recruited individuals are same, the presence of long-lived individuals could result in underestimates of individual turnover when the average values of the pairwise indices were used. Extending individual-based indices to multiple-temporal datasets may allow accurate assessment of the temporal variability of community composition based on individual turnover.

A null model is a pattern-generating model based on randomization of ecological data and is designed to produce a pattern that would be expected in the absence of a particular ecological mechanism (Gotelli and Graves 1996; Gotelli and McGill 2006). Null models for temporal patterns in ecological communities have been evaluated previously (e.g., cyclic shift permutations; Hallett et al. 2014, 2016). Such models use a matrix of the abundance, presence or absence, or biomass of species over multiple-temporal units (e.g., years) as an input. In each realization of the algorithm, the time series of every species is shifted forward by a random number of units, independently of other species (Kalyuzhny 2020). However, no null model approach using individual randomization has been developed for community dynamics, and no research has been conducted on the patterns of individual persistence using either beta-diversity indices or null models.

Therefore, the purpose of this paper was to develop new multiple-unit indices that can quantify compositional variability across time according to the speed of individual turnover, as well as novel random permutation criteria to evaluate the deviation in individual persistence from the perspective of ecological drift. In the present study, I extended individual-based beta-diversity indices to multiple units after briefly reviewing the history of pairwise and multiple-unit indices and also introduced a newly developed null model. Then, I applied these indices to data from a 50-ha forest dynamics plot on Barro Colorado Island (BCI) in Panama (Condit et al. 2019), as a case study to demonstrate the usage and properties of indices compared to average pairwise measures. Using both novel indices and a null model provides opportunities to answer the following questions pertaining to the BCI plot. Are differences in the degree of compositional variability explained only by differences in the degree of individual turnover? Does the degree of deviation in individual persistence from each expected value of ecological drift lead to differences in compositional variability among forest plots?

Materials and methods

Conventional incidence-based and abundance-based indices and novel individual-based indices

Two types of dissimilarity indices have been proposed in previous studies: incidence-based and abundance-based dissimilarity indices (Baselga 2013b; Nakadai, 2020). Specifically, the widely used Bray–Curtis dissimilarity index is an abundance-based extension of the Sørensen index (Legendre and Legendre 2012; Baselga 2013b). The Bray–Curtis dissimilarity index was originally proposed as a percentage difference index by Odum (1950), yet it is commonly misattributed to Bray and Curtis (1957). The incidence-based and abundance-based methods target the number of species and the number of individuals, respectively. Conventional beta-diversity indices are based on intersection (A, a) and relative complement (B, b and C, c) components. For example, the Sørensen dissimilarity index (dsor; Sørensen 1948) is formulated as follows:

$${d}_{sor}= \frac{b+c}{2a+b+c},$$
(1)

where a is the number of species common to both units (i.e., sites or temporal units), b is the number of species that occur in the first but not the second unit, and c is the number of species that occur in the second but not the first unit. Similarity, the Bray–Curtis dissimilarity index (Odum 1950) is formulated as follows:

$${d}_{BC}= \frac{B+C}{2A+B+C}.$$
(2)

The abundance-based components are as follows:

$$A= \sum_{i}\mathrm{min}({x}_{ij},{x}_{ik}),$$
(3a)
$$B= \sum_{i}\left[{x}_{ij}-\mathrm{min}({x}_{ij},{x}_{ik})\right],$$
(3b)
$$C= \sum_{i}\left[{x}_{ik}-\mathrm{min}({x}_{ij},{x}_{ik})\right],$$
(3c)

where xij is the abundance of species i at unit (i.e., site or time) j, and xik is the abundance of species i at unit (i.e., site or time) k. Therefore, A is the number of individuals of each species common to both units j and k, whereas B and C are the numbers of individuals unique to units j and k, respectively.

Therefore, it should be possible to distinguish the community dynamics component of persistence (P) from those that change over time, i.e., mortality and recruitment (M and R, respectively). M refers to individuals occurring at time 1 (abbreviated T1) that die before time 2 (abbreviated T2), thereby representing the death of individuals during the time between T1 and T2. R identifies individuals that were not alive at T1 (or were not yet in the size class counted, e.g., diameter breast height ≥ 1 cm) but were subsequently alive and counted at T2, thereby representing R between T1 and T2. Here, the term “recruitment (R)” is used because I have tree communities in mind. “Colonization” in the sense of Tatsumi et al. (2021) could be employed when the approach introduced in this study is applied to studies on taxa other than trees (e.g., bird ringing data). Component P refers to the persistence of individuals from T1 to T2. Although P, M, and R relate to components A, B, and C, respectively, here I have used different characters to emphasize the processes represented by each component. The respective abundances are calculated as follows:

$$P= \sum_{o}\mathrm{min}({z}_{o1},{z}_{o2}),$$
(4a)
$$M= \sum_{o}\left[{z}_{o1}-\mathrm{min}({z}_{o1},{z}_{o2})\right],$$
(4b)
$$R= \sum_{o}\left[{z}_{o2}-\mathrm{min}({z}_{o1},{z}_{o2})\right],$$
(4c)

where zo1 is individual o at time T1, and zo2 is individual o at time T2. Both zo1 and zo2 take a value only of 1 or 0 (presence and absence, respectively). Therefore, P is the total number of individuals present at both T1 and T2, whereas M and R are the numbers of individuals unique to T1 and T2, respectively (Nakadai 2020).

T1 and T2 can be generalized as unit j (Tj) and unit k (Tk), and generalized versions of P, M, and R for units j and k are simply formulated as follows:

$$P= \sum_{o}\mathrm{min}({z}_{oj},{z}_{ok}),$$
(5a)
$$M= \sum_{o}\left[{z}_{oj}-\mathrm{min}({z}_{oj},{z}_{ok})\right],$$
(5b)
$$R= \sum_{o}\left[{z}_{ok}-\mathrm{min}({z}_{oj},{z}_{ok})\right],$$
(5c)

which are used for the main analysis of this study.

This formulation of the individual-based temporal beta-diversity index for Bray–Curtis (dMR) dissimilarity can be expressed as follows:

$${d}_{MR}= \frac{M+R}{2P+M+R},$$
(6)

In some cases, there is turnover of different individuals belonging to the same species even when the community composition is stable over time, which contributes to a dynamic compositional equilibrium across time periods. Therefore, the components that change over time (M and R) are re-organized into two further components: compositional equilibrium (E) and shift (B and C), respectively (Nakadai 2020; Fig. 1). Component M can be partitioned into lost individuals that contribute to compositional change in a community (B) and lost individuals replaced by conspecifics, contributing substantially to equilibrium in community composition (Eloss) (Fig. 1). Similarly, component R can be partitioned into gained individuals that contribute to compositional change in a community (C) and gained individuals replacing conspecifics, thereby contributing substantially to equilibrium in community composition (Egain) (Fig. 1). By definition, the numbers of Eloss and Egain individuals are identical; hence, I replaced both Eloss and Egain with the coefficient E, i.e., the component that contributes to dynamic compositional equilibrium (Nakadai 2020). Furthermore, components M and R can also be analyzed in the same manner as community data (Nakadai 2020). Therefore, the Bray–Curtis dissimilarity between communities can be calculated based on M and R as follows:

$${v}_{s}= \frac{B+C}{2E+B+C}=\frac{B+C}{M+R},$$
(7)

where vs indicates the speed of compositional shifts in a community relative to the total speed of individual turnover associated with M and R. If vs is 0, then individual turnover would not contribute to compositional shifts in a community; i.e., M and R are identical, and an equilibrium state exists. If vs is 1, all individual turnovers contribute to compositional shifts, and M and R are completely different. It is also possible to calculate vs as dBC divided by dMR, as follows:

$$d_{BC} \div d_{MR} = \frac{B + C}{{2A + B + C}} \div \frac{M + R}{{2P + M + R}} = \frac{B + C}{{2\left( {P + E} \right) + B + C}} \div \frac{M + R}{{2P + M + R}} = \frac{B + C}{{2P + M + R}} \div \frac{M + R}{{2P + M + R}} = \frac{B + C}{{2E + B + C}} = { }v_{s} .$$
(8)
Fig. 1
figure 1

Comparison of species composition and individuals (20 individuals of four species total) among temporal units, showing the components (A, B, C, M, R, P, Eloss, and Egain) used for assessing compositional variability over time. The squares, circles, pentagons, and hexagons indicate species 1, 2, 3, and 4, respectively. The numbers under each element (e.g., 12) represent the corresponding time step. For components A and P, the number 123 is also shown, indicating the shared species and individuals across the three time steps, respectively. Detailed explanations of each component are provided in Table 1. The dashed lines with numbers (i.e., a denominator and a numerator) indicate that the numerator chosen from among candidates of the same species (i.e., the number in the denominator) will be Eloss or Egain because the components of Eloss and Egain do not require individual identity information

Full size image
Table 1 Explanations of the components used in the conventional and novel indices of beta-diversity
Full size table

Conventional multiple-unit beta-diversity indices

Multiple-site dissimilarity indices were first proposed by Diserud and Ødegaard (2007), which made it possible to consider species information shared among more than three units using the inclusion–exclusion principle (Erickson 1996). Since then, various multiple-unit beta-diversity indices have been proposed (Baselga et al. 2007; Baselga, 2013b), including several new multiple-unit indices based on pairwise dissimilarity indices (Baselga et al. 2007; Baselga 2010, 2012, 2017). Specifically, these multiple-unit dissimilarity indices share the A (a), B (b), and C (c) components of indices based on abundance (incidence). amu, bmu, and cmu correspond to a, b, and c, respectively, of incidence-based pairwise dissimilarity indices.

Although both bmu and cmu are simply sums of all components divided into two parts in the pairwise indices, amu is based on the inclusion–exclusion principle described above and is difficult to understand intuitively. For example, in the case of three units, the component amu can be calculated as follows (Baselga et al. 2007; Baselga, 2010):

$${a}_{12}+{a}_{23}+{a}_{13}-\mathrm{min}\left({a}_{12},{a}_{23},{a}_{13}\right)=\sum_{j<k}{a}_{jk}- {a}_{123},$$
(9)

where ajk is the number of species common to both units. Considering a general case with n units, the component amu can be formulated as follows:

$$\sum_{j<k}{a}_{jk}- \sum_{j<k<l}{a}_{jkl}+ \sum_{j<k<l<m}{a}_{jklm}-\bullet \! \, \bullet \! \, \bullet .$$
(10)

Equation 9 can be converted based on the inclusion–exclusion principle, and the resulting components are as follows:

$${a}_{mu}=\sum_{j}{S}_{j}-{S}_{\mathrm{T}}$$
(11a)
$${b}_{mu}=\sum_{j<k}\mathrm{min}\left({p}_{jk},{p}_{kj}\right)$$
(11b)
$${c}_{mu}=\sum_{j<k}\mathrm{max}\left({p}_{jk},{p}_{kj}\right)$$
(11c)

where Sj is the total number of species at site j, ST is the total number of species at all sites considered together, and pjk and pkj are the numbers of species exclusive to units j and k, respectively, when compared pairwise (Baselga 2010). Specifically, pjk is the number of individuals included in unit j but not in unit k, and vice versa.

Similar to the case of amu described above, the component Amu can be calculated for n units as follows (Baselga 2010):

$${\sum }_{i}^{S}\left[\sum_{j<k}\mathrm{min}({x}_{ij},{x}_{ik})- \sum_{j<k<l}\mathrm{min}\left({x}_{ij},{x}_{ik},{x}_{il}\right)+\sum_{j<k<l<m}\mathrm{min}\left({x}_{ij},{x}_{ik},{x}_{il},{x}_{im}\right)-\sum_{j<k<l<m<n}\mathrm{min}\left({x}_{ij},{x}_{ik},{x}_{il},{x}_{im},{x}_{in}\right)\bullet \! \, \bullet \! \, \bullet ,\right].$$
(12)

Applying the inclusion–exclusion principle allows the Amu component to be simplified. Amu, Bmu, and Cmu correspond to A, B, and C, respectively, from abundance-based pairwise dissimilarity indices:

$${A}_{mu}={T}_{A\mathrm{B}}- {\sum }_{i}^{S}\mathrm{max}({x}_{ij},{x}_{ik},{x}_{il},{x}_{im},{x}_{in}\dots )$$
(13a)
$${B}_{mu}=\sum_{j<k}\mathrm{min}\left({\sum }_{i}^{S}\left[{x}_{ij}-\mathrm{min}({x}_{ij},{x}_{ik})\right], {\sum }_{i}^{S}\left[{x}_{ik}-\mathrm{min}({x}_{ij},{x}_{ik})\right]\right)$$
(13b)
$${C}_{mu}=\sum_{j<k}\mathrm{max}\left({\sum }_{i}^{S}\left[{x}_{ij}-\mathrm{min}({x}_{ij},{x}_{ik})\right], {\sum }_{i}^{S}\left[{x}_{ik}-\mathrm{min}({x}_{ij},{x}_{ik})\right]\right)$$
(13c)

where TAB is the total abundance in the data set (details see Baselga, 2017).

Multiple-unit indices for Sørensen and Bray–Curtis are formulated as follows (Baselga 2010, 2013b):

$${d}_{sor.mu}= \frac{{b}_{mu}+{c}_{mu}}{2{a}_{mu}+{b}_{mu}+{c}_{mu}}$$
(14)
$${d}_{BC.mu}= \frac{{B}_{mu}+{C}_{mu}}{2{A}_{mu}+{B}_{mu}+{C}_{mu}}.$$
(15)

Novel individual-based multiple-unit beta-diversity indices

In this study, I developed individual-based components and novel indices for multiple-unit data sets. In accordance with existing methods for the calculation of multiple-unit dissimilarity, the persistence (Pmu) component can be calculated for three time steps as follows:

$$\mathrm{min}\left({z}_{o1},{z}_{o2}\right)+\mathrm{min}\left({z}_{o2},{z}_{o3}\right)+\mathrm{min}\left({z}_{o1},{z}_{o3}\right)-\mathrm{min}\left({z}_{o1},{z}_{o2},{z}_{o3}\right).$$
(16)

Considering a general case with n units, the component Pmu can be formulated as

$${\sum }_{o}^{I}\left[\sum_{j<k}\mathrm{min}({z}_{oj},{z}_{ok})- \sum_{j<k<l}\mathrm{min}\left({z}_{oj},{z}_{ok},{z}_{ol}\right)+\sum_{j<k<l<m}\mathrm{min}\left({z}_{oj},{z}_{ok},{z}_{ol},{z}_{om}\right)-\sum_{j<k<l<m<n}\mathrm{min}\left({z}_{oj},{z}_{ok},{z}_{ol},{z}_{om},{z}_{on}\right)\bullet \! \bullet \!\, \bullet ,\right].$$
(17)

Using the inclusion–exclusion principle, Pmu can be simplified and the components Pmu, Mmu and Rmu reformulated as follows:

$${P}_{mu}=\sum_{j}{I}_{j}-{I}_{\mathrm{T}}$$
(18a)
$${M}_{mu}=\sum_{\mathrm{j}<k}\mathrm{min}\left({q}_{jk},{q}_{kj}\right)$$
(18b)
$${R}_{mu}=\sum_{\mathrm{j}<k}\mathrm{max}\left({q}_{jk},{q}_{kj}\right).$$
(18c)

where Ij is the total number of individuals in a unit j, IT is the total number of individuals in all units considered together, and qjk and qkj are the numbers of individuals exclusive to units j and k, respectively, when compared pairwise.

The multiple-unit individual-based index is formulated as follows:

$${d}_{MR.mu}= \frac{{M}_{mu}+{R}_{mu}}{2{P}_{mu}+{M}_{mu}+{R}_{mu}}.$$
(19)

Furthermore, as an extension of Eq. 8, I also introduced vs.mu, which indicates the compositional variability in a community relative to the individual turnover, as follows:

$$v_{s.mu} = d_{BC.mu} \div d_{MR.mu} .$$
(20)

Novel randomization criteria to construct a null model of individual-tracked monitoring data

Here, I introduce a novel method for constructing a null model to evaluate the degree of individual persistence over time. Figure 2 is a conceptual diagram outlining the basis of the new null model. Cases (i) and (ii) shown in Fig. 2 illustrate the patterns of individual persistence of two individuals. The degree of individual persistence differs between the two cases, even if the numbers of individuals involved in recruitment and mortality were equal across the time steps. In both cases, one individual was recruited, and one individual died. Specifically, in case (i), the individual z1 persists over time, whereas individual z2 was recruited between T1 and T2 and died between T2 and T3. On the other hand, in case (ii) described below, individual z1 died between T2 and T3, and z2 was recruited between T1 and T2 and persisted, at least until the end of observation at T4. Such differences in the patterns of individual persistence lead to differences in the value of the newly developed individual-based temporal multiple-unit dissimilarity index (dMR.mu), even when the numbers of individuals undergoing mortality and recruitment are equal. The component Pmu can account for differences in the pattern of individual persistence. Specifically, Pmu decreases, and the index dMR.mu increases, when the number of long-lived individuals (i.e., more than three time steps) increases, and vice versa. The null model approach introduced here uses those properties of the component and the index.

Fig. 2
figure 2

Conceptual diagram showing the basis of the new null model. Both (a) case (i) and (b) case (ii) showed patterns of individual persistence for two individuals. The pattern of individual persistence differed between the two cases, even if the numbers of individuals affected by recruitment and mortality were the same at each time step. Here, in both cases, one individual was recruited, and one individual died. Specifically, in case (i), individual z1 persisted over time, but individual z2 was recruited between times 1 and 2 and died between times 2 and 3. On the other hand, in case (ii), individual z1 died between times 2 and 3, and individual z2 was recruited between times 1 and 2 and continued to persist at least until time 4. The average value of the individual-based temporal beta-diversity index (\(\overline{{d_{{{\text{MR}}}} }}\)), which indicates the degree of individual turnover, may be the same in both cases, but the newly developed multiple-unit dissimilarity index will have different values (dMR.mu) because the component Pmu can account for differences in the pattern of individual persistence. Specifically, Pmu decreases, and the index dMR.mu increases, when the number of long-lived individuals (i.e., those persisting for more than three time steps) increases. The sprout and the skull and crossbones symbols indicate recruitment and mortality, respectively

Full size image

In the novel null model, the simulations use the original numbers of recruited and mortal individuals, but mortal individuals are randomly selected at each time step. The simulations were repeated (e.g., 100 times) for each plot. The index dMR.mu can be calculated from each repeated simulation. By comparing the observed and simulated values of the index dMR.mu, the standardized effect size of the degree of deviation in individual persistence from the null model (Pses.all) can be calculated as the observed value minus the mean of the null distribution, divided by the standard deviation in the null distribution. Calculating the deviation from the null model is commonly performed to express biological differences regardless of the units of the indices (e.g., McCabe et al. 2012; Nakadai and Kawakita 2016). The value of Pses.all is close to 0, indicating random deaths of individuals (i.e., ecological drift). In cases in which the value of Pses.all is greater than 0, the proportion of long-lived individuals (more than three time steps) is larger than that in the null model.

Case study: the permanent 50-ha forest dynamics plot on BCI in Panama

To demonstrate the use of the new individual-based multiple-unit indices, I analyzed the data from Condit et al. (2019). The permanent 50-ha forest dynamics project plot on BCI in Panama was established in 1981 by Robin Foster and Stephen Hubbell (Hubbell and Foster 1983). Data from this plot have been used in many scientific papers, and descriptions of the area and survey methods are widely available (e.g., Legendre and Condit 2019). This 50-ha plot has been subjected to eight detailed surveys since its inception: in 1982, 1985, 1990, 1995, 2000, 2005, 2010, and 2015. The BCI data were divided into 1,250 quadrats (20 m × 20 m). Harms et al. (2001) classified quadrats into six habitat zones and a zone of mixed habitats (Fig. 3a, Table S1). The six habitat zones are grouped 1–6 and the mixed habitat zone is group 7 (Harms et al. 2001). Detailed information on each group and a map are shown in Fig. 3a.

Fig. 3
figure 3

Comparison of the conventional multiple-unit Bray–Curtis dissimilarity (dBC.mu; b and f), individual-based temporal beta-diversity dissimilarity (dMR.mu; c and g), composition variability (vs.mu; d and h), and the degree of deviation in individual persistence from the null model (Pses.all; e and i) among habitat zones (a) using the pairwise Wilcoxon rank sum test, which allows for pairwise comparisons between groups with correction for multiple testing. Habitat types labeled with the same letter did not significantly differ in terms of that index based on the Wilcoxon rank sum test with Holm’s method for multiple comparisons. Detailed information about the habitat zones and results of the Wilcoxon rank sum test are shown in Tables S1 and S3, respectively. In (a–e), the horizontal and vertical axes indicate easting and northing, respectively

Full size image

First, I calculated the values of the two individual-based beta-diversity indices for all possible pairs of time points and then determined their average values (\(\overline{{d }_{MR}}\), \(\overline{{v }_{s}}\)) for each quadrat. In addition, I calculated the newly developed multiple-unit indices (dMR.mu, vs.mu) for each quadrat and tested their relationships using simple linear regression to characterize the novel indices. I predicted that the calculated average values would always be lower than the novel multiple-unit indices based on the results of Baselga (2013a), who compared the patterns of the average values of conventional beta-diversity indices and multiple-unit indices. Second, because high individual turnover can result in high compositional variability within a community, to explore how much the values deviate from predictions based on individual turnovers, I tested the relationship between newly developed multiple-unit indices (dMR.mu) and two previous ones (dBC.mu, dsor.mu). Third, I calculated the degree of deviation in individual persistence from the null model (Pses.all) based on 100 randomizations via a novel random permutation method (details described above) for each quadrat. I analyzed differences among habitat types in compositional variability (dBC.mu), the degree of individual turnover (dMR.mu), compositional variability in individuals undergoing turnover (vs.mu), and the degree of deviation in individual persistence from the null model (Pses.all) using the Wilcoxon rank sum test, which allows for pairwise comparisons between groups with correction for multiple testing. The p-values were adjusted by the Holm method. This analysis seeks to afford novel insights on how BCI data can be interpreted using indices based on individual identity information. Finally, because the patterns of individual turnover and persistence may affect community variability, to confirm the effect of individual persistence on community compositional variability, I tested the relationship between the degree of deviation in individual persistence from the null model (Pses.all) and the conventional multiple-unit dissimilarity index (dBC.mu) using a generalized additive model fitted by the restricted maximum likelihood (REML) method. I predicted that the conventional multiple-unit dissimilarity index (dBC.mu) would decrease with increasing deviation of individual persistence from that of the null model (Pses.all).

All analyses were conducted in R software (ver. 4.0.3; R Development Core Team, 2020). The betapart R package (Baselga and Orme 2012) was used to calculate multiple-unit versions of the Bray–Curtis and Sørensen dissimilarity indices. The R package mgcv (Wood 2017) was used to fit the generalized additive model to the relationship between the degree of deviation in individual persistence from the null model (Pses.all) and the conventional multiple-unit dissimilarity index (dBC.mu). The function pairwise.wilcox.test was used to compare the differences among the habitats.

Results

Simple linear regressions showed that dMR.mu increased significantly with increasing \(\overline{{d }_{MR}}\) (adjusted R2 = 0.971, p < 0.0001; Fig. S1a) and vs.mu increased significantly with increasing \(\overline{{v }_{s}}\) (adjusted R2 = 0.899, p < 0.0001; Fig. S1b). Both multiple-unit indices were higher than the average values of the pairwise indices at all times. Both dBC.mu and dsor.mu increased significantly with increasing dMR.mu (adjusted R2 = 0.756, p < 0.0001; Fig. S2a and adjusted R2 = 0.444, p < 0.0001, respectively; Fig. S2b).

The value of dBC.mu was significantly higher in group 4 (swamp) than in the other groups (Fig. 3b,f, Table S2a). In addition, no significant difference was found between groups 1 and 3, while the value of group 2 differed significantly from those of groups 1 and 3 (Fig. 3b,f, Table S2a). The dMR.mu values showed a similar pattern to that of dBC.mu. dMR.mu was significantly higher in group 4 (swamp) than in the other groups (Fig. 3c,g, Table S2b). No clear difference was observed between groups 1 and 3, but the value of group 2 differed significantly from those of groups 1 and 3 (Fig. 3c,g, Table S2b). The vs.mu value was considerably and significantly higher in group 4 (swamp) than in all other habitats (Fig. 3d,h, Table S2c). In addition, no significant difference in vs.mu was found between groups 1 (low plateau) and 2 (high plateau), but the value of group 3 (slope) differed from those of groups 1 and 2 (Fig. 3d,h, Table S2c). There were no significant differences in Pses.all, among groups 4 (swamp), 5 (streamside), and 6 (young forest), but the values were significantly lower in groups 4 and 5 than in the other groups (i.e., groups 1, 2, 3, and 7) (Fig. 3e,i, Table S2d). Moreover, no significant difference was found between groups 1 (low plateau) and 2 (high plateau), but group 3 (slope) differed significantly from groups 1 and 2 in terms of Pses.all (Fig. 3e,i, Table S2d).

In the generalized additive model, compositional variability (dBC.mu) showed a significant response to the degree of deviation in individual persistence from the null model (Pses.all: edf = 3.951, F = 31.59, adjusted R2 = 0.112, p < 0.0001), with the values of the compositional variability index constant from Pses.all of approximately 12–25 and tending to increase from Pses.all of approximately 12–0 (Fig. 4).

Fig. 4
figure 4

The relationship between the degree of deviation in individual persistence from the null model (Pses.all) and the conventional multiple-unit dissimilarity index (dBC.mu). The solid black line was drawn using a generalized additive model fitted using the restricted maximum likelihood method (adjusted R2 = 0.112, p < 0.0001, edf = 3.951, F = 31.59). The blue dashed lines represent the 95% confidence intervals

Full size image

Discussion

Individual turnover is a major driver of compositional change and variability within a community, although it has rarely been discussed. Higher individual turnover clearly facilitates higher variability in community composition; thus, compositional variability over time can occur even when considering only ecological drift. Hubbell’s unified neutral theory of biodiversity (Hubbell 2001) postulates ecological neutrality among species within communities and emphasizes the role of ecological drift in shaping biodiversity patterns. To elucidate compositional variability, the impact of individual turnover (i.e., ecological drift) should be excluded. However, almost no methods exist to quantify the impact of individual turnover on compositional variability. The individual-based indices introduced herein can assess compositional variability considering the contribution of individual turnover (Table 2). In addition, information about individual identity across multiple time steps provides opportunities to test individual persistence. Within and among species, some individuals live longer (i.e., individuals persist in the community for a longer period) than others. The patterns of individual persistence can be tested for their degree of deviation from expected patterns due to stochastic drift using the novel null model introduced here. This model can add a new assessment axis to empirical studies of community dynamics.

Table 2 Summary of indices, formulae, and references citing both the indices developed in previous studies and the novel indices proposed in this study
Full size table

In the present study, I applied both conventional and novel indices, as well as the degree of deviation based on the null model to the BCI plot. The average values of both \(\overline{{d }_{MR}}\) and \(\overline{{v }_{s}}\) across all pairs are always lower than multiple-unit indices (i.e., dMR.mu and vs.mu; Fig. S1), in accordance with the findings of Baselga (2013a). The results of the case study clearly demonstrated that the novel individual-based index was significantly correlated with conventional incidence- and abundance-based indices (Fig. S2); however, variations among communities still exist. The value of vs.mu, which represents compositional variability after considering individual turnover, differed among habitat zones, especially for group 4 (swamp), which differed significantly from all other habitats (Fig. 3a,c, Table S1). Both the degree of compositional variability (dBC.mu) and individual turnover (dMR.mu) were significantly higher in group 4 than in all other habitats. There was significant compositional variability in the swamp even after controlling for individual turnover. Thus, individual turnover is a major driver, but a large component of variability remains unexplained by individual turnover. The environmental stability of each plot is a possible source of the unexplained effect. For example, Nakadai (2020) identified compositional changes after consideration of the individual turnover that was clearly associated with the degree of climate change (i.e., rate of temperature change) across the Japanese Archipelago. Compositional change over time in this swamp habitat has also been reported in terms of responses to a drought event at the BCI plot (Legendre and Condit 2019). In addition, although dBC.mu and dMR.mu were significantly higher in groups 1 (low plateau) and 3 (slope) than in group 2 (high plateau), vs.mu did not follow the same pattern. Instead, vs.mu was significantly lower in groups 1 (low plateau) and 2 (high plateau) than in group 3 (slope). These results suggest that the high plateau (group 2) has lower individual turnover, probably because it is isolated and stable and consequently has lower compositional variability. However, when considering individual turnover, the slope (group 3) differs from both the low and high plateau (groups 1 and 2).

The degree of deviation in individual persistence was lower in group 4 (swamp) than in groups 1, 2, 3, and 7, but not group 5 (streamside) or 6 (young forest), suggesting that the three habitat types (i.e., swamp, streamside, and young forest) may be approaching their expected mortal patterns due to ecological drift. Moreover, based on the generalized additive model, compositional variability was associated with the degree of deviation in individual persistence from the null model (Fig. 4). This suggests that an increase in the proportion of the long-lived individuals (i.e., higher values of Pses.all) reduces compositional variability, even if the level of individual turnover remains the same. Both the results and the interpretation thereof indicate that the novel approach introduced in the present study affords new insights into the forest dynamics of community composition. In previous studies, individual persistence was recognized in community dynamics as a factor facilitating stable coexistence among species and buffering population growth (i.e., the storage effect) (Chesson and Warner 1981; Chesson 2000). The effect of buffered population growth is expected to be low compositional variability and turnover of individuals over time throughout the community. This prediction is consistent with the results related to Pses.all. In other words, habitat types with higher Pses.all (except groups 4, 5, and 6) were better stabilized by the storage effect.

The novel approaches introduced in the present study focus on datasets collected via individual tracking. I chose one of the most representative and widely employed datasets (the BCI forest plot; Hubbell and Foster 1983) as an example. It is difficult to apply the new methods to taxa in which individuals are not identified (e.g., aquatic plankton communities) unless they are all experimentally tracked. However, I emphasize that my approach can be applied in many other scenarios. For example, the Forest Global Earth Observatory (ForestGEO; https://forestgeo.si.edu/) is a global network of scientists and forest research sites including the BCI plot; data from 72 sites across 27 countries are available (Davis et al. 2021). On the regional scale, the Monitoring Sites 1,000 Project launched by the Ministry of the Environment, Japan (Ishihara et al. 2011) includes 60 forest plots across the Japanese archipelago (Ishihara et al., 2011; Nakadai 2020). Both the ForestGEO and the monitoring 1000 projects feature individual tracking; my methods are applicable. In addition, recent developments in remote sensing techniques have rendered it possible to identify individual trees (e.g., Guillén-Escribà et al. 2021; Weinstein et al. 2021), although not (yet) to track individuals. In the near future, remote sensing will yield individual tracking datasets on a global scale. It is sometimes difficult to identify clonal plant species or individuals with multiple stems. However, if the target shifts from the individual to the stem, my methods are applicable. Therefore, depending on the purpose and status of a study, it may be possible to change the analytical resolution and the target. Finally, the concept of individual-based beta-diversity is applicable not only to plant community datasets but also to animal datasets (Nakadai 2020). Specifically, datasets on bird ringing, other classical mark–recapture procedures, and biologging, which record animal ranges and movements, can be analyzed using my methods because multiple time-series datasets collected from the same place can track individual turnover across time and space (Nakadai 2020).

Although temporal beta-diversity has attracted much attention, especially at large spatial scales (Gotelli et al. 2017; Brice et al. 2019; Magurran et al. 2019), it is simply an extension of spatial beta-diversity. Both methodological and conceptual frameworks rooted in community dynamics should be used. Above all, information from “individuals” has been absent from empirical studies of temporal dynamics in community ecology and biodiversity. On the other hand, recent advances in modern coexistence theory (Chesson 2000) suggest that the impact of individuals across generations (e.g., individual persistence of adult trees) is a storage effect that facilitates species coexistence and stabilizes community dynamics. In this context, the methods proposed here are expected to be useful for addressing a wide range of research questions related to temporal changes in biodiversity, especially studies using forest monitoring data to track individuals, by bridging gaps between recent theoretical advances and accumulated empirical data.