A temporal beta‐diversity index to identify sites that have changed in exceptional ways in space–time surveys

Abstract Aim This paper presents the statistical bases for temporal beta‐diversity analysis, a method to study changes in community composition through time from repeated surveys at several sites. Surveys of that type are presently done by ecologists around the world. A temporal beta‐diversity Index (TBI) is computed for each site, measuring the change in species composition between the first (T1) and second surveys (T2). TBI indices can be decomposed into losses and gains; they can also be tested for significance, allowing one to identify the sites that have changed in composition in exceptional ways. This method will be of value to identify exceptional sites in space–time surveys carried out to study anthropogenic impacts, including climate change. Innovation The null hypothesis of the TBI test is that a species assemblage is not exceptionally different between T1 and T2, compared to assemblages that could have been observed at this site at T1 and T2 under conditions corresponding to H0. Tests of significance of coefficients in a dissimilarity matrix are usually not possible because the values in the matrix are interrelated. Here, however, the dissimilarity between T1 and T2 for a site is computed with different data from the dissimilarities used for the T1–T2 comparison at other sites. It is thus possible to compute a valid test of significance in that case. In addition, the paper shows how TBI dissimilarities can be decomposed into loss and gain components (of species, or abundances‐per‐species) and how a B–C plot can be produced from these components, which informs users about the processes of biodiversity losses and gains through time in space–time survey data. Main conclusion Three applications of the method to different ecological communities are presented. This method is applicable worldwide to all types of communities, marine, and terrestrial. R software is available implementing the method.


| INTRODUC TI ON
Community ecology is the scientific study of the interactions among species in natural communities, their distribution through space and their evolution through time, and of the relationships between the species and their environment. Changes in community composition through time are at the center of community ecology research (McEwan, Dyer, & Pederson, 2011;Pickett, Collins, & Armesto, 1987;Vellend, 2016), including the nature of these changes (e.g., gains and losses of species) and their quantitative importance.
In the spatial context, the variation in community composition among sites in a region of interest has been called beta diversity by Whittaker (1972), who defined the well-known concepts of alpha, beta, and gamma diversities. In recent years, interest of ecologists and managers has turned to the study of the temporal variation in community composition, either at a single site or at a series of sites repeatedly surveyed across time. This temporal variation was called temporal beta diversity by Legendre and Gauthier (2014) and Shimadzu, Dornelas, and Magurran (2015). Temporal variation can be the result of gradual or abrupt changes in environmental conditions, including man-induced alterations such as the present worldwide climate warming.
Statistical inference methods have been proposed for the analysis of temporal changes in community composition. For example, (a) the temporal convergence or divergence in composition of a set of communities can be studied by testing for differences in multivariate dispersion among surveys (Anderson, 2006); (b) shifts in mean composition of monitored communities can be statistically tested using multivariate analysis of variance procedures (Anderson, 2001;Legendre & Anderson, 1999), including null models (Schaefer, Gido, & Smith, 2005); (c) the interaction between the factors space and time and other complex spatio-temporal structures can be studied and tested for significance (Angeler, Viedma, & Moreno, 2009;Legendre, Cáceres, & Borcard, 2010;Legendre & Gauthier, 2014).
In several application fields, researchers want to compare observations made at several sites and at two different times. The question of interest is as follows: are there sites where the difference between survey times seems exceptionally large? These sites would be worth examining in more detail to identify and compare the causes of the differences. Here are some examples. (a) In paleoecology, comparison of ancient and modern diatom communities preserved in lake sediment cores may indicate areas where acute anthropogenic processes have singularly changed the surrounding land use (e.g., Winegardner, Legendre, Beisner, & Gregory-Eaves, 2017).
(b) When a strong natural or man-made or environmental impact has taken place at a known point in time and an ecological community had been surveyed ahead of the impact, ecologists may survey that community again to determine how it was affected by the impact, and then how it may have recovered in later surveys (e.g., Legendre & Salvat, 2015). (c) In community ecology, when studying a permanent stem-mapped forest dynamics plot divided into quadrats, examining surveys made at two different times may indicate sections of the forest that have been exceptionally affected by a disturbance, for example, a climatic or anthropogenic event (e.g., Legendre & Condit, 2019). This paper describes a method to test, for several sampling units (objects), the differences between data vectors corresponding to observations made at time 1 (abbreviated T1) and time 2 (T2). I will refer to these objects as sites in this paper, although they may be of other natures, for example, experimental enclosures. A dissimilarity D computed between times T1 and T2 for a site, using community composition (occurrences, frequencies or biomass) or gene frequency data, is called a temporal beta-diversity Index (TBI); it measures the change in community composition (or temporal beta diversity) from T1 to T2. A change through time is directional; species presences, species abundances, or gene frequencies may have been gained and/or lost between T1 and T2. So, it will be of interest to examine the loss and gain components of the TBI indices, in addition to the TBI index values and their significance.
The observed data, assembled in matrices Mat.1 for time T1 and Mat.2 for T2 (Figure 1), may be of different kinds; in landscape ecology and genetics, the data are community composition or population gene frequencies observed at different sites, the same in the two surveys. (a) The null hypothesis (H 0 ) to be tested in the statistical F I G U R E 1 Schematic representation of the first step of the method. Data in matrices Mat.1 (for Time 1) and Mat.2 (for Time 2) are used to compute a vector of TBI dissimilarities D i for all sites i (data rows). For example, for site i = 1, vectors T1 1 and T2 1 are compared to compute D 1 , the dissimilarity between data at time 1 (T1) and time 2 (T2) testing part of the method is that a site is not exceptionally different in community composition between T1 and T2, for presenceabsence or abundance data, compared to assemblages that could have been observed at this site at T1 and T2 under conditions corresponding to H 0 . An exceptionally different site is a site with an index of dissimilarity between T1 and T2 that has an extreme value in the upper tail of the distribution of TBI index values. This value may also have been produced by a different (e.g., ecological) process than the one having generated most other values in the distribution. To determine how extreme a value has to be in the reference distribution before it is considered extreme, we will rely on the usual significance levels (e.g., 5%, 1%, etc.) of statistical analysis. (b) In the data representation part of the method, we will use the species losses and gains between T1 and T2 (for presence-absence or abundance data) to uncover the main ecological changes that have taken place between T1 and T2 in the study area or in subsets (geographic or environmental) of that area.
The paper is organized into sections as follows: Methods (Section 2), Numerical simulations (Section 3), TBI analysis of environmental or trait data (Section 4), Ecological applications (Section 5) where the method is illustrated by three case studies, and Concluding remarks and prospective applications (Section 6).

| ME THODS
The proposed method consists basically in the following analyses of temporal variation, detailed in the following subsections. (2.1) Computation and testing: (2.1.1) a temporal beta-diversity Index (TBI i = D i ) is computed for each site i between the data vectors corresponding to T1 and T2 (Figure 1), using an appropriate dissimilarity index (D). Then (2.1.2), when it is pertinent to the problem at hand, the indices can be tested for significance using a permutational procedure. (2.2) Partitioning dissimilarities into losses and gains: four of the D indices that can be used in this type of analysis (two indices for abundance and the corresponding two indices for presence-absence data) also allow the computation of species losses and gains at each site between T1 and T2. These statistical elements provide users with detailed information, at the site level, about the response of the community to the event or change that occurred between T1 and T2. (2.3) The loss and gain statistics can be used together to draw B-C plots, which illustrate whether the temporal changes at the various sites are dominated by gains or by losses. (2.4) The Software subsection lists the R software available to carry out the analyses. Since Whittaker (1972), ecological dissimilarities have been used to measure beta diversity among sampling units. Koleff, Gaston, and Lennon (2003) reviewed 24 beta-diversity indices proposed in the literature while Legendre and De Cáceres (2013) described 14

| Temporal beta-diversity Indices: computation and testing
properties of 11 dissimilarity indices that are appropriate for betadiversity studies. The Ružička dissimilarity (1958) was identified by Podani, Ricotta, and Schmera (2013) and by Legendre (2014) as another appropriate index for beta-diversity studies.
Some of these indices are used as TBI indices in this paper: the percentage difference (Odum, 1950;sometimes  in the present section of the paper; this index is also known as the Bray-Curtis dissimilarity; see Legendre & De Cáceres, 2013,

| Computation of temporal betadiversity Indices
Consider two data matrices, Mat.1 and Mat.2, about the same objects observed at times T1 and T2; each matrix has n sites (with in- The percentage difference dissimilarity (D %diff ; method "%difference" in the R function TBI.R, also known as the Bray-Curtis index in other computer packages), and the Ružička dissimilarity (D Ruz ; method "ruzicka" in the R function; see Software in Section 2.4 below) can be used for beta-diversity assessment. They are obtained by computing a dissimilarity function (equations shown below). With presence-absence data, the percentage difference produces the (1 -S Sørensen ) dissimilarity whereas the Ružička dissimilarity produces (1 -S Jaccard ), where S designates a similarity index.
The chord, Hellinger, and log-chord distances are members of the Box-Cox family of distances (Legendre & Borcard, 2018). They are classical indices for beta-diversity studies (Legendre & De Cáceres, 2013). Their calculation involves two steps: first the calculation of a transformation of each data row (i.e., the chord transformation, the Hellinger transformation, or the transformation log(y + 1) followed by the chord transformation; Legendre & Borcard, 2018), followed by calculation of the Euclidean distance. These indices, as well as the Euclidean distance itself, are also implemented in the TBI.R function and will be used in the simulations and in Ecological application 1 below, although they are less interesting for the comparison of community composition matrices than the percentage difference and Ružička indices, which provide additional information about losses and gains of species. The indices in the Box-Cox family are available in the computer function for TBI analysis to ensure compatibility with other multivariate analyses that users may want to do using these popular indices.
When the percentage difference or the Ružička dissimilarity is used as TBI indices, one can compute two derived indices to study the directionality of the change through time at each site, as proposed by Legendre and Salvat (2015). Consider data vectors y 1 and y 2 corresponding to the multi-species observations at T1 and T2 for a focal site of interest. The following calculations can be done on these vectors: • A j is the part of the abundance of species j that is common to the two survey vectors: A j = min(y 1j , y 2j ). A is the sum of the A j values for all species. It represents the unscaled similarity between two surveys.
• B j is the part of the abundance of species j that is higher in survey 1 than in survey 2: B j = (y 1j -y 2j ) if y 1j > y 2j ; B j = 0 otherwise. B is the sum of the B j values for all species. It is the unscaled sum of species losses between T1 and T2.
• C j is the part of the abundance of species j that is higher in survey 2 than in survey 1: What are the ecological applications of D loss and D gain ? For each site, one can explore which process, between D loss and D gain , shows the largest contribution to the temporal D %diff or D Ruz dissimilarity; in other words, which process is dominant at each site. The means of the D loss and D gain components across the sites express the dynamics of the community over all sites. For observations across a large number of sites within a region, or (e.g.) in all quadrats of a stem-mapped dynamics forest plot, the B/den and C/den statistics can be mapped, plotted as species losses and gains in the new B-C plots described in Section 2.3 (see example in Ecological application 3), or studied in other ways (e.g., in Ecological application 2).

| Testing procedure
When it is pertinent to the problem at hand, TBI indices can be tested for significance. An example is found in Ecological application 1. The data are permuted at random in both matrices, as described below, and the index is recomputed. This procedure is repeated a large number of times, and a p-value is calculated for the TBI difference between T1 and T2 at each site i. A detailed description of the permutation method follows. are also indices about individual sites.

3.
The data are permuted columnwise, species by species, and in exactly the same way in matrices T1 and T2 ( Figure 2). To accomplish that in a computer function, a given permutation of the two matrices is started with the same random seed in both, and that seed is changed at the beginning of each new permutation. With this method, the abundance values of a given species at site i (e.g., site 1) in matrices T1 and T2 of the original data will be shifted to another site position (e.g., site 9) in both matrices in the permuted data. What is permuted is then a series of differences between T1 and T2, for each species separately.

4.
In the TBI test, we are looking for site vectors whose dissimilarities between T1 and T2 would be exceptionally large. We are not interested in a systematic difference that would affect all sites F I G U R E 2 Random permutation of the community composition data is done separately for each species (column), in matrices T1 (left) and T2 (right). This figure shows an example where a permutation of species j brings value y 1j to position y 9j and value y 5j to position y 1j . The exact same permutation, involving all values in column j, is done in matrices T1 (left) and T2 (right to obtain a correct experiment-wise error rate.
Note 1: Ecological interpretation of the permutation models-The permutation procedure described in paragraphs 2 and 3 above, whereby the values are permuted separately for each species, follows a concept of species assemblages that considers the species in an assemblage to be under the influence of a variety of processes; they do not form a pseudo-organismic entity that would react as a unit to these processes. The first formulation of this concept is due to Gleason (1926) (Hubbell, 2001), and interactions among species, all of which generate spatial structures in community data (Legendre, 1993;Legendre & Legendre, 2012 That permutation method was used in additional simulations.

Results of these simulations are not reported in detail in Supporting
Information Appendix S1 but they are summarized in the second paragraph of section "b. Permutation methods" in that Appendix.
Note 3-It is also possible to identify the species that have changed in a significant way between T1 and T2, either in the whole study area or in separate habitat types.

| Species losses and gains: the B-C plot
We can also use the B/den and C/den statistics as coordinates of points (representing sites) in bivariate graphs with B/den in the abscissa and C/den in the ordinate. We call these graphs B-C plots.
They display visually the relative importance of the loss and gain processes across the study sites, informing researchers about the detailed and global structure of the species losses and gains.
A B-C plot is presented in Ecological application 3 (Chesapeake Bay benthos data). In that B-C plot, a diagonal green line, with slope of 1, was drawn through the origin; it represents the the- In B-C plots, the points representing sites can be labeled with colors or symbols representing the types of environment, the geographic areas where they come from, or any other independent classifier of interest. Separate B-C plots can be drawn for sites surveyed in different types of environment, although all sites may have been analyzed in the same TBI analysis. Comparison of these separate plots will immediately show which types of environment have produced mostly losses or gains in species occurrences or abundances.

| Software
These calculations are implemented in the TBI.R function available in the R package adespatial (Dray et al., 2019). Function plot.TBI.R is also available in adespatial to draw B-C plots. Examples of output files of the TBI function are shown in Appendices S3 and S4, which complement Ecological applications 1 and 3. Function tpaired.
krandtest.R can be used to identify the species whose abundances have changed in a significant way between T1 and T2. That function is described in Supporting Information Appendix S5; it is also available in the R package adespatial on CRAN.
Another R package, codyn (Hallett et al., 2018(Hallett et al., , 2016, computes a variety of metrics for temporal diversity analysis. The percentage difference D is used as one of the metrics for temporal comparison of communities observed at different times at a focal site. Function turnover of that package includes an option to compute species losses and gains divided by the %difference denominator, as in the TBI.R function; losses and gains are the indices called B and C in the present paper.

| NUMERIC AL S IMUL ATIONS
Numerical simulations were used to check the type I error rate and power of the permutation method described in Section 2.1.2 above.
The data simulation methods and results are described in detail in Supporting Information Appendix S1. A summary of these results is presented here, with recommendations to users.

| Simulation to estimate type I error rates
The simulation results reported in Supporting Information Appendix S1 show that the TBI tests had correct rates of type I error for the two community-like data generation methods (Poisson and lognormal) and all dissimilarity indices available in the TBI.R function, and this for all significance levels (α) considered, from α = 0.01 to α = 0.50. The testing method is thus valid in all these circumstances (Edgington, 1995).

| Simulations to compare power of D indices
For the analysis of community composition data, the percentage difference and Ružička indices produced tests with the highest power, followed by the indices in the Box-Cox family: the chord, Hellinger, and log-chord distances. The Euclidean distance alone produced TBI tests with extremely low power with community composition-like data. This distance should not be used for TBI tests of community composition data (Figures S1.5 and S1.6, Supporting Information Appendix S1). However, power was high enough to recommend the Euclidean distance for test involving standardized environmental data.
In summary, the best combination to obtain TBI tests of community data with maximum power is to use the percentage difference or the Ružička indices. These two dissimilarities can also be decomposed into species losses (B/den) and gains (C/den), which can be used to examine the processes of losses and gains at the site level and produce B-C plots.
Additional simulations involving different numbers of sites with an effect and different total numbers of sites showed that power of the test was high as long as the proportion of sites with an effect to be detected (i.e., sites made to be exceptional) was smaller than n/2, where n is the total number of sites in the study ( Figure S1.7, Supporting Information Appendix S1).

| TB I ANALYS IS OF ENVIRONMENTAL OR TR AIT DATA
It could be interesting to identify the sites where the changes in environmental conditions (e.g., land use) were the most important. The TBI method can be used to compare two matrices containing the same environmental variables observed at T1 and T2 and determine, for example, if these sites are also those where the community has changed the most. The analysis of environmental variables is a situation where the Euclidean distance would be appropriate as a basis for computing a TBI index.
• If all environmental variables are quantitative, they should be standardized using the same parameters (means, standard deviations) for matrices T1 and T2, before they are used as input in TBI analysis. How to do that is described and implemented in a function provided in Supporting Information Appendix S2.
• If the data are factors or a mix of quantitative and factor variables, one should join matrices T1 and T2 one on top of the other, as described in the explanation paragraph of Supporting Information Appendix S2, and then compute a Gower dissimilarity matrix D.
Apply principal coordinate analysis (PCoA) to the square-rooted Gower dissimilarities because a Gower D matrix is non-Euclidean.
Square rooting should make the Gower matrix Euclidean, which is necessary before PCoA in this case; we have to recuperate and use all PCoA axes for TBI analysis; hence, the values should be real and not complex numbers. Then, separate the two transformed matrices and use them as input into TBI analysis.
• For community trait matrices with mixed precision levels (quantitative and qualitative traits), use the same method as in the previous paragraph: compute a Gower dissimilarity matrix, as recommended by Laliberté and Legendre (2010), then PCoA of the square-rooted dissimilarities; split the principal coordinates into two matrices and compute TBI analysis using the Euclidean distance. No application of TBI analysis to environmental or trait data is presented in this paper to save space.

| ECOLOG IC AL APPLIC ATIONS
The three ecological applications that follow use multivariate community data. They were chosen to illustrate different facets of TBI analysis, not to draw ecological conclusions about these three particular ecosystems. Application 1 illustrates the importance of carrying out TBI analysis using a dissimilarity index designed for the analysis of community composition data; the Euclidean distance produced nonsignificant and uninterpretable results. In application 2, The B and C components of the TBI dissimilarity are used to analyse the effects of a major disturbance (El Niño) on communities; the analysis is complemented with a standard canonical analysis of the community data. Application 3 illustrates the construction and interpretation of a B-C plot.

| Ecological application 1-Insecticide treatments in mesocosms
Observations on the abundances of 178 invertebrate species (macroinvertebrates and zooplankton) subjected to insecticide treatments in aquatic mesocosms (called "ditches") were used by van den Brink and ter Braak (1999) as an application example in their paper describing Principal Response Curves (PRC) analysis. The authors agreed to make the data available to researchers in the CANOCO program documentation and in the R package vegan (Oksanen et al., 2017).
The experiment involved twelve mesocosms, which were surveyed on eleven occasions. Four mesocosms served as controls (dose = 0) and the remaining eight were treated once with the insecticide chlorpyrifos, with dose levels of 0.1, 0.9, 6.0, and 44.0 μg/L in two mesocosms each. The data are log-transformed species abundances, y tr = log e (10y + 1). In their paper, the authors used the logtransformed invertebrate data in PRC analysis; PRC preserved the Euclidean distance among the observations.

Results
On the contrary, the Euclidean distance is known to be inappropriate for such studies (Legendre & Legendre, 2012;Orlóci, 1978) and, indeed, TBI tests based on that distance did not produce significant differences in community composition between surveys #4 and #11 in any of the mesocosms, including M11 and M12.
Detailed analysis of the species losses (B/den) and gains (C/den), obtained from TBI analysis computed with the percentage difference (Supporting Information Appendix S3, first section), showed that in the eight treated mesocosms (called Sites 5 to 12 in Supporting Information Appendix S3), the changes in community composition (abundance data) always consisted of species gains; that is, statistic C/den (gains) was always larger than B/den (losses).
The mean values of B/den and C/den for these eight mesocosms showed that gains (C/den) represented 56% of the dissimilarities, as expected in a study of recovery after an insecticide treatment.
The permutational paired t test showed a highly significant difference (p = 0.0066) between losses and gains across the eight treated mesocosms.

TBI calculations using the Sørensen D (Supporting Information
Appendix S3, Section 2) indicated that, in addition to mesocosms #11 and 12, mesocosm #10, treated with 6 μg/L of insecticide, also displayed a significant difference between T1 and T2 at significance level 0.05. This result indicates that in the insecticide experiment, the reappearance of species (positive change) gave a clearer signal of community recovery than the increase in species abundances (Supporting Information Appendix S3, Section 1). Brown and Suharsono (1990) surveyed coral communities (75 species) at 10 sites in the island of South Tikus, Indonesia, in the years 1981, 1983, 1984, 1985, 1987, and 1988. An El Niño event occurred in 1982-1983, which caused coral bleaching and death of coral colonies, and triggered changes in the composition of coral communities.

| Ecological application 2-South Tikus Island coral communities
They reported that "as many as 80-90% of corals died on the reef F I G U R E 3 Tikus Island coral data. (a) Changes in dissimilarity D computed from the quantitative coral community compositions between years, and its components B/den (losses) and C/den (gains); den is the denominator of the dissimilarity index D, (2A + B + C) in this figure.
The 1981 survey, before the El Niño event, is compared in turn to the 1983, 1984, 1985, 1987, and 1988 surveys. (b) Same for the species occurrence (i.e., presence-absence) data flats at the study sites, with the major casualties being branching species in the genera Acropora and Pocillopora." Coral forms colonies which occupy surfaces, so that the data are not in numbers of individuals but in areal cover of each species. The sum of the species areal covers at a site may exceed 100% because coral colonies may overlap one another vertically. The Brown and Suharsono (1990) data have been used in several papers to demonstrate the application of multivariate methods for the analysis of beta diversity and the comparison of surveys across time, for example, by Warwick, Clarke, and Suharsono, (1990), Anderson et al. (2011), and Chao and Chiu (2016). Following these papers, the data in the present application were treated as if they were species abundances. They were obtained from Appendix S1 of the Anderson et al. (2011) article.
The data were analyzed in two complementary ways. (a) Figure 3 presents an analysis of the species loss (B/den) and gain (C/den) components of the dissimilarity D between the 1981 survey, before the El Niño event, and the five following surveys: 1983, 1984, 1985, 1987, and 1988, for abundance and presence-absence data. (b) This analysis is completed with a redundancy analysis (RDA) biplot shown in Figure 4, showing the changes in community composition with time. This biplot was produced as follows: first, a percentage difference matrix was computed among all years and sites; the dissimilarities were square-rooted to make the matrix Euclidean, and that matrix was subjected to principal coordinate analysis (Gower, 1966).
The entire matrix of principal coordinates was used as the response data in a RDA against a factor representing the six survey years of the study. This form of canonical ordination is called distance-based redundancy analysis (dbRDA, Legendre & Anderson, 1999).
We will examine the changes in community composition between the 1981 survey, before the El Niño event, and the following five surveys: 1983, 1984, 1985, 1987, and 1988. This study is not meant to identify sites that were exceptionally different between two years or test specific hypotheses about them. Instead of testing the TBI statistics for sites, we will carry out a detailed study of the species loss (B/den) and gain (C/den) statistics described in the Methods. These statistics were computed with the denominator (den) of the percentage difference index, D %diff ; they decompose the percentage difference D into additive components, losses (B/den) and gains (C/den).
First, we will plot the mean values of the B/den, C/den, and D statistics computed across the sites, in comparisons with the 1981 survey (before the El Niño event) with all successive surveys in turn (1983, 1984, 1985, 1987, and 1988, after El Niño) (Figure 3  The TBI function was run over all year pairs over the 10 study sites. Results showed that dominance of B/den (losses) over C/den (gains) was significant for all year pairs, as shown by the overall paired t tests of the asymmetry, described in the Section 2.
Does that mean that some of the species that had disappeared had recovered, or that only the species that remained had increased their abundances-per-species? The answer is found in Figure 3b, which displays the same statistics computed for species occurrence data. In that graph, the B/den and C/den lines would be horizontal at value 0 if all species had remained present and the change after El Niño was only in the abundances. Instead, that graph shows that many species disappeared at first from the surveyed sites after El seem to represent random variation due to observed random losses and gains of species, which may be due in part to sampling variation and in part to random species losses and gains.

The communities found in South Tikus Island after the natural El
Niño event strongly differed in species composition from the structure they had in 1981, and they kept changing, apparently randomly, in later years. A similar phenomenon was shown in the Legendre and Salvat (2015) study, where the disturbance of marine mollusk communities was due to a strong man-made disturbance. In both cases, the observed changes are compatible with the neutral theory of generation of biodiversity and changes in communities (Hubbell, 2001).

| Ecological application 3-Chesapeake Bay data
The data set used in this example was extracted from the Maryland Data Sets of the Chesapeake Bay Benthic Monitoring Program A simple classification of the sites by an environmental factor, water temperature during the 2005 fall survey, was used to separate the sites in two groups, providing an example of the kind of information that can be derived from displaying different habitat groups as symbols or colors in B-C plots. These two groups of sites could also be drawn in separate B-C plots. These separate plots would show TA B L E 2 Number of species in subsets of the Chesapeake Bay fauna data surveyed during 13 years, spring and fall. In total, 205 benthic species were found at the 27 survey sites They are collecting data over land, in lakes and in the oceans to assess the effects of climate change on natural communities and other types of biodiversity data. Researchers would like to identify the sites where important changes have taken place. They can then focus their attention onto these sites and seek what has been going on there, and why community composition has changed in an exceptional way at these sites. The TBI method was designed for this type of research.
The present paper is the first description of the ecological theory and statistical developments behind TBI analysis, describing the method in its present state of development. It should provide opportunities to researchers to apply this new method of analysis to a broad palette of ecological and genetic questions.
The paper has shown that it is possible to compute a valid test of significance for dissimilarity indices, which are used to compare data about sites collected at different times. Additionally, the paper has shown how four of the TBI indices can be decomposed into loss and gain components (of species, or abundances-per-species) and how these components can be used to produce B-C plot, a new type of plot that informs users about the processes of biodiversity losses and gains through time found in space-time survey data.
The following indices were found to be appropriate for computation and testing of TBI indices for community composition data: the percentage difference, Ružička, Sørensen, Jaccard, chord, Hellinger, and log-chord indices. The first four present the advantage that they can be decomposed into species losses and gains at each site, which makes them the preferred choices in most studies. In studies where community data must be projected into Euclidean space, for instance before PCA, RDA, or other forms of linear analysis, data can be transformed using the chord, Hellinger, and log-chord transformations, and the TBI analysis can be performed using the corresponding D indices, that is, the chord, Hellinger, or log-chord distances. The simulation study has also shown that the Euclidean distance is inappropriate for TBI analysis of community composition data. It can be used, however, for TBI analysis of environmental or species trait data.
These authors have shown that the Jaccard/Sørensen and Ružička/ percentage difference D indices can be decomposed into replacement (or spatial turnover), 2 min(B,C) and richness/abundance difference, |B-C|, for either presence-absence or abundance data.
Because they developed their indices for spatial analysis, Baselga and Podani did not emphasize the asymmetry implicit in a temporal comparison between T1 and T2. In temporal studies, where processes are directional, comparison of the loss (B) and gain (C) components, which are central to the TBI method described in the present paper, is informative.
Analysis of the B and C components brings us to the heart of the mechanisms by which communities change through time: losses (b) and gains (c) of species, losses (B) and gains (C) of individuals of the various species. B-C analysis is especially interesting in species-rich communities where researchers cannot examine the changes in each species individually.
B-C analysis can also be applied to subgroups of sites, for example, habitat types. In addition, it can be used to compare the changes that occurred in specific groups of species that are known to react differently to environmental stressors, for example, different age or size classes, or species of different origins, for example: the temperate, transitional, and boreal trees found together in the forested southern portion of Canada.
Different ecological applications were worked out with co-authors during the development of the TBI method. Some of them have already been published. Working on these papers provided opportunities to develop the TBI theory and software, through the analysis of pertinent application questions, hypotheses, and data.
These applications provide examples that ecologists may find useful as guides for the analysis of their own data, in addition to the ecological applications summarized in the previous section of the paper: • Impact of a field experiment-The loss (B/den) and gain (C/den) statistics were first analyzed by Legendre and Salvat (2015) to compare community composition data (marine mollusks) during 30 years, before and after a man-made disturbance on an atoll in the Pacific. This disturbance to the mollusk community was the atmospheric test of a Hydrogen bomb in 1968. In a particular study, researchers may be mostly interested in identifying the sites with high and significant TBI indices. In other studies, interest may be in a fine analysis of the changes in the loss and gain components of the dissimilarity in community composition, compared to a pre-disturbance situation. One can look at these components in graphs that allow researchers to compare, for example, different subsets of the data. Ecological examples have been shown in the paper for these different situations.

ACK N OWLED G M ENTS
I am thankful to Amanda Winegardner and Lucie Kuczynski who experimented with this testing procedure and computed it on real data during the development phase of this paper and to Daniel Borcard and Marie-Hélène Brice who provided useful comments on the manuscript before submission. Lucie Kuczynski suggested the acronym TBI for the new method. I am also thankful to Daniel Borcard who revised a preliminary version of the manuscript and to Joseph R.
Bennet and three anonymous reviewers who provided helpful comments on the submitted manuscript. This research was supported by Natural Sciences and Engineering Research Council of Canada (NSERC) research grant no. 7738 to P. Legendre.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R CO NTR I B UTI O N S
P. Legendre designed the methods described in this paper, programmed the software, conducted the analyses of the numerical examples and the simulation study, and wrote the manuscript.

DATA ACCE SS I B I LIT Y
The data used in the ecological applications are publicly available in the references provided.