In karinorman/biodivTS: What the Package Does (One Line, Title Case)

library(kableExtra)
library(magrittr)
library(tidyverse)
#library(printr)
library(rticles)
library(knitr)

knitr::opts_chunk$set(cache=FALSE, message=FALSE, warning=FALSE, out.width = "\\textwidth")
options(knitr.kable.NA = '')

# no sci notation integers pleases
options(scipen=999)

Significance Statement {.unnumbered}

How are ecological communities responding to anthropogenic change? Functional trait approaches allow us to assess this questions by directly measuring shifts in the characteristics that define a species' ecological role. Here, we pair vertebrate community time series with functional trait data to take a first look at temporal trends in functional composition across biomes, realms, and taxa. We found no evidence of systematic functional loss across communities, with maintenance of all measured aspects of functional structure. We also highlight the need for targeted data collection and methodological expansion to further assess functional trends.

Introduction

Ecological communities are experiencing unprecedented change as a result of anthropogenic pressures such as climate change, land use change, and invasive species. Impacts of these pressures are well documented at a global scale by an accelerating global extinction rate [@barnosky2011], and fundamental changes in some of the most well-studied systems [e.g. coral bleaching, @sully2019]. At the local scale however, species diversity tells a different story. Recent syntheses of local trends in biodiversity over time have found no net change in local species diversity despite ongoing turnover [@brown2001; @dornelas2014; @vellend2013; @vellend2017] and evidence of significant shifts in community composition underlying consistent species richness [@brose2016; @gotelli2017; @li2020]. While communities are clearly changing, our most common species-based approaches do not fully capture the nature of that change.

Functional diversity offers a potentially powerful alternative for detecting and describing community change by providing a mechanistic link between species' response to environmental change (response traits) and the processes they perform (effect traits) [@lavorel2002; @mcgill2006; @suding2008]. By describing the functional trait space, functional diversity metrics capture the disproportionate impact of losses or gains of functionally unique species. Functional diversity metrics are therefore particularly well suited for assessing community shifts underlying even constant species richness trends.

Beyond simply characterizing changes in community structure, trends in functional composition also have important implications for ecosystem stability, function, and resilience. There is increasing evidence functional diversity is a better predictor of ecosystem function than species-based metrics [@cadotte2011; @gagic2015], and that different facets of functional diversity play essential roles in maintaining ecosystem stability [@craven2018; @morin2014]. Indeed, almost all hypothesized mechanisms underpinning the relationship between species diversity and ecosystem function are trait-dependent [@hillebrand2009]. Determining functional trends therefore gives a more fundamental picture of potential trends in critical ecosystem processes.

It is critical to establish whether or not functional loss is prevalent across communities. While functional loss is frequently cited as one of the most pressing concerns of the anthropocene [@cardinale2012; @dirzo2014; @young2016], it is not necessarily inevitable even in scenarios of species loss [@daz2001]. Forecasts of functional loss range from negligible [@gallagher2013] to dire [@petchey2002; @pimiento2020]. And while some observed trends show significant functional loss [@flynn2009] others document no loss even in some of the most heavily impacted communities [@edwards2013; @matuoka2020]. On paleoecological time scales functional composition shows mixed responses to environmental change and extinction events [@jackson2015; @dornelas2018], with significant impacts of species extinctions on functional diversity in some taxa and not others [@pimiento2017]. Some losses of functional diversity are indisputable on both paleocological and contemporary timescales such as continued trophic downgrading due to loss of large-bodied mammals, but implications of those losses for local diversity patterns are less clear [@estes2011; @smith2018] .

Assessments of broad-scale temporal change in functional diversity have previously been limited by a lack of functional trait data. The majority of work has therefore focused largely on system-specific studies with traits collected in situ. Ongoing efforts to assemble functional traits for a variety of taxa have made synthesis of existing community assemblage data and functional traits possible for the first time, providing initial insights into the ways functional diversity changes on a broad scale for specific taxa [e.g. fish, @trindade-santos2020; birds, @jarzyna2016; @barnagaud2017]. However, to date there has been no cross-taxa assessment of temporal functional change for a broad geographic and taxonomic extent.

Here we perform the first multi-taxa, multi-system assessment of functional diversity change through time. We focus on mammal, bird, and amphibian species as a significant subset of the world's biodiversity heavily impacted by anthropogenic change. While examining trends in plants, invertebrates, and other vertebrate species is of equal interest, trait data for those taxa raise additional challenges such as limited and biased species coverage [@fitzjohn2014], a lack of accepted species-level means, and differences in the types of traits collected. In order to ensure comparability across taxa in trait type and data quality we therefore focus on mammals, birds, and amphibians. Traits were intentionally selected to be representative of a species' Eltonian niche, thereby summarizing the functional role they play in the community [@wilman2014].

We assess thousands of mammal, bird, and amphibian functional diversity time series to determine whether or not there is a general trend of functional change, both in observed metrics and in metrics corrected for changes in species richness. We distinguish between three possible scenarios of functional change: 1) significant loss of functional diversity with accompanying shifts in other functional metrics, 2) no functional diversity loss, but significant shifts in other functional metrics, 3) maintenance of functional diversity and composition. Based on expectation due to human impacts, we expect to find a significant functional loss with further restructuring indicated by the additional metrics.

Results

We found no significant overall trend in species richness or summary functional diversity metrics (observed or standardized) (Fig \ref{fig:timeseriesPlot}). We did find a significant overall decrease in Jaccard similarity, indicating increasing turnover through time. Non-significant overall trends indicate that although some studies experience increasing or decreasing trends, the average trend across studies was plausibly 0 (Table \ref{tab:resultsTab}). Trends for different taxa, biomes, or realms were also non-significant with the exception of a significantly increasing trend for functional evenness of global studies (characterized by having samples on multiple continents), and a significantly decreasing standardized functional richness slope for freshwater studies. However, with only two global studies and two freshwater studies these results are not general trends. We found similar results for the CWM models, with a significant trend for only one trait, percentage of aerial foraging utilization, which showed a significant positive trend.

knitr::include_graphics(here::here("figures/3met_long.jpeg"))

read_csv(here::here("figures/man_model_table.csv")) %>%
  select(metric, everything()) %>%
  filter(group != "Residual" | is.na(group)) %>%
  mutate(metric = ifelse(metric == "Jaccard_base", "Jaccard", metric),
         metric = ifelse(metric == "S", "Species Richness", metric),
         effect = ifelse(effect == "ran_pars", "random", effect),
         group = case_when(
           group == "rarefyID:study_id" ~ "time series within study",
           group == "study_id" ~ "study",
           group == "Residual" ~ "residual",
           TRUE ~ group),
         term = case_when(
           term == "(Intercept)" ~ "Intercept",
           term == "year_scaled" ~ "Year",
           term == "sd__(Intercept)" ~ "SD Intercept",
           term == "cor__(Intercept).year_scaled" ~ "Corr(Intercept, Year)",
           term == "sd__year_scaled" ~ "SD Year",
           term == "sd__Observation" ~ "SD Observation",
           TRUE ~ term
         )) %>%
  mutate(across(c(estimate, std.error, p.value), round, 2),
         p.value = ifelse(p.value == 0, "<0.001", p.value)) %>%
  arrange(metric != "Species Richness", metric, !is.na(group), group, term == "Corr(Intercept, Year)", term) %>%
  dplyr::rename(grouping = group) %>%
  kable(format = "latex", caption = "Model estimates and statistics for general trend models for species richness, Jaccard similarity, and standardized functional diversity metrics. Additional model estimates, including CWM models, can be found in the supplement.") %>%
  column_spec(3, width = "2cm") %>%
  collapse_rows(columns = 1:3, valign = "middle") %>%
  kable_paper() %>%
  kable_styling(latex_options= c("scale_down", "hold_position"))

At the study level, 4 studies experienced a significant trend in species richness and only 10 of 54 studies for observed metrics and 3 of 54 studies for standardized metrics experienced a significant trend for a metric other than Jaccard similarity (Table \ref{tab:trendTab}). Most significant trends for observed functional metrics are in functional richness and disappeared after standardization, indicating that richness increases were likely due to changes in the number of species. However, the majority of the functional richness trends were negative, indicating that when change is occurring it is more often a loss. Hypothesis testing for study-level trends is likely affected by multiple testing issues and some trends identified as significant are therefore potentially spurious. Rather than interpreting changes in specific studies, we present these results as a general picture of how few studies experienced a significant trend and highlight that even that count is likely an overestimate.

library(tidyverse)

data.frame(sign = c("+", "-"), S = c("1", "3"), `Jaccard Similarity` = c("0", "37"), FRic = c("2", "6"), FEve = c("1", "0"), FDiv = c("0", "1"), `SES FRic` = c("0", "1"), `SES FEve` = c("0", "2"), `SES FDiv` = c("0", "1")) %>%
  column_to_rownames("sign") %>%
  kable(caption = "Number of studies that experienced a significant trend in each calculated metric out of 53 total studies.")

Study-level slopes for multiple metrics were significantly related to the duration and start year of studies. Slopes for species richness were significantly more negative with later start date and more positive shorter duration studies. Jaccard similarity and functional evenness both had significantly more negative slopes with more recent start year, whereas functional divergence was significantly more positive. Slopes for functional evenness were also significantly more positive for longer duration studies. Results were consistent between standardized and observed metrics with exception of functional evenness, which was negatively related to duration for observed data and positively related for standardized data. See supplement for estimates and p-values for all models.

Discussion

Our study represents the largest broad-scale multi-taxa assessment of functional change through time to date, giving a first look at aggregate and local trends in functional diversity in mammal, bird, and amphibian communities. Surprisingly, we did not detect an overall trend in any of the chosen functional diversity metrics. As with previous species-based syntheses, we also found no overall trend in species richness accompanied by increasing turnover through time [@dornelas2018], indicating that non-significant trends in functional metrics may be consistent with similar well-documented species derived trends. We found no evidence of systematic functional richness loss or functional change. A lack of trend for almost all realms, biomes, and taxonomic groups gives further evidence that directional functional change is absent from all systems observed in our dataset. Additionally, results from CWM models show that there were very few general shifts in functional trait values.

This striking result could be a product of two possible processes, one ecological and one methodological. Null trends appear to give strong evidence of systematic maintenance of functional structure due to common ecological processes, however multiple limitations of current approaches in synthesis could potentially be obscuring a true underlying global trend. We discuss both options further here.

Evidence of Ecological Processes

Communities demonstrated almost universal maintenance of functional composition. While the majority of the studies (\~70%) included in our data experienced significant species turnover, only three (for standardized metrics) experienced a significant shift in any functional dimension. This suggests certain characteristics of the functional space are maintained even in the face of significant change in species identity, specifically the size of the functional space occupied by the community (FRic), the distribution of species and individuals within that space (FEve and FDiv), and the mean of individual trait distributions (CWM). On average, species additions have similar functional characteristics as lost species and therefore maintain the structure of the functional space.

These results challenge assumptions that functional loss is the default state of all or even many communities. And while there is some evidence that functional loss may be more common than functional gains for communities experiencing functional change, there was no significant trend for even the longest running and most heavily impacted studies. The North American Breeding Survey for example is considered an authoritative dataset on the state of bird populations on the continent and underpins policy decisions about bird conservation [@pardieck2020; @rosenberg2019; @sauer2017]. No more robust dataset exists to capture North American avian community change, yet we detected no general shifts in functional structure across the dataset. Further, none of the 5 included studies that experienced a manual manipulation (e.g. burning, grazing exclosure, etc) experienced any significant functional trends.

Results are also seemingly inconsistent with predictions for trait shifts under global change. For example, mean body size is predicted to decrease as a result of climate change impacts and megafaunal loss [@sheridan2011], a phenomena which has already been well documented empirically and experimentally in multiple taxa [@caruso2014; @forster2012; @huss2019; @tseng2018]. While the species-level trait means used here are not appropriate for assessing intraspecific body size shifts, we would expect to see shifts in community-wide means due to local losses of large-bodied species. Instead, we found no evidence of a trend in CWM body size.

While we did not directly measure changes in rare species, our results further contradict likely scenarios of loss predicted due to rare species extinction. Rare species, defined by small populations and geographic restriction, are simultaneously more likely to be functionally distinct and at higher risk for extinction [@davies2004; @harnik2012; @loiseau2020]. Locally, communities losing functionally rare species should exhibit strong functional shifts as lost species can eventually no longer be replaced by functionally similar species [@leitão2016]. Observed patterns were instead consistent with species replacement by functionally redundant species from the species pool. Still, for many time series we likely did not have a large enough time window to capture community and species pool impoverishment due to extinction.

What does local maintenance of functional structure mean for ecosystem function? The vast majority of experimental and observational work links declines in function to declines in functional or species diversity [@brose2016; @cadotte2011; @duffy2007]. By those criteria very few communities in our dataset are in a state of concern for loss of functionality. However, shifts in metrics are only relevant if the underlying traits are those most critical for ecosystem function. We were limited in this analysis to the traits available rather than those with strong empirical links to function. Similarly, the dimensions of functional space most important for ecosystem function are still a topic of ongoing debate, and at least some known aspects important for multifunctionality were not measured here [e.g. dispersion, rarity, abundance of dominant species, @bagousse-pinguet2021]. Still, the fact that we observed so many communities maintaining structure across the most commonly used metrics for linking biodiversity and function calls into question how previous work translates to natural communities. Metrics need to be both closely linked to changes in ecosystem function and also experiencing shifts in natural communities to be meaningful.

Potential Methodological Limitations

Here we approach the question of functional change using the best available data and biodiversity synthesis approaches. However, a number of gaps in best practices may be obscuring a true underlying trend. First, the BioTIME database, while the most comprehensive data source of time series available, is limited in temporal and geographic scope. Most time series span only a few years (Fig \ref{fig:taxaMap}) and may not provide the statistical power necessary to detect trends. The database is also not a representative sample of the world's biodiversity or areas of greatest threat [@gonzalez2016; @vellend2017], and the subset of data in this study exhibits a strong Northern Hemisphere bias. We may simply not have data from those areas experiencing the greatest perturbation [@hughes2021], particularly scenarios of conversion to urban, human-dominated landscapes. While evidence from other work shows even disturbed communities can maintain functional structure [@edwards2013; @matuoka2020], these results should not be interpreted as evidence of low functional impact in areas of heavy human disturbance.

Second, despite using the most comprehensive trait databases for these taxa, we were still limited to species-level means of the traits deemed important by database creators. The importance of intraspecific variation is well documented [@desroches2018; @violle2012], however individual-level traits are rarely collected alongside monitoring data, especially for the longest running efforts. Species-level traits may be obscuring more subtle shifts in the trait space happening within species. Likewise, available trait data may not capture the traits experiencing the greatest change.

Third, while we use here the most common metrics for describing functional diversity, they do not measure some potentially important aspects of the functional space. Most notably, the summary metrics we calculated do not capture shifts in the location of the functional space as a whole. For example, two communities could have very similar metric values but no overlap in their trait spaces. This is especially relevant in the context of biodiversity change as a species loss could be replaced by a species with very different functional attributes, but the replacement would go undetected if the new species expanded the trait space by the same degree and had similar abundance. This scenario may be common in communities tracking changing environmental conditions. While trait CWM's capture axis shifts, approaches for assessing multidimensional shifts in functional space are still relatively new [@barros2016; @blonder2018; @mammola2019] but could shed critical insight into functional composition changes of this nature.

Policy Implications

Our results should not be interpreted as an indication that the ongoing biodiversity crisis is less severe than previously described, or that there is no concern for functional change as a result of anthropogenic impact. These findings do not negate a substantial body of work linking functional degradation to direct human intervention in the form of land use change and intensification or habitat fragmentation [@flynn2009; @magioli2021; @tinoco2018], but rather illustrate trends for communities experiencing background levels of environmental change. Rather than assuming functional structure will be maintained in areas of concern, our work indicates that when measurements of functional diversity show significant shifts, it should be considered an indication of substantial community change and outside the normal expectation.

Future Work

Here we make a significant first step in establishing a general trend for functional diversity through time across a variety of taxa and systems. We present the conclusion best supported by available data and acknowledge that it is still too early to confidently distinguish between true ecological pattern and methodological limitations. The most pressing next step is for intentional and targeted data collection efforts. We join others in the call for increased monitoring in under sampled areas and continued efforts to centralize existing data sources [@gonzalez2016; @hughes2021; @vellend2017]. Data that fill geographic, taxonomic and trait gaps should be prioritized over further collection of data that replicate existing biases. One relatively low-cost high-reward data investment is collation of additional species-level trait means. Intentional trait selection is critical for linking functional patterns to ecological processes [@zhu2017], however synthesis is constrained to the traits in a few taxa-specific databases. Trait collection should explicitly consider existing frameworks for linking traits to processes [e.g. the response and effect framework @lavorel2002] to facilitate clear ecological interpretation of potential functional changes.

Material and Methods

Data

We obtained mammal, bird, and amphibian time series from the BioTIME database, a global repository of high quality assemblage time series. All studies included in the database follow consistent sampling protocols and represent full assemblages rather than populations of single species [@dornelas2018]. Following best practices for the database [@blowes2019], studies with multiple sample locations were split into individual time series following a standardized spatial scale. Scale was set by a global grid with cell size determined based on the sample extent of studies with only a single location [see @dornelas2018 for details on how sample extents were defined], with the area of each cell set to one standard deviation away from the mean of the single extent locations. All samples from a study within a single cell were considered to be a single time series, and species abundances were combined for all samples.

tibble(Taxa = c("Mammals", "Birds", "Amphibians"),
       Timeseries_Count = c(48, 2380, 11),
       Species_Count = c(184, 700, 184),
       Trait_Source = c("Elton Traits", "Elton Traits", "Amphibio"),
       Traits = c("body mass, diet, active diel period", "body mass, diet, nocturnality, \nforest foraging strata, pelagic specialist", "habitat, diet, active diel period, \nactivity seasonality, body mass, body length, \nmin maturation size, max maturation size, \nmin offspring size, max offspring size, \nreproductive output, breeding strategy")
) %>%
  mutate_all(linebreak, align = "l") %>% 
  kableExtra::kable(format = "latex", col.names = c("Taxa", "Number of \nTime Series", "Number of \nSpecies", "Trait Source", "Traits"), escape = FALSE, 
               caption = "Summary of the data in the final trait database.") %>%
  kable_classic() %>% 
  kable_styling(latex_options= c("scale_down", "hold_position")) #%>%
   #kableExtra::column_spec(2, width = "5cm") #%>%
  # kableExtra::column_spec(4, width = "10cm")

We gathered trait data from the Elton Trait Database [mammals and birds, @wilman2014] and Amphibio [amphibians, @oliveira2017]. These databases include species-level means for traits that partially represent species' multifaceted function in the community including body size, diet, and behavioral characteristics. For the full list of traits included in the analysis for each taxon see Table \ref{tab:traitSources}. Multiple traits (i.e. diet, foraging strata, activity seasonality, active diel period) were broken down into percentage or binary use for each level.

In order to ensure taxonomic consistency across datasets, BioTIME species were paired with trait data based on their species identifier from the Integrated Taxonomic Information System database (retrieved 09-15-2020 from the on-line database, https://doi.org/10.5066/F7KH0KBK), obtained through the taxadb R package [@norman2020; @rcoreteam2021]. If more than one species in the assemblage data resolved to the same identifier, observations were considered the same species. For trait data, traits for all species of the same identifier were averaged. Only studies with at least 75% trait coverage were included and observations for species with no trait data were excluded. In order to have a sufficient number of species to calculate functional diversity metrics, years with fewer than 5 species observed were also excluded.

Many studies had a variable number of samples within years. To account for this inconsistency in sampling effort we used sample-based rarefaction by bootstrap resampling within years for each time series based on the smallest number of samples in a year for that time series.

Our final dataset included 2,443 time series from 53 studies in 21 countries and 15 biomes and 13 different traits (Fig \ref{fig:taxaMap}). The earliest sample was in 1923 and the most recent was in 2016. Only four studies (consisting of 11 time series) came from Amphibian studies due to the limited availability of amphibian time series and low species richness values for assemblages (Table \ref{tab:traitSources}). Amphibians are of particular concern due to impacts of habitat loss and pollution [@gibbons2000], so we include data while acknowledging that general inference for amphibians as a clade is not possible with the time series available. For a full breakdown of studies and their characteristics, see the supplement.

knitr::include_graphics(here::here("figures/study_map_hist.jpeg"))

Diversity Metrics

We calculated yearly metrics of functional and species diversity for each time series. Species-based metrics include species richness (S) and Jaccard similarity (J) as a measure of turnover. Jaccard similarity was calculated relative to the first observed year for a time series. A negative trend in J would therefore indicate increasing turnover.

Functional diversity metrics were calculated using the dbFD function from the FD R package [@laliberté2010]. Here we report functional richness (FRic), functional evenness (FEve), and functional divergence (FDiv) which together describe three complementary characteristics of the functional space [@mason2005; @hillebrand2009]. FRic assesses the volume of the trait space occupied by species in the community, with higher values indicating communities with species of more extreme trait values. FEve describes how species are distributed across the trait space and how abundance is distributed across species. Higher values of FEve indicate more even spacing of species in the trait space and individuals across species. FDiv measures the degree to which species and their abundances maximize differences in the functional space. Higher values of FDiv therefore correspond to communities where many highly abundant species are on the edges of the trait space. We also calculated the community-weighted mean (CWM) of included traits to examine shifts in the distribution of each trait.

All available trait data for each study were included in functional diversity calculations with the exception of traits that were the same value for all observed species in the study. For variables with multiple levels each level was included as a separate trait axis. Continuous traits were z-score scaled to give each trait equal weight in the trait space [@leps2006; @schleuter2010]. The number of trait axes was limited to the maximum number of traits that fulfills the criteria $s >= 2^t$, where $s$ is the number of species and t is the number of traits. This restriction allows for a sufficient number of axes to capture the trait space while maintaining computational feasibility [@blonder2018]. Metrics incorporated weighting based on species abundance where available (three studies were presence only).

Null Models

To assess functional change independent of species richness we calculated the standardized effect size (SES) for each of the three summary functional diversity metrics (FRic, FEve, FDiv) from null estimates [@swenson2012]. Null model corrections allow us to assess the degree to which the observed functional diversity metric deviates from the value expected by chance in a randomly assembled community. Null estimates were calculated for each rarefied sample by randomly sampling species from the species pool for each year and randomly assigning observed abundances to species. Species pools included all species observed for a time series. This process was repeated 500 times to get an estimate and standard deviation of the null expectation for the metric for each rarefaction sample for that time series. We used these values to calculate SES using the following formula: $SES = [F_{obs} - mean_{(F_{null})}]/SD_{(F_{null})}$. We then calculated the median SES estimate for each metric from all the rarefaction samples for a time series. SES estimates can be interpreted as how much of the functional characteristic (richness, evenness, divergence) was observed beyond what was expected by chance for a community of that species richness.

Analysis

We estimated general trends for each diversity metric using a linear mixed effects model with a random slope and intercept for each study and each time series nested within the study. All time series with data for a given trait were included in the model. Trends could not be estimated for the CWM's of three traits due to limited data: maximum age at maturity, minimum age of maturity, and minimum litter size. We fit 24 individual CWM models.

For all other species and functional diversity metrics (S, J, FRic, FEve, FDiv) we obtained study-level estimates of temporal change from the Best Linear Unbiased Predictors (BLUPs) for each overall trend model. BLUP's provide estimates for the conditional mean and variance of a random effect from which we calculated 95% confidence intervals to determine significance of study-level slopes.

To test for trends within and between different levels of taxa, biome, and realm we fit separate models with each of those factors added as a predictor to the original model structure. We estimated within-level slopes and calculated between-level contrasts using the emmeans package [@lenth2021]. We assessed the impact of time series duration and start year on study-level trends using general linear models with duration and start year as predictors. All models were executed using the lme4 package in R and p-values were calculated by Satterthwaite's degrees of freedom method using the lmerTest package with a significance level of $\alpha = 0.05$ [@bates2015; @kuznetsova2017; @rcoreteam2021].

References

karinorman/biodivTS documentation built on Dec. 23, 2024, 10 p.m.