Spatial Scaling Patterns of Functional Diversity

Ecology, biogeography and conservation biology, among other disciplines, often rely on species identity, distribution and abundance to perceive and explain patterns in space and time. Yet, species are not independent units in the way they interact with their environment. Species often perform similar roles in networks and their ecosystems, and at least partial redundancy or difference of roles might explain co-existence, competitive exclusion or other patterns reflected at the community level. Therefore, considering species traits, that is, the organisms’ functional properties that interact with the environment, might be of utmost importance in the study of species relative abundances. Several descriptive measures of diversity, such as the species-area relationship (SAR) and the species abundance distribution (SAD), have been used extensively to characterize the communities and as a possible window to gain insight into underlying processes shaping and maintaining biodiversity. However, if the role of species in a community is better assessed by their functional attributes, then one should also study the SAR and the SAD by using trait-based approaches, and not only taxonomic species. Here we merged species according to their similarity in a number of traits, creating functional units, and used these new units to study the equivalent patterns of the SAR and of the SAD (functional units abundance distributions - FUADs), with emphasis on their spatial scaling characteristics. This idea was tested using data on arthropods collected in Terceira island, in the Azorean archipelago. Our results showed that diversity scales differently depending on whether we use species or functional units. If what determines species communities’ dynamics is their functional diversity, then our results suggest that we may need to revaluate the commonly assumed patterns of species diversity and, concomitantly, the role of the underlying processes.


INTRODUCTION
Species diversity, or simply biodiversity, encompasses several scales, from the genetic (phylogenetic diversity at the species level) to that of species (taxonomic diversity), populations and ecosystem functions (functional diversity) (Wilcox, 1984;United Nations Environment Programme [UNEP], 1992). One of the major goals of ecology is to describe this diversity at different spatial and temporal scales and seek for the processes shaping and maintaining it. However, species do not exist in isolation, they interact with each other and are influenced, and influence, abiotic processes in the environments where they exist. Therefore, a major challenge is to connect the description of taxonomic diversity with the processes by which species interact with their environments. Functional diversity has the potential to link these two views of a community (Asner et al., 2017). Our main purpose here is to study patterns of functional diversity using tools commonly used to assess patterns of species taxonomic diversity with an emphasis on its scaling properties.
By functional diversity we mean: "the values and range in the values, for the species present in an ecosystem, of those organismal traits that influence one or more aspects of the functioning of an ecosystem" (Tilman, 2001). As (Petchey and Gaston, 2006) pointed out, this is a relatively narrower definition than that found in original studies, and that emphasizes the importance of "measuring functional trait diversity, where functional traits are components of an organism's phenotype that influence ecosystem level processes" (Petchey and Gaston, 2006). An important aspect highlighted by functional diversity is that species often overlap in their traits, leading communities with different species (taxonomic) composition to having a similar functional composition . Therefore, because the emphasis is no longer on taxonomic differences, as recently highlighted by Wong et al. (2019).
The main novelty of our work is its emphasis on the spatial scaling attributes of functional diversity by looking at the equivalent to the species-area relationship (hereafter SAR) that we will call the functional units-area relationship (hereafter FUAR), and the equivalent to the species abundance distribution (hereafter SAD), that we will call the functional units abundance distribution (hereafter FUAD) and where we adopted the following definition of functional units: a set of organisms sharing the same unique combination of attributes. Concerning the FUAD, as we elaborate later, our purpose is not to discuss which distribution gives the best fit at one spatial scale, but how the distributions change as a function of scale. To do this we will use the (raw) moments of the distributions. This is done with an aim to providing new tools and, importantly, to discerning new quantitative patterns associated with functional diversity.
One the most studied patterns of species taxonomic richness is the species-area relationship (e.g., Rosenzweig, 1995). The species-area relationship describes how the number of species changes as a function of area size. However, as shown by Scheiner (2003), there are different types of SARs depending on how data are collected and presented. For example, some SARs describe how the number of species change as a function of the size of islands in an archipelago (Whittaker et al., 2014), while others describe how the number of species accumulates as sampled (nested or non-nested) areas increase (Matthews et al., 2016); in this study we look at the SAR of non-nested areas. An important aspect of SAR studies is that they implicitly assume a scaling approach and try to identify a pattern. Which function best describes this scaling pattern has been a topic of contention among ecologists, as reviewed by Rosenzweig (1995). One of the most popular functions is a power law of the form, S = cA z , where S is the number of species, A the area size, and c and z constants to be determined. In our analyses we will compute the equivalent to the SAR but using species functional units instead of species and will compare functional units-area relationships, the FUAR, and taxonomic-based SAR curves of the same dataset. Smith et al. (2013) have already introduced the concept of functional-diversity-area relationship, to explore the scaling properties of individual traits alone, and when extending the concept to multiple traits they used the concept of convex hull model introduced by Cornwell et al. (2006), while here we group taxonomic species based on their traits to form functional units.
Although species richness and how it changes as a function of sample or area size are important attributes of a community, they do not provide information on the species relative abundances. There are several ways to account for species abundance, but probably the most intuitive, and the one we will use here, is the number of individuals. Typically, species relative abundances are depicted using a histogram, called a species abundance distribution, relating the number of species (the y-axis) with the number of individuals (the x-axis), with the latter usually expressed on a logarithmic scale of base 2 (e.g., McGill et al., 2007); logarithms are used in order to accommodate the wide range of abundances usually encountered in a community sample. In other words, a SAD answers the question of how many species exist within a given range of number of individuals. In line with our previous delineated approach of functional units, we will analyze SADs based on species functional units, the functional units abundance distributions, or FUAD, and, again, we will compare these distributions with the taxonomic SADs of the same dataset.
In contrast with the SAR, that looks at the number of species at different scales, the SAD is often computed and studied at one single scale (e.g., McGill et al., 2007). However, Borda-de-Água et al. (2012Borda-de-Água et al. ( , 2017 argued that SADs should also be studied as a function of scale, because histograms (i.e., the distributions) change considerably depending on the size of a sample (or area size). The visual description by Preston (1948) of a veilline progressively revealing more classes of the distribution as sample size increases, provides a good visual representation of the way the distributions change: for small sample sizes SADs are usually monotonically decreasing functions with the maximum occurring for the singleton class, while for larger sample sizes the distributions become bell-shaped. Once we realize that the SADs change as a function of sample size, an important task is to characterize that change, that is, the scaling properties of the distributions. Borda- de-Água et al. (2012, 2017) used the raw moments of the distribution to describe the scaling of the SADs (see section "Materials and Methods" for the definition of raw moments). The characterization of the distributions at different scales based on their moments, enables the SADs to be extrapolated to larger areas, whose sizes pose sampling problems because of economical or other practical reasons. Here, and for the same reasons, we will study how SAD based on functional groups instead of species (i.e., functional units abundance distributions) change as a function of scale and will use the raw moments to describe their scaling properties.
In summary, the rationale of our study is the following. If functional units provide a more accurate description of how species interact among them and the environment, then the above scaling patterns of species richness and relative abundance should be reformulated in terms of the species' functional units and not of the species taxonomical identities. As we will see, using data on arthropods collected in the Azores archipelago, describing the scaling proprieties of the diversity of a community using species' functional units leads to quantitatively different results from those obtained when using taxonomic species.

Study Sites and Sampling Strategy
We used data on arthropod species collected in Terceira island, Azores (Figure 1). The Azores archipelago (37 • -40 • N; 25 • -31 • W) consists of nine islands and it is one of the world's most isolated archipelagos. Terceira is the island for which we have the most data (see Borda-de-Água et al., 2017), hence we will focus our analysis on this island. Terceira is the third largest island of the archipelago with an area of 402 km 2 and it is currently estimated to be 0.40 million years old (Calvert et al., 2006). We gathered data in native forest reserves from 39 transects, each 150 m long and 5 m wide. We sampled all transects with the same sampling effort using standardized methods; for more details see Borges et al. (2005) or Rigal et al. (2018). We have identified 181 arthropod species for Terceira island, and for each species we have also recorded the number of individuals sampled.

Species Identification
The arthropods collected were identified in the laboratory at the species level for the taxa Araneae, Opiliones, Pseudoscorpiones, Diplopoda, Chilopoda and Insects (excluding Collembola, Diplura, Diptera and Hymenoptera). Taxonomic identification was performed in two steps: (i) trained parataxonomists sorted samples into morphospecies (i.e., recognizable taxonomic units, sensu Oliver and Beattie, 1996) using a non-complete reference collection; (ii) experienced taxonomists assisted in the identification of the morphospecies. All species were classified as indigenous or exotic. Indigenous species comprise Azorean endemics and other native non-endemics. Exotic species are those considered to have colonized via human mediation, many of which have a cosmopolitan distribution (Borges et al., 2010). As in some of our previous studies (e.g., Borges et al., 2005;Gaspar et al., 2008), we dealt with unidentified morphospecies as follows. When other species in the same genus, subfamily or family were present in the archipelago and all belonged to the same colonization category (according to Borges et al., 2010), the unknown morphospecies were classified similarly. If no information was available, we assumed the species to be native since exotics are usually widespread and easier to identify (Borges et al., 2010).

Functional Data
Different species interact with the environment and other organisms in different ways. Yet, standardization across taxa is not only possible but provides important insights (e.g., Chichorro et al., 2019, Chichorro et al., 2020. Often a trait-based approach has to be adapted to accommodate such differences, by using different proxies that are equivalent for contrasting taxa. In our case, as we only studied arthropods, which share many characteristics, standardization was relatively easy. In fact, in this last decade recent efforts have been made to develop trait database for arthropods as a whole (See Schweiger et al., 2005;Simons et al., 2016;Rigal et al., 2018). Functional traits were then selected based on the known life-history of the taxa analyzed as representing different ways the organisms can be affected in their probability of survival and hence abundance in the studied forests. This way we are implicitly studying mechanisms that affect species abundances in the native forests of the Azores and how such abundances change depending on the spatial scale.
For all arthropod species, we collated data for five traits namely, body size, dispersal abilities, trophic level, microhabitat and origin. Body size was measured for the individuals sampled in this study and was coded using ordered size classes delimited by the intervals (0, 0.5), (0.5, 1.5) . . . (40.5, 41.5). Both dispersal abilities and trophic level were collected from an extensive literature search, including manuscripts with the first descriptions of the species, first species records for the Azores, brief notes and ecological studies. Information was also obtained from experts who have identified the specimens or from experts of a given taxonomic group when information for a particular species was not available. Functional information was assigned to each species according to their adult characteristics, except for Lepidoptera, where traits were assigned with reference to the larval stage. For the unidentified morphospecies, we assigned functional traits of the nearest taxonomic resolution (genus, family), except for body size which was measured directly from the individuals.
Specifically, for the trait "dispersal ability, " the species were categorized into high and low dispersal classes based on ecological attributes and morphological characteristics. This could be for instance, the presence of active wings for beetles (Coleoptera) and bugs (Hemiptera), ballooning for spiders and evidence of flying ability for endemics and general natural history guides for the other species. To be considered as a good disperser, a species has to be able to disperse between fragments of native forest and surpass the current matrix of man-made habitats (e.g., pastures). For origin we identified three classes: "endemic, " "introduced, " or "native non-endemic." One can argue that the "origin" of a species is not a relevant trait to how a species functions within a community. We use this category because introduced species may have (new) traits that are not present among the native species and, for a similar reason, we identify endemic species because they may have unique traits adapted to particular conditions of the archipelago or, in the case of single island endemics, of an island (see also Rigal et al., 2018).
Thus, each species was categorized based on five traits: • Average body size classes; • Dispersal ability, with two classes: "high" or "low"; • Trophic level, with two classes: "herbivore" or "predator"; • Microhabitat, with two classes: "canopy" or "pitfall"; • Origin, with three classes: "endemic, " "introduced, " or "native non-endemic." For each species we determined the combination arising from the five above traits. For instance, a species could have the following combination (average body size = 3.5, dispersal = "low, " location = "canopy, " origin = "native, " trophic level = "herbivore"). From a functional perspective, we think of all species with the same traits' combination as being just one functional unit. Thus, we identified all species with the same traits' combination and added their number of individuals to form a functional unit; obviously, the number of functional units is necessarily smaller or equal to the number of species.

Statistical Procedures
In order to obtain the species and functional units-area relationships (FUAR) we started from a random transect and added the species of the nearest transect, and then of the next nearest transect until reaching the desired number of transects. By adding the nearest transect we better replicate a spatial sampling process that avoids overinflating the effect of beta diversity that would happen if sites were added, irrespectively, of spatial distance. The final SAR and FUAR curves are the average of the curves obtained from different starting transects. This procedure leads to one of the six types of speciesarea relationship described by Scheiner (2003), and it has the advantage of retaining information of the relative distance between the transects.
We construct histograms of functional units abundance distributions (FUAD) as the classical taxonomic species abundance distributions (SADs). Figure 2 shows some examples. On the y-axis we plot the number of functional units and on the x-axis the number of individuals on a logarithmic scale of base 2, with the classes of the number of individuals defined as 1, 2-3, 4-7, et seq. When a species abundance distribution resulted from the accumulation of data from two or more plots, we always added the data from the nearest transect, as we did for the SARs. Also, as before, the final SAD and FUAD curves are the average of the distributions obtained from different starting transects.
In order to characterize the scaling properties of the FUADs, we use their raw moments (occasionally, when there is not risk of confusion, we call them simply "moments") as we have done previously for SADs (Borda-de-Água et al., 2012(Borda-de-Água et al., , 2017. If a community is made of S FG functional units and x j = log 2 (X j ), where X j is the number of individuals of the j functional unit, then the raw moment of order n, M n , is calculated from Frontiers in Ecology and Evolution | www.frontiersin.org notice that M 0 is equal to 1, and M 1 corresponds to the average. Other, more familiar moments, such as the variance (the second central moment), can be obtained as combinations of the raw moments.

The Extrapolation Procedure
Here we used an improved method to extrapolate the distributions based on that described in Borda-de-Água et al. (2012,2017). The basic idea consists of estimating the moments, using formula (1), up to a given order n, plotting these values and fitting the curves of each order moment and then, using these curves, estimating the moment values for the new areas. In order to reproduce the abundance distribution based on the extrapolated moments, we used an improved version (Mukundan, 2004)  We start by presenting the SAR and the FUAR (Figure 3); because all transects have the same dimensions, we equate the number of transects with area size. The obvious difference between the two curves is that the FUAR in a log-log plot starts exhibits a plateau, showing clear signs of starting to "saturate, " while the traditional species-area curve reveals the typical (almost) straight line in a log-log plot, that is the number of taxonomic species keeps increasing until the very last transect considered. In a sense, given that the number of functional units is a subset of the number of species, one could expect this to happen, but notice that the possible plateau only starts to appear at a large number of transects.
We showed in Figure 2 the SADs and FUADs: plots a and c were obtained after averaging the distributions of individual transects and plots b and d when we added the data from all the transects. Notice that at the level of one transect (plots A and B) the SAD and FUAD are very similar, being almost monotonically decreasing curves. From here on, we will describe curves that have such a monotonically decreasing shape as "logseries-like, " but we do not ascertain that the logseries distribution would provide the best fit, it is, hence, just a terminology to describe the general shape of the curve. On the other hand, when all transects are added together they are very different (plots C and D). Specifically, while the SAD of the 39 transects retains a logseries-like shape, the FUAD has its maximum at an intermediate class and, although it still has a large number of singletons, overall it has a bell shape. In other words, the FUAD evolves faster from a logseries-like distribution to a bell-shaped distribution than the species-based SAD. From here on, we will describe curves with a bell-shape distribution, thus with the maximum at intermediate classes, as "lognormal-like" without implying that the lognormal would provide the best fit; as with the terminology for the "logserieslike" distributions, it just a general description of the shape of the distribution.
We estimated the moments of the FUAD using Eq. 1, Figure 4. The overall behavior of the moments of the functional units as a function of the number of transects is not very different from those of taxonomical diversity (Borda-de-Água et al., 2017), that is, above a certain number of transects (or area) the moments are approximately linear in log-log scales. Using the procedure described in the Supplementary Material, we extrapolate these moments for a larger number of transects and used them to reconstruct the FUAD (Figure 5); for sake of comparison we applied the same procedure to the traditional taxonomical SAD. The SAD evolves from an almost monotonically decreasing function to a curve with a clear peak when we extrapolate the distribution for 70 transects (the green line in Figure 5A). However, the appearance of a peak for the SAD is slow compared to that of the FUAD, whose peak for intermediate classes is already clear in the FUAD for 39 transects. When we increase the number of transects the maximum peak for the extrapolated curves shifts to higher abundance classes, and there is a reduction in the number of singletons.

DISCUSSION
Under the premise that functional units have a more meaningful correspondence with the underlying processes of maintenance and generation of species diversity (Tilman, 2001), we applied several procedures commonly used on studies of communities but using, instead of taxonomical species, functional units, that is, entities resulting from aggregating species with identical traits. Common criticisms often leveled at studies of functional diversity can also be raised here, namely, that the results depend on the choice of the traits, on the way traits are coded (e.g., number of categories identified) and on the scales used to identify traits (e.g., Ricotta, 2005). Undoubtedly, we would have obtained quantitatively different results if we had chosen a different set of traits or a different scale for the average body size, for instance. However, we describe this work as a first exploratory attempt to determine how diversity scaling changes depending on whether we use species or functional units.
These exploratory results showed that the scaling behavior of the species and of the functional units is quantitatively different. For instance, the FUAR relationship does not have a linear relation when plotted on log-log scales and shows a clear plateau, and contrary to the SAR, when area increases the FUAD develops a peak earlier than the corresponding SAD. Concerning the FUAR we acknowledge that its shape is a function of the number of traits chosen. For instance, if we had chosen fewer traits, the plateau would have been reached even earlier, because there would have been fewer functional units. Although it is outside the scope of the work, a pattern that is worth analysing in the future is how the FUAR shape changes in a community as a function of the number of traits used. Furthermore, to overcome the problem of trait dimension in driving FUAD patterns, we could also use a null model approach to test whether FUAR was simply a product of species richness (see Whittaker et al., 2014, for example).
To understand the faster transition in a FUAD, it is important to understand the transitions in a typical SAD. Using computer simulations, Borda-de-Água et al. (2007) conclude that the SADs undergo three regimes when area increases (their Fig. 17.9). For very small sample sizes, the distributions are monotonically decreasing curves, thus the maximum occurs for the singleton class, and they can be described as logseries-like distributions. When the sample size increases the number of singleton species, and those of other rare classes, decreases and the distributions develop a peak for intermediate abundance classes, having a lognormal like shape, although with some singletons and other rare species still being present. This shift in the shape of the species abundance distribution can be observed in real data, such as the data on tropical tree species of a 50 ha plot in Barro Colorado Island, Panama; see for example Borda-de-Água et al. (2012; Figure 1). The explanation for such a transition is the following: (i) when the number of individuals is small and several species are present, the number of individuals per species is small and most species are rare; (ii) once more individuals are collected some, if not most, belong to species already present, thus, species accumulate more individuals and they move away from the rarity classes toward those classes corresponding to more abundant species, therefore, overall, intermediate abundant classes start having more species and the curve becomes more bell shaped. When sample size increases further, the distribution becomes again monotonically decreasing, resembling again a logseries distribution. The explanation for such transition is the following: most species are rare for very small sample sizes, but only a few are really rare in the (meta)community, therefore, those that are not truly rare move to more abundant classes when sample sizes increase, but those that are truly rare remain in the classes of the rare species. For very large samples, we collect species from different regions of the community, each region containing rare species, thus at the level of the entire community there are a large number of rare species that are added to the rare abundance classes. The presence of a large number of rare species can also be the result of the temporal evolution of community. In fact, hyperdominance (i.e., few species having most of the individuals) may happen due to few winners dominating the community and the addition of mostly rare species to the regional pool with speciation operating on long timescales (McGill et al., 2019). As new species are added at a regional scale through evolution, hyperrarity and hyperdominance are generated. Consequently, in time, the number of such (truly) rare species can be very high and the SAD is a monotonically decreasing function. We suggest that the latter transition from a lognormal to a logseries is difficult to observe in real data because it requires large sample sizes. However, in addition to the evidence from simulations (Borda- de-Água et al., 2007), theoretical results based on Neutral theory (Hubbell, 2001) and data on Amazon tree genera do suggest that the SAD for very large communities is logserieslike (Ter Steege et al., 2006;Hubbell et al., 2008;Hubbell, 2013;Ter Steege et al., 2013).
The previous description for the three regimes observed for SADs may not apply when we look at the abundance distributions of functional units. The first transition from the logseries-like to the lognormal-like may still be present, as revealed by Figure 5, though it occurs faster than the corresponding SAD. This is not surprising given that functional units are collections of species, thus a faster increase in the number of intermediate abundance functional units and a decrease in the number of rare functional units is not surprising. However, we speculate that the transition from lognormal-like to logseries-like for large scales may never occur for FUADs, because, though we expect more rare species to enter the sample, these may have the same set of traits of other species already present in the community, thus the number of functional units in the rare classes will not increase. In fact, we expect the number of functional units in rare abundance classes to decrease. In other words, although the number of species is likely to increase when the number of sampled individuals increase, the number of trait combinations, and thus that of functional units (but not their abundance as measured by the number of individuals), is likely to remain approximately constant. If this is indeed the case, and if the functioning of the community is (partially) determined by the relative abundance of the functional units, then communities at large scales are mainly characterized by lognormal-like distributions and not logserieslike distributions. This has implications for the investigation of the mechanisms and patterns that govern communities. The important mechanisms shaping a community may lead to FUADs that are inherently lognormal-like.
We finalize with an observation that points to future work. Once we group species based on their similarity, it is implicit that some species may be redundant. Such a view, however, should be taken carefully because we only identified a subset of all possible traits. It is possible that some non-identified traits may correspond to a species role in the ecosystem that is irreplaceable, and whose absence could lead to profound transformations in the community, including the loss of other species. In this respect it is worth noting that several works have shown that communities exhibit approximately the same number of species, S, and individuals, N, over time, that is, there is a regulation of temporal trends of species diversity (Brown et al., 2001;Gotelli et al., 2017;Blowes et al., 2019). This is a remarkable result, although the mechanisms of such regulation are not (fully) understood. One question that remains is whether such regulatory trends also occur if temporal community diversity is analyzed in terms of functional diversity. In fact, some studies have shown that functional diversity changes over time in some communities (e.g., Mendez et al., 2012;White et al., 2018;Dolezal et al., 2020). However, if regulation is even stronger when we consider functional units, this suggests that one possible way to identify functional units, would be to consider several trait combinations, and identify which trait combinations lead to the smallest variation of the number of functional units and their abundances over time. This approach could provide an objective choice of the traits and respective functional entities. Finally, we define functional units using categorical traits but some traits are usually recorded as continuous variables (e.g., morphological traits), therefore requiring alternative measures of functional units such as clustering methods (Villéger et al., 2012).

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found here: https://doi.org/10.3897/BDJ.4.e10948.

AUTHOR CONTRIBUTIONS
LB-Á, PC, and PB designed the original research. SA did most of the analytic work. All the authors contributed to writing up the article.

ACKNOWLEDGMENTS
We thank H. M. Pereira for discussions related to functional diversity metrics.