Spatial and temporal analysis of beta diversity in the Barro Colorado Island forest dynamics plot,Panama

Pierre Legendreand Richard Condit

Abstract Background:Ecologists are interested in assessing the spatial and temporal variation in ecological surveys repeated over time.This paper compares the 1985 and 2015 surveys of the Barro Colorado Forest Dynamics plot(BCI),Panama,divided into 1250(20m×20m)quadrats.Methods,spatial analysis:Total beta diversity was measured as the total variance of the Hellinger-transformed community data throughout the BCI plot.Total beta was partitioned into contributions of individual sites(LCBD indices),which were tested for significance and mapped.Results,spatial analysis:LCBD indices indicated the sites with exceptional community composition.In 1985,they were mostly found in the swamp habitat.In the 2015 survey,none of the swamp quadrats had significant LCBDs.What happened to the tree community in the interval?Methods,temporal analysis:The dissimilarity in community composition in each quadrat was measured between time 1(1985)and time 2(2015).Temporal Beta Indices(TBI)were computed from abundance and presence-absence data and tested for significance.TBI indices can be decomposed into B=species(or abundances-per-species)losses and C=species(or abundances-per-species)gains.B-C plots were produced;they display visually the relative importance of the loss and gain components,through time,across the sites.Results,temporal analysis:In BCI,quadrats with significant TBI indices were found in the swamp area,which is shrinking in importance due to changes to the local climate.A published habitat classification divided the BCI forest plot into six habitat zones.Graphs of the B and C components were produced for each habitat group.Group 4(the swamp)was dominated by species(and abundances-per-species)gains whereas the five other habitat groups were dominated by losses,some groups more than others.Conclusions:We identified the species that had changed the most in abundances in the swamp between T1 and T2.This analysis supported the hypothesis that the swamp is drying out and is invaded by species from the surrounding area.Analysis of the B and C components of temporal beta diversity bring us to the heart of the mechanisms of community change through time:losses(B)and gains(C)of species,losses and gains of individuals of various species.TBI analysis is especially interesting in species-rich communities where we cannot examine the changes in every species individually.

Keywords:Beta diversity,B-C plots,BCI forest dynamics plot,Space-time analysis,Temporal beta diversity

Background

Whittaker(1972)defined beta diversity as spatial differentiation,or the variation in species composition among sites within a region of interest.Different ways of measuring beta have been proposed by Whittaker himself and by other authors.Several authors have argued that the variance of a community composition data table(sites x species)is an appropriate measure of beta diversity across the sites(Koleff et al.2003;Legendre et al.2005;Anderson et al.2011).It has also been shown that the total variance can be computed from a dissimilarity matrix D obtained from an appropriately chosen dissimilarity index D.Legendre&De Cáceres(2013)described properties of dissimilarity indices that are necessary for beta diversity studies and identified 11 indices that were appropriate for such studies;two more indices were added to that list by Legendre(2014)and Legendre&Borcard(2018).

The concept of beta diversity can be extended to time.It was called temporal beta diversity by Legendre and Gauthier(2014)and by Shimadzu et al.(2015).Temporal beta variation can be the result of gradual or abrupt changes in environmental conditions,including man-induced alterations such as the ongoing worldwide climate warming.Research on temporal variation in communities is conducted,for example,at hundreds of research sites affiliated to the International Long Term Ecological Research(ILTER)network(https://www.ilter.network),and in the dynamics forest plots affiliated with the Forest Global Earth Observatory(ForestGEO;see https://forestgeo.si.edu/).

A statistical method for the analysis of temporal beta diversity has recently been described by Legendre(2019).The objective is to compare observations made during two separate surveys through time,involving several sites.The method has two distinct parts:(1)estimate the change in each geographic sampling unit(site)between time 1(abbreviated T1)and time 2(T2),using an appropriate dissimilarity index,called a temporal beta-diversity index(TBI),and test the significance of that change to identify the sites where the change has been exceptionally important;these sites may be worth examining to identify the causes of the differences.And(2)partition the dissimilarity information into finer indices of losses and gains of species,or of abundances-per-species,which may tell us something about the processes at work in the system,which may have generated these changes.Applications of that method have already been published in palaeoecology(Winegardner et al.2017)and in the study of freshwater(Kuczynski et al.2018)and marine animal communities(Legendre&Salvat 2015).Three other ecological applications of the method are presented in the Legendre(2019)paper.

The TBI method will be illustrated in this paper by an analysis of a 50-ha permanent forest dynamics plot.The precise questions we will try to answer about the plot are the following:(1)Is there a region in the plot where the changes in community composition have been especially important between 1985 and 2015?(2)If there is such an area,what happened in the community?Was the change characterized by species(orabundances-per-species)losses or gains?(3)And finally,what are the species that changed more strongly?If the change was mostly a gain of species,was this the result of a migration or an invasion,and where did these species come from?

Methods

Plot and censuses

The Barro Colorado Forest Dynamics plot(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,including a recent one by Condit et al.(2017b)in this journal,where the plot and the survey method are described.Some statistics about the plot are provided in Table 1.This 50-ha plot has been subjected to 8 detailed surveys since its inception:in 1982+,1985,1990,1995,2000,2005,2010 and 2015.The present paper is the first application of the TBI method of analysis to the data of a permanent forest dynamics plot.It will serve as a model for the analysis of other forest plots.At the moment,the ForestGEO network includes some 70 such forest plots throughout the world.The census data can be downloaded from Condit et al.(2012,2017a).

The BCI data were divided into 1250 quadrats of size(20m×20 m)before analysis;this was deemed the most appropriate scale of study to answer the ecologicalquestions addressed in the present paper,which compares the surveys conducted in 1985 and 2015.The union of the species identified in the 1985 and 2015 surveys was 328.Figure 1 shows a habitat classification map of the BCI forest proposed by Harms et al.(2001).

Table 1 Basic information about the Barro Colorado Forest Dynamics plot(BCI)in Panama.For more information,see Condit et al.(2017b)and https://forestgeo.si.edu/sites/neotropics/barro-colorado-island

Statistical methods for spatial analysis

A preliminary spatial analysis of the forest plot will be done first.Total beta diversity,computed as the total variance of the Hellinger-transformed multivariate community composition data throughout the BCI plot,will be partitioned into Local Contributions to Beta Diversity(LCBD indices)in order to determine where the quadrats are found that are the most different from the multivariate centroid of all quadrats.LCBD indices can be tested for significance.The computation method and interpretation have been described by Legendre&De Cáceres(2013).Calculations were done using the beta.div function of the adespatial package in R(Dray et al.2019).

Statistical methods for temporal analysis

The main part of the paper will consist in a temporal beta diversity analysis focussing on the three questions listed at the end of the Introduction.The methods are the following.

Space-time interaction

For censuses repeated through time,as it is the case in Forest Dynamics Plots,the interaction between space(S)and time(T)is of great interest.A significant interaction would indicate that the spatial structure of the response data has changed through time,and conversely,that the temporal variations differed significantly among the sites(quadrats),indicating,for example,an effect of climate change.The space-time interaction(STI)can be tested on univariate or multivariate data using the method proposed by Legendre et al.(2010),despite the lack of replicated observations.That method is implemented in function stimodels()of the adespatial package in R(Dray et al.2019).The community data were transformed using the Hellinger transformation(Legendre&Gallagher 2001)before this analysis.

Temporal beta diversity indices(TBI)

Fig.1 Map of the BCI forest dynamics plot habitat classification proposed by Harms et al.(2001).The plot is divided into 1250 quadrats of size(20 m×20 m)

The dissimilarity in community composition was measured for each of the 1250 quadrats of the BCI plot between time 1(T1=1985 survey)and time 2(T2=2015 survey).The analysis was repeated separately for the 30 quadrats of group 4(the swamp).

The index used for community composition data will be the percentage difference(or%difference;Odum 1950),also known as the Bray-Curtis index.The S?rensen index will be used for occurrence(presence-absence)data;this is the binary form of the percentage difference index(Legendre&Legendre 2012).In the context of a comparison through time,these indices are called Temporal Beta Indices(TBI).These indices can be tested for significance,as shown in Legendre(2019).Each index,which compares data from a quadrat at T1 and T2,is composed of two parts:B=species(or abundances-per-species)losses and C=species(or abundances-per-species)gains.

The B and C statistics will be used to produce B-C plots,with B(losses per quadrat)in the abscissa and C(gains per quadrat)in the ordinate,as described in Legendre(2019).B-C plots display visually the relative importance of the loss and gain processes in a study area,informing researchers about details of the processes of biodiversity losses and gains across the sites,through time,in space-time surveys.B-C plots will be shown for the whole BCI forest plot and for six separate habitats within the plot.

The mean of the differences between the B and C statistics is also computed across all sites in a study.A positive value of(C-B)indicates that the study area was dominated by gains,whereas a negative value indicates overall losses of species or abundances-per-species.The(C-B)difference across all quadrats was tested for significance using a paired t-test computed for the C and B statistics from all sites in the study.The calculation involves the mean of the differences(C/den-B/den)at all sites;in the present study,den was the denominator of Odum's percentage difference index,(2A+B+C).

The calculations are implemented in the TBI()and plot.TBI()functions,available in the R package adespatial(Dray et al.2019).

Changes in species composition from T1 to T2

In order to understand in more detail the demographic changes at BCI in the quadrats of the swamp habitat in the 1985-2015 interval,we analysed the changes in the 209 species found in the swamp from 1985 to 2015,across its 30 quadrats,using paired t-tests.The tests were carried out with 9999 random permutations of the values,in each quadrat,between T1 and T2.A Holm(1979)correction for multiple testing(209 simultaneous tests had been conducted)was applied to the computed p-values.The calculations are implemented in function tpaired.krandtest.Rpaired.t.test.spec(),available in the R package adespatial(Dray et al.2019).

Results

Preliminary spatial analysis

The quadrats that have high and significant LCBD indices(quadrats with red rims in Fig.2a)are the most different in species composition from a hypothetical quadrat found at the multivariate centroid of a PCA ordination.The significant quadrats are found mostly in the central portion of the map,in an area called the swamp by Harms et al.(2001).This is an old forest habitat where water accumulates during the rainy season and which harbours a different flora than the rest of the BCI plot.The oil palm,Elaeis oleifera(Kunth)Cortés,is a conspicuous indicator species of the swamp community(Plate 1).

Comparison of Fig.2a and b shows that no quadrat with a red rim remained in the swamp area in the 2015 survey data.This preliminary analysis allows us to ask a new question:Why is it that quadrats with highly significant LCBD indices are no longer found in the swamp in 2015?

Temporal beta diversity analysis

The number of stems declined substantially during that period 1985-2015(Fig.3):there was a net loss of 35,638 stems between 1990 and 2005,representing about 15%of the stems that were found in 1990(244,015 stems).

Space-time interaction

A significant space-time interaction was reported across 4 censuses(from 1982 to 83 to 1995)by Legendre et al.(2010)for the whole community of 315 species in the(20m×20 m)quadrats and for 43%of the individual species tested separately.

We checked that this was still the case for the seven censuses conducted between 1985 and 2015.Indeed,a significant space-time interaction was detected for the(20m×20 m)quadrats.This significantinteraction means that the spatial structure of the multivariate data has changed significantly among the censuses.We will now study these differences in more detail between censuses#2(1985)and#8(2015)for the whole plot,then for each habitat type separately,using temporal beta diversity indices and B-C plots.

Temporal beta diversity indices(TBI)

The results for the group 4 quadrats,abundance data,are mapped in Fig.4;the TBI()function output is shown in Supporting information,Additional file 1:Appendix S1.Two quadrats had significant TBI indices at the 0.05 significance level after correction for 30 simultaneous tests(underscored values in output vector$p.adj,Additional file 1:Appendix S1);3 more swamp quadrats were significant at the same level before correction for multiple tests.These 5 swamp quadrats had the highest TBI values among the 1250 quadrats of the BCI forest plot,indicating that they were the most different in composition between the 1985 and 2015 surveys.

Fig.2 The BCI plot divided into 1250 quadrats.The LCBD index values are represented in shades of grey;darker is higher.The quadrats with indices significant at the 0.05 significance level have a blue rim.After Holm correction for multiple testing(1250 simultaneous tests),the quadrats that were still significant at the 0.05 level are represented with a red rim.a 1985 survey;b 2015 survey

Plate 1 The swamp near the centre of the BCI forest plot.The swamp is seasonal;it is inundated during the wet season.Stepping stones installed by researchers are visible in the picture.The oil palm,Elaeis oleifera,is a conspicuous indicator species of the swamp community.Photo P.Legendre,4 February 2015

The mean of the differences between of C/den(losses)and B/den(gains)over the entire BCI forest plot was negative,indicating dominance of losses across all 1250 quadrats(detailed output file not shown).The test of significance of that difference was highly significant.In the separate analysis of the group 4 data(detailed output file in Additional file 1:Appendix S1),this same comparison showed a positive difference of the mean of(C/den-B/den)over the swamp quadrats,indicating thatgains ofindividuals-per-species had dominated group 4 in the time interval under study;see output element$t.test_B.C in Additional file 1:Appendix S1.The test of significance of that difference was also highly significant.

B-C plots

A first B-C plot is shown in Fig.5 for the 1250 quadrats of the BCI forest.In this overall plot,the green line is above the red line,materializing in a figure the statistical resultfound in theprevioussubsection:lossesof abundances-per-species dominated the changes in the entire BCI forest from 1985 to 2015.

Separate B-C plots were then produced for habitat groups 1 to 6,abundance data(Fig.6).In all habitat types except group 4,the green line is above the red line,indicating that species losses dominated gains.In group 4(the swamp),however,the red line is above the green line,indicating that gains dominated losses in this habitat.The same result was found when analysing occurrence data in the 6 habitat types;see Supporting information,Additional file 1:Appendix S2.

Changes in species composition from T1 to T2

The paired t-tests computed separately for the 209 species of group 4(the swamp)over its 30 quadrats showed that 7 species had significant increases in abundances(number of stems)in these quadrats from 1985 to 2015(Table 2a)whereas one,Bactris major,decreased in abundance.The species that increased in abundances also increased in occurrences,and conversely.

Fig.3 Number of tree stems in the BCI forest plot recorded during 7 surveys from 1985 to 2015

Fig.4 Map of TBI indices,group 4(swamp).Circle radii are proportional to the TBI index values.Circle fill colours:gold,C＞B(gains);white,B=C;grey,B＞C(losses).Rim colours:red,adjusted p-value≤0.05;green,unadjusted p-value≤0.05 but adjusted p-value＞0.05;black,p-value not significant at 0.05 level

For comparison,six of the seven species that increased in group 4 also increased in abundances in group 1(Table 2b).Only Trichilia tuberculata decreased in group 1 while increasing in group 4.

Of course,these 209 individual tests are not proper paired t-tests of comparison of means because the data are from contiguous quadrats,hence they are strongly autocorrelated.Yet,the tests allowed us to identify the species where the differences in abundances between T1 and T2 were the strongest in group 4(Table 2a).B.major decreased in abundance and in occurrences in both habitats.

The species we were looking for here are those that changed the most in the T1-T2 interval.This is different from the search for indicator species of the six habitat groups.For example,we mentioned above that the oil palm,Elaeis oleifera,was a most conspicuous indicator species of the swamp habitat.This species has not noticeably changed in abundance in the swamp between 1985 and 2015,going from 16 to 17 individuals.

Discussion and conclusion

Results of the individual analyses presented in the previous section combine into a coherent story of what happened in the BCI forest plot.

Preliminary spatial analysis

Analysis of the local contributions to beta diversity(LCBD)at the quadrat level showed that in 1985,the swamp area(habitat group 4)had a highly distinctive community,with many quadrats significantly different from the other habitat types of the BCI plot.In 2015,this community was still different in composition but the importance of the difference had subsided.This observation called for a more detailed analysis of the difference between the 1985 and 2015 surveys in the swamp area.

Temporal beta diversity analysis

At BCI,the period 1982-1992 included several extreme dry seasons;since then,there have been few such drought events.High tree growth rates and death rates were observed during the drought but not since then(Condit et al.,2017b).A simple indicative statistic is the loss of 36,761 tree stems between 1990 and 2010,or 15%of the stems surveyed in 1990(Fig.3).

Fig.5 B-C plot for species abundance data,1250 quadrats(dots)of the BCI forest.Green line with slope of 1:line where gains equal losses.The red line is drawn parallel to the green line(i.e.with a slope of 1)and passing through the centroid of the points.Its position below the green line indicates that,on average,species losses dominated gains between 1985 and 2015.Dot fill colours correspond to habitat types;habitat colour codes are the same as in Fig.1.Red arrow:direction where D=(B+C)increases

The annual amount of rain has a strong influence on tropical forest composition.It determines the height of the water table and has consequences on many other environmental factors.Changes in soil wetness and dryness is an important factor determining changes in species composition in the long run,considering the fact that tropical trees are long-lived.However,one cannot predict precisely which species will change their distributions following changes in rainfall.Different species react differently to the drought.Then,neutral theory(Hubbell 2001)predicts that changes in species composition have a strong random component,with the constraint that the species adapted to greater wetness are more likely to expand if rainfall increases whereas those more adapted to dry conditions are more likely to regress.This question was discussed in depth by Condit et al.(1996).Legendre et al.,2010,Fig.4)illustrated the changes to two species at BCI,Poulsenia armata and Beilschmiedia pendula,between 1983 and 1995.The first one decreased in most quadrats where it was present whereas the second increased in abundance.Chisholm et al.(2014)analysed changes in tree abundance in 12 forests across the world over periods of 6-28years and pointed out the importance of the variance in environmental drivers in determining the changes in community composition.

A significant space-time interaction was found in the analysis of the seven censuses conducted at BCI between 1985 and 2015,meaning that the spatial structure of the multi-species community has changed significantly among censuses,confirming statistically the observation of the field scientists.

Among the Temporal beta diversity indices(TBI)computed over all 1250 quadrats of the BCI forest plot,the 5 quadrats that had the highest indices were found in the swamp habitat.What is so special in this habitat?It seems to have been more strongly affected by the drought of the 1980's and 1990's than the other habitats.

Fig.6 B-C plots for species abundance data,six habitat types.Symbol diameters are proportional to the TBI statistics.In all habitat types except group 4,the green line is above the red line,indicating that species losses dominated gains from 1985 to 2015.In group 4(the swamp)on the contrary,the red line is above the green line,indicating that gains dominated losses in this habitat.See Fig.5 where the construction of the green and red lines is described.Symbol fill colours correspond to habitat types;habitat colour codes are the same as in Fig.1

From 1985 to 2015,35,697 tree stems(or 14.9%of the 1985 survey counts)were lost in the BCI forest plot excluding the swamp,whereas in habitat group 4(the swamp),1354 stems were gained(or 48.4%compared to the 1985 survey).The tree community in the swamp reacted very differently from the rest of the BCI plot to the drought this region of Panama suffered from 1982 to 1992.This is a major difference of the swamp community compared to the remainder of the BCI forest plot.The B-C plots in Fig.6 show these losses and gains in detail in the various habitat groups.

Of the 1354 tree stems that were gained in the swamp,595 belonged to the 7 species that showed significant increases in that habitat(Table 2a).Most of these species had also increased in the surrounding habitat group 1(Table 2b).The drought created gaps in the forest canopy;the few species that were favoured by the drought moved into theswamp and modified itsspecies composition.

B-C analysis:What have we learned?

Analysis of the B and C components of the TBI dissimilarity between surveys has brought us to the heart of the mechanisms by which communities change through time:losses(b)and gains(c)of species,losses(B)andgains(C)of individuals of the various species.B-C analysis is especially interesting in species-rich communities like the BCI forest,where we cannot examine the changes in each species individually.B-C analysis can be applied to subgroups of sites,as was done in this study for the 6 main habitat types of the BCI dynamics forest plot.It provided detailed information on difference in dynamics between habitat types.

Table 2(a)The 8 species that changed significantly in abundance in the 30 quadrats of group 4(the swamp)from T1(1985)to T2(2015).(b)Changes in abundances of these 8 species in 500 quadrats of group 1 during the same time interval.A taxonomic data table archived in Condit et al.(2017a)includes the family names and taxonomic authorities for all 455 species of the BCI data base

With this paper,we are hoping to encourage ecologists analysing other dynamics forest plots to use TBI analysis to answer questions about the changes that have taken place in their forests,either after gradual changes in environmental(e.g.climatic)conditions or as the result of punctual events that may have affected the plots,e.g.hurricanes or ice storms.

Additional file

Additional file 1:Appendix S1.TBIANALYSIS OF GROUP4(THE SWAMP)USING THETBI()FUNCTION INR.Appendix S2.B-CPLOTS,OCCURRENCE(PRESENCE-ABSENCE)DATA.

Acknowledgements

We acknowledge Robin Foster as the lead field scientist setting up the Barro Colorado 50-ha plot,whose botanical expertise was essential over the first three censuses and whose guidance was essential for all of us still working in the plot.The censuses have been made possible through support of the U.S.National Science Foundation(awards 8206992,8906869,9405933,9909947,0948585 to S.P.Hubbell),the John D.and Catherine D.McArthur Foundation,and the Smithsonian Tropical Research Institute.We also thank the dozens of field assistants and botanists who have measured and mapped the trees over the past 35years.This research was also supported by research grant#7738 from the Natural Sciences and Engineering Research Council of Canada(NSERC)to P.Legendre.

Funding

Statements about funding are included in the Acknowledgements paragraph.

Availability of data and materials

The BCI census data can be downloaded from Condit et al.(2012,2017a).See the References section of the paper.

Authors'contributions

PL developed the methods of analysis,wrote the software,carried out the analyses of the BCI forest data and wrote most of the paper.RC directed the field work at BCI since 1991,organized all aspects of the database,helped with the interpretation of the results,and wrote portions of the Discussion in the paper.Both authors read and approved the final manuscript.

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Author details

1Département de sciences biologiques,Université de Montréal,C.P.6128,succursale Centre-ville,Montréal,Québec H3C 3J7,Canada.2Field Museum of Natural History,1400 S.Lake Shore Dr,Chicago,IL 60605,USA.3Morton Arboretum,4100 Illinois Rte.53,Lisle,IL 60532,USA.

Received:29 September 2018 Accepted:30 January 2019