The importance of grain and cut-off size in shaping tree beta diversity along an elevational gradient in the northwest of Colombia

Johanna Andrea Martínez-Villa,Sebastián González-Caro and álvaro Duque*

Abstract Background: Species turnover (β-diversity) along elevational gradients is one of the most important concepts in plant ecology. However, there is a lack of consensus about the main driving mechanisms of tree β-diversity at local scales in very diverse ecosystems (e.g., Andean mountains), as well as how the sampling effect can alter β-diversity estimations. Recently, it has been hypothesized that patterns of change in β-diversity at local scales along elevational gradients are driven by sampling effects stemming from differences in the size of the species pool rather than by underlying community assembly mechanisms.Thus, we aim to evaluate the relative extent to which sampling effects, such as species pool size, grain size, and tree size cut-off, determine species sorting, and thus, the variability of β-diversity at local scales along elevational gradients in the northwest of Colombia.Results: Using 15 1-ha permanent plots spread out along a 3000 m elevational gradient,we used standardized βdeviation to assess the extent to which either sampling effects or the community assembly mechanisms determine the changes in species composition at local scales. Standardized β-deviation was measured as the difference between the observed and null β-diversity divided by the standard deviation of the null β-diversity. We found that the magnitude of change in local β-deviation along the elevational gradient was significant and dependent on the employed spatial grain size and tree size cut-off. However, β-deviation increased with elevation in all sampling designs, which suggests that underlying community assembly mechanisms play a key role in shaping local βdiversity along the elevational gradient.Conclusions: Our findings suggest that grain size enlargement and the inclusion of trees with small diameters will improve our ability to quantify the extent to which the community assembly mechanisms shape patterns of βdiversity along elevational gradients. Overall, we emphasize the scale-dependent nature of the assessment of βdiversity. Likewise, we call for the need of a new generation of enlarged forest inventory plots along gradients of elevation in tropical forests that include small individuals to improve our understanding about the likely response of diversity and function to global change.

Keywords: Andean forests, Null models,Species pool, Species sorting, Sampling effect

Background

Spatial turnover in community composition (β-diversity)along elevational gradients has been one of the most striking and studied patterns in ecology (Whittaker 1960;Lomolino 2001; Rahbek 2005). In tropical mountain systems, β-diversity is expected to decrease with elevation(Tello et al. 2015) due to the influence of different community assembly mechanisms that could vary along the elevational gradient (Laiolo et al. 2018). Overall, different assembly mechanisms,such as dispersal limitation(Condit et al.2002),species sorting(Qian and Ricklefs 2007),habitat specialization(Janzen 1967;Jankowski et al.2009),and priority effects (Chase 2010; Fukami 2015), have been thought to explain the spatial turnover in the composition of plant communities. However, sampling effects associated with the size of the species pools and the regional abundance distributions have recently been proposed as the main cause of the observed decreased in β-diversity along elevational gradients (Kraft et al. 2011). In other words, the observed variation in β-diversity along steep elevational gradients may be primarily driven by differences in the size of the species pools and the number of individuals per specie generated by biogeographical or regional processes(Ricklefs 1987) rather than by the underlying mechanisms of community assembly described above. Disentangling the relative importance that species pool size(Kraft et al.2011)or community assembly mechanisms have on determining β-diversity at different scales along elevation gradients in the tropics is paramount for developing robust forest conservation plans capable of maintaining diversity(Lomolino 2001;Rahbek 2005).

The spatial scale at which vegetation studies are developed is a key factor that can strongly influence β-diversity gradients (Stier et al. 2016). The concept of scale involves two factors:i)extent,the geographical area where comparisons are made; and ii) grain size, the unit of measurement at which data are collected or aggregated for analysis(Whittaker et al. 2001). In a fixed extent,a variation in grain size implies a variation in the sampled relative species abundances and,subsequently,in the spatial patterns of aggregation (Crawley and Harral 2002). Directly related to βdiversity, when the spatial grain size of local communities increases, species present in the regional species pool will be better represented, generally lending to a decline in βdiversity(Barton et al.2013).Along an elevational gradient,the use of 0.1-ha plots with grain sizes of 0.01-ha has been widely used to assess and detect fine-grained environmental variation effects on determining β-diversity at a local scale(Kraft et al.2011;Mori et al.2013;Tello et al.2015).However, in species-rich communities, smaller grain sizes may lead to the undersampling of individuals, an issue that can artificially enhance β-diversity (Condit et al. 2005). Comparative studies of β-diversity at contrasting grain sizes along elevational gradients are needed to help disentangle the extent to which either sampling effects or community assembly mechanisms shape β-diversity patterns.

Along elevational gradients, another largely unexplored issue pertains to the likely effect that different diameter at breast height (DBH) cut-off sizes can have in β-diversity assessments(Mori et al.2013).Overall,reducing the minimum size, or DBH, of the sampled individuals increases the community size, potentially increasing floristic diversity measurements as well (Stier et al. 2016). In tropical mountains, the most popular DBH cut-off size utilized to assess changes in β-diversity along elevational gradients are individuals with DBHs varying from ≥2.5 cm (Kraft et al. 2011; Myers et al. 2013; Tello et al. 2015) to ≥10 cm DBH (Girardin et al. 2014). However, none of these studies have evaluated the likely comparative effect that tree cut-off size variation can have on shaping β-diversity.The sampling effect of keeping the grain size constant and decreasing the DBH cut-off will cause a change in species relative abundance; and whereby this difference in abundance may lead to changes in the extent to which underlying ecological mechanisms can explain the overall pattern of diversity (Powell et al. 2011; Chase and Knight 2013). In other words, sampling not only has a potential effect on the diversity patterns, but also on our ability to identify the underlying community assembly mechanisms that drive these observed patterns.For example,in tropical lowlands, several studies have proposed that enhancing community size by including smaller individuals (e.g.shrubs and juveniles)may lead to a higher influence of deterministic processes,such as soil fertility,on defining species sorting (Duque et al. 2002; Comita et al. 2007).Understanding the effect of different tree cut-off sizes in determining the magnitude of β-deviation at a local scale along elevational gradients will help to distinguish sampling constructs from true ecological signals. This is essential in helping researchers to identify the underlying drivers of species distribution and forest function in the tropical Andean mountains.

In order to identify the likely influence of local community assembly mechanisms on shaping β-diversity along elevational gradients, we first need to determine whether β-diversity deviates from null (stochastic) processes (Kraft et al. 2011). Null models help to disentangle ecological assembly mechanisms by quantifying random processes in the ecological community and making comparisons among regions with different species pool sizes possible (Chase and Myers 2011). A positive standardized difference between the observed β-diversity and the expected β-diversity obtained from a null model divided by the standard deviation of the null model (defined here as β-deviation), indicates a higher β-diversity than expected by chance due to the influence of local processes that cause an aggregated non-random spatial pattern of species distribution (Mori et al. 2013; Tello et al. 2015). However, a positive and systematic increase of β-deviation along the elevational gradient, after removing sampling effects and differences in the size of species pools among sites, is not enough and fails to identify the underlying community assembly mechanism (e.g species sorting or dispersal limitation) responsible for an aggregated non-random pattern along the whole elevational gradient (i.e. Tello et al. 2015).Mirroring the magnitude of the operating species assembly mechanisms found along the latitudinal gradient (Myers et al. 2013), we might expect the relative importance of biological processes, such as dispersal limitation, to decrease with elevation;an opposing effect to species sorting,which can be positively correlated with elevation.

In this study, we employed a nested sampling design using a series of 15 1-ha plots scattered in wet forests located in northwestern Colombia, where the Andean mountain ranges end,to examine the role that species pool size, grain size and tree cut-off size played in determining β-diversity along elevational gradients. For this study, we had three main hypotheses: i) under the assumption that local variation in species composition primarily depends on the size of the species pool, we do not expect any significant relationships between β-deviation and elevation to occur after controlling for the species pool (Kraft et al.2011). In contrast, if ecological mechanisms (e.g. species sorting) determine a non-random spatial species distribution, the variation on β-deviation may show a systematic change with elevation as a result of the harsh conditions imposed by highlands (after Tello et al. 2015). ii) The increase of grain size within a fixed extent increases the floristic similarity among samples (hereafter grain size hypothesis),and thus,decreases β-diversity.We expect the magnitude of the relationship between elevation and βdeviation (the slope of the line) to decrease with the increase of grain size at a local scale along the elevational gradient. iii) The reduction of the selected tree cut-off size will increase the local community size and will reduce the compositional differences between samples.We also would expect a reduction in the β-deviation of each plot along the elevational gradient.

Methods

Study area

The study area was located in the northwest region of Colombia between 5°50′ and 8°61′ North and 74°61′and 77°33′ West. This region encompasses a highly variable elevational gradient in terms of its topography, climate, and soils. The study was conducted using data collected from 15 permanent 1-ha (100 m × 100 m) forest inventory plots which were established between 2006 and 2010. The permanent plots were established across a large geographic area that covers approximately 64,000 km2, mostly within the province of Antioquia (Fig. 1) and span an elevation gradient of 50 to 2950 m asl.The average distance between plots was 160.5 km (ranging=26.1-419.5 km). The Andean region in Colombia contains only approximately 34% of its original natural cover primarily due to historical deforestation (Duque et al. 2014; Cabrera et al. 2019). Thus, at least in some of the surveyed locations, we expected to find some previous human disturbances, specifically in the El Bagre,Carepa and Necoclí plots (Fig. 1), which were located in small forest fragments (≈50 ha). These plots may have experienced human disturbance and elevated tree mortality along forest edges (Duque et al. 2015).

Plot censuses

In each 1-ha plot, all shrubs, trees, palms, and tree ferns with a diameter at breast height (DBH)≥10 cm (hereafter “l(fā)arge trees”) were mapped, tagged, and measured.Additionally, all of the plants with a DBH ≥1 cm (hereafter“all trees”)were also mapped,tagged and measured in a 40 m×40 m subplot(1600 m2)located near the center of each plot(Additional file 1:Figure S1).Voucher specimens were collected for each potentially unique species in each plot.We collected vouchers in all cases where there was any doubt as to whether an individual plant was the same species as another individual that was already collected within the same plot. Taxonomic identifications were made by comparing the specimens with herbarium material and with the help of specialists for some plant groups. Vouchers are kept at the University of Antioquia’s Herbarium (HUA).The plants that could not be identified to the species level were classified into morphospecies based on differences in the morphology of their vegetative characters. Approximately 3.5% of individuals were excluded from the analysis due to low-quality vouchers resulting from a lack of clear botanical characters,earlier stages of development,or incorrect enumeration. In total, we identified 26,222 individuals,112 families,428 genera and 1707 morphospecies.

Sampling effects DBH cut-off and species pool size effect

We divided the dataset into three DBH cut-off sizes: i)large trees: represented by all individuals with a DBH ≥10 cm tallied in the entire 100 m × 100 m plots (1-ha); ii)small trees: represented by all individuals with a 1 cm ≤DBH ＜10 cm, which were measured only in the 40 m × 40 m subplot inserted within the 1-ha plot (Additional file 1: Figure S1); iii) all trees: represented by all individuals with a DBH ≥1 cm tallied in the 40 m × 40 m subplot (0.16-ha) described above. In order to assess the effect of species pool size for each one of the tree DBH cut-off sizes employed to generate our three sampling communities (large, small and all trees), we used the species richness corresponding to each data set. For large trees, we used the species richness from each 1-ha plot but only including trees with a DBH ≥10 cm. For the small and all trees categories, we used their respective species richness from each 0.16-ha plot (40 m×40 m) (see Table 1).

Fig.1 Location of 15 1-ha plots in Antioquia on a regional map(inset)to show its location within Colombia.The elevation range of the plots are presented in grayscale:white color for plots located between 0 and 1000 m asl;gray for plots located between 1000 and 2000 m asl,and black for those located between 2000 and 3000 m asl

Grain size effect

The grain size hypothesis was assessed by employing three different grain sizes. For large trees, we used 10 m × 10 m (0.01-ha), 20 m × 20 m (0.04-ha) and 50 m × 50 m (0.25-ha). The grain size used to analyze the influence of the spatial scale for small and all trees were 5 m×5 m (0.0025-ha), 10 m × 10 m (0.01-ha) and 20 m × 20 m (0.04-ha). The differences in the spatial grain size among large versus small and all trees are due to individuals with a DBH ≥1 cm were only measured in the 40 m ×40 m subplot.

Environmental features

The elevation of each plot was calculated using a GPS.Each elevation point corresponds to the 0,0 point located in the lower-left corner of each plot along the gradient (Additional file 1: Figure S1). Samples of the soil A horizon (mineral soil after removing the organic layer) from five points in each 20 m×20 m quadrat were collected (N=25 composite samples per 1-ha plot). At each point, a 500 g soil sample was taken from a depth of 10-30 cm; the five samples from each quadrat were then combined, and a 500 g composite sample was taken and air-dried after removing macroscopic organic matter. pH, Ca, Mg and K concentrations were analyzed at the Biogeochemical Analysis Laboratory at the National University of Colombia in Medellín. Exchangeable Ca,Mg, and K were extracted with 1 mol·L-1ammonium acetate and analyzed using atomic-absorption. Soil pH was measured in water as one-part soil to two parts water. Other soil cations, such as N and P, were not measured due to logistical constraints of sampling at this spatial resolution and scale.

We used geostatistical methods to obtain spatial predictions of soil variables at spatial scales smaller than 20 m×20 m (5 m×5 m and 10 m×10 m). We first computedempirical variograms to test the likely spatial structure of each soil variable (pH, Ca, Mg, and K) within the 1-ha plot. The variograms for the four variables did not show any spatial significant trend.Therefore, we used a bilinear interpolation method based on resampled soil data to obtain values of soil variables at different grain sizes in each plot. This method employs the distance-weighted average of the nearest pixel values to estimate the values of no measured points (Hijmans 2016). We calculated soil variables at the 50 m×50 m grain size using the mean of the soil variables at the 20 m×20 m scale. Spatial analyses were conducted using the geoR (Ribeiro and Diggle 2001)and raster(Hijmans 2016)packages.

Table 1 Description and location of the 15 1-ha permanent plots in the northwest of Colombia.Latitude (North) and Longitude(West) are presented in geographical coordinates (degrees). N:total number of individuals. S: species richness. The columns of 0.16-ha contain the information about N(number of individuals) and S(species number) by different DBH cut-off size in the 40 m×40 m subplot inside the plot.The column of 1-ha has information about N and S for the large trees in the whole plot

Estimations of β-diversity

We calculated the observed β-diversity (BDobs) based on abundance data (Legendre and Gallagher 2001; De Cáceres et al. 2012). Taking into account all living trees by species in each one of the plots, for every grain size,we built a matrix (X=[xij]) with dimension n×p (quadrat × species), where X is the community matrix of each plot and xijcontains the number of individuals of species j in the quadrat (grain) i (De Cáceres et al. 2012). For each matrix X=[xij], β-diversity was estimated in two steps. First, we transformed the abundances of each species by grain size using the Hellinger transformation.This transformation consists in standardize the abundance of each species by rows. It means, to standardize the abundance of each species by the total abundance of the site (in this case, species by grain), in each plot.Then, the square root of these values is taken (Legendre and Gallagher 2001). Thus, data set express species abundance as square-root transformed proportionate abundance in each grain by site (Jones et al. 2008). The Hellinger transformation is given by:

where Yijis the transformed matrix, xijis the value of species j in site i,k is the species index and p is the number of species in a given grain with row and column indices i and j(Tan et al.2017).The Hellinger transformation standardizes species abundance and reduces the weight of the most abundant species in the analysis. The use of the Hellinger transformation makes community compositional data containing many zeros (“double zero”) suitable for analysis by linear methods (Legendre and Gallagher 2001; Legendre 2007). Secondly, we estimated BDobsas the variance of Y(De Cáceres et al.2012),which is calculated as follows:

where SS(Y) is the sum of squares and n is the number of quadrats. BDobsis 0 when all quadrants have exactly the same composition and 1 when they do not share any species.

Null model

We used a null model to quantify the extent to which the variation in the size of species pool (different species number due to the DBH cut-off size) and scale (different grain size) account for variation in β-diversity (Kraft et al.2011). The species pool for large, small and all trees was defined as the observed number of species in either the 1-ha or the 0.16-ha plots (after Kraft et al. 2011). The null model randomizes the location of trees among grains within the plot, creating communities that vary in relation to the location of individuals, but fixing the community size(number of individuals),and thus,the observed relative species abundance of each species pool (Tello et al. 2015).This null model removes the local ecological mechanism that creates non-random patterns,such as aggregation and intraspecific co-occurrence (De Cáceres et al. 2012). The Hellinger transformation is then applied to the randomized matrix and expected β-diversity(BDexp)is calculated using the formula presented above.This process is repeated 1000 times per plot,for each grain size,and for each predefined DBH cut-off size. The BDexpis calculated as the mean of 1000 iterations of the null model.

β-deviation (BDdev) was defined as the standardized effect size (SES) calculated using the difference between BDobsand BDexpdivided by the standard deviation of the frequency distribution of the null model (SDexp).

Positive values in the slope of the variation between BDdevalong elevational gradients indicate a significant effect of community assembly mechanisms on determining the rate of change in species composition at local scales(Chase and Myers 2011; Tello et al. 2015). Contrarily,values of the slope of the variation in BDdevalong elevational gradients non-significantly different from zero (0)are primarily due to sampling effects that come up along with the variation in the size of the species pool (Kraft et al.2011).

Data analysis

We used linear mixed regression models (LMM; Zuur et al. 2009) to identify the main determinants of change in BDobs, BDexp, and BDdevalong the elevational gradient. Variables included in the LMM as fixed effects were:grain size, size of the species pool, elevation (m asl) and soil heterogeneity. Soils heterogeneity was assessed for each grain size using the interpolated values from 20 m×20 m subplots described above. To represent soils heterogeneity at a local scale,we used the variance of the subplot scores on the first axis of a principal component analysis(PCA). PCA was applied to pH, Ca, Mg, and K concentrations. PCA analyses were performed for each grain size and DBH cut-off size (Additional file 1:Methods). Soils heterogeneity was modeled as a continuous variable. Finally, plot identity (or plot name) was included as a random effect to control for particular conditions of each site (Zuur et al. 2009). The interaction term between grain size and elevation was included to directly assess the combined effect of these variables on shaping the β-diversity(BDobs,BDexp,and BDdev).

In LMMs, the marginal explained variation (R2marginal) is associated with fixed effects, while the conditional explained variation (R2conditional) associated with random effects. Because individuals with DBH ≥1 cm and with 1 cm ≤DBH ＜10 cm were not sampled at the 50 m×50 m scale, we were unable to include the three tree size categories in the same model. Therefore,separate models were used for large trees, small and all trees. The best model for each DBH cut-off size was chosen using the backward stepwise model selection based on the Akaike information criterion (AIC) (Crawley 2007). In order to assess the likely spatial autocorrelation in our models, we extracted the residuals for each model (BDobs, BDexp, and BDdev, for large, small and all trees), separating them by grain size, and assigning the respective spatial coordinate to each one. Then, we estimated a semi-variogram based on 100 draws to define an envelope for the significance of the observed spatial structure of the residuals. This analysis was performed with the geoR package (Ribeiro and Diggle 2001).

All analyses were performed in R 3.3.0 (Core Team 2016).

Results

Elevation and species pool

As we expected, BDobsand BDexpdecreased with elevation independent of the grain size and DBH cut-off size(Fig. 2). In contrast, BDdevincrease with elevation, also in all grain sizes, regardless of the DBH cut-off size (Fig.2). After controlling for the regional species pool effect,BDdevstill showed an increase with elevation. Overall,the standardized local BDdevincreased from lowlands to highlands, which suggests a differential effect from the underlying species assembly mechanism in accordance to elevation.

Grain size

Both BDobsand BDexpdecrease with grain size independent of the tree DBH cut-off size (Fig. 2). The slopes among grain size, or the relationship BDdev-elevation,were significantly different for large trees, but small and all trees did not show any significant difference among grains (Additional file 1: Figure S2).

Fig.2 Observed(BDobs), expected(BDexp),and standardized (BDdev)patterns of variation of β-diversity along the elevational gradient.β-deviation(BDdev)is taken as (BDobs- BDexp)/SDexp. Upper panel(A,B,C): large trees(DBH ≥10 cm). Middle panel(D,E,F):small trees(1 ≤BDH ＜10 cm)and Lower panel(G,H,I): all trees(DBH ≥1 cm).Large trees are taken into account in an area of 1-ha.Small and all the trees are taken into account in 0.16-ha plot

Determinants of local scale changes in tree β-diversity along the elevational gradient

According to the LMMs, the BDobswas significantly associated with grain size, the size of the species pool and elevation for the three size-classes employed (large trees,small trees, all trees). The interaction between grain size and elevation was only significant for large trees. The BDexpwas significantly associated with grain size and elevation for the three DBH cut-off size employed, while the size of species pool was significant for large and all the trees but only marginally significant for small trees.The BDdevwas significantly associated with grain size for all the three DBH cut-off size. The interaction between grain size and elevation was significant for large and small trees, but not for all the trees. Finally, the marginal explained variation (R2marginal) by the models was almost always the same than that explained by the conditional variation (R2conditional) for observed and expected β-diversity and for BDdevin large trees. However, the marginal and conditional explained variation for BDdevfor small and all trees had differences, which indicates greater relative importance of random effect for the last two tree sizes (Table 2). Model residuals showed no evidence of spatial autocorrelation (Additional file 1: Figs. S3-S5).

Discussion Sampling effects

In this study, we assessed three hypotheses regarding the influence of sampling effects (size of species pool, grain size, and tree cut-off size) on the variation of local βdiversity along elevational gradients in the northern region of the Andean mountains of Colombia. Overall, we found that observed and expected β-diversity decreased with elevation, but that the standardized β-deviation followed an increasing trend with elevation after controlling for the effect of species pool size. The systematic increase in the β-deviation with elevation was independent of the grain size employed, indicating that alternative underlying community assembly mechanisms had a significant role in shaping tree β-diversity along this elevational gradient. Our finding contradicts the claim of sampling effects due to the species pool size as the key determinant of changes in β-diversity (sensu Kraft et al.2011). Therefore, our results emphasize the importance that different community assembly mechanisms have on shaping the observed decrease in local β-diversity along elevational gradients in tropical forests (Mori et al. 2013;Tello et al. 2015), rejecting our first hypothesis.

Following some studies on tree β-diversity along latitudinal gradients (De Cáceres et al. 2012; Sreekar et al.2018), our second hypothesis predicted and confirmed a decrease in both the observed and expected tree βdiversity with the increase in grain size along an elevational gradient. Regarding the β-deviation, our findings were dependent on the DBH cut-off tree size as predicted by the third hypothesis, similar to other studies along elevational gradients (Mori et al. 2013). Mori et al.(2013) claimed that the overall β-diversity decreases in response to the DBH cut-off size, contrary to βdeviation. Therefore, for large trees (DBH ≥10 cm), we accept the hypothesis that changes in grain size have a significant effect on the assessment of the standardized β-deviation, and conclude that the larger the grain size,the lower the observed β-diversity, but the higher the βdeviation. In other words, especially for large trees, and along elevational gradients, the probability of detecting the influence of community assembly mechanisms increase positively at larger grain sizes (Fig. 2). A likely explanation for this pattern could be that large trees are those that survived self-thinning and their spatial distribution, at smaller spatial scales (e.g. 0.04-ha), are more random than at larger scales, which indicates that the degree of aggregation does not vary much at such small grain sizes.

Table 2 Results from the best-fit linear mixed models for large (＞10 cm DBH),small(1 cm ≤BDH ＜10 cm) and all trees (DBH ＞1 cm). BDobs: observed β-diversity.BDexp: expected β-diversity. BDdev: β-deviation (BDobs - BDexp)/SDexp. Conditional R2 takes into account both fixed and random effects to measure the goodness of adjustment and prediction power,while marginal R2 only has the fixed effects part. NS ＞0.05, *p ≤0.05, ** p ≤0.01,***p ≤0.001

Table 2 Results from the best-fit linear mixed models for large (＞10 cm DBH),small(1 cm ≤BDH ＜10 cm) and all trees (DBH ＞1 cm). BDobs: observed β-diversity.BDexp: expected β-diversity. BDdev:β-deviation (BDobs -BDexp)/SDexp. Conditional R2 takes into account both fixed and random effects to measure the goodness of adjustment and prediction power,while marginal R2 only has the fixed effects part. NS ＞0.05, *p ≤0.05, ** p ≤0.01,***p ≤0.001 (Continued)

When assessing the β-deviation for the small and all individuals size classes (DBH ≥1 cm), the interaction between grain size and elevation included in the LMMs was significant for small trees but not for all trees. This contrasting result, stemming from similarly nested datasets (see Table 1), hampers our capacity to make conclusions as to the effect of grain size on the local βdeviation for the small and all individuals along the elevational gradient. In fact, when using an independent Analysis of Covariance (ANCOVA) to evaluate the grain size- elevation interaction term, only large trees were significant(Additional file 1:Table S1;Figure S2).The low sampling size (4) used to assess tree β-diversity at the largest grain size (4) may be a reason for the high variance observed when we included individuals with DBH ≥1 cm. In the Andean mountains, the lack of sampling schemes of plots ≥1-ha that include individuals with DBHs ≥1 cm,such as those available for tropical lowlands (i.e Anderson-Teixeira et al.2015),prevents us from concluding about the expected trend of the β-deviation at larger grain sizes along the elevational gradient in tropical forests.

Tree community assembly mechanisms along the elevational gradient

The increase of β-deviation in relation to elevation indicates that in colder regions, the extent to which species assembly mechanisms operate is higher compared to warmer areas. One important conclusion to note is that low temperatures may impose constraints to plant establishment and functioning,and play a key role in determining species distribution (Kitayama and Aiba 2002;Girardin et al. 2014). For example, changes in species composition could be associated with changes in species richness along elevational gradients in very diverse understory families,such as Rubiaceae(r=-0.58,p=0.02).

Soil variation has been shown to be a key community assembly mechanism which shapes species sorting at local scales in some tropical forests (Russo et al. 2005; John et al. 2007). However, in this study, we did not find soil variation to be significantly associated with the local βdeviation along the elevational gradient. This result did not support the idea of an increase in plant habitatassociation of juveniles and shrubs (Duque et al. 2002;Comita et al.2007; Fortunel et al. 2016).Nonetheless,our soil variation index focuses primarily focuses base content,hindering our ability to understand the likely influence of other very important soil cations,such as P and N,which,in tropical lowland forests (Condit et al. 2013), have been identified as key elements for tree species distribution.Furthermore,soil sampling was only carried out at the 20 m×20 m scale, which might have obscured processes operating at smaller spatial scales. Additional studies testing the likely influence of topographic and edaphic variables,not considered here,will shed new insights on the still unanswered question about the extent to which environmental filtering locally shapes species sorting, and thus,the gradient of β-diversity at local scales along elevational gradients in tropical forests.

The lack of significance of soil variation on shaping species sorting implies that other community assembly mechanisms, rather than environmental filtering, are likely drivering the observed change in β-diversity at a local scale with elevation. Mirroring the latitudinal gradient (Myers et al.2013),a systematic decrease in the importance of dispersal limitation(sensu Hubbell 2001)with elevation seems the first likely alternative assembly mechanism to explain the increase in β-deviation observed in this study. Another possible explanation for the positive deviations of βdiversity is the hypothetic positive increase of densitydependence with the size of the species pool (Lamanna et al.2017),which suggests that the stronger the conspecific and heterospecific the negative dependence is, the higher the diversity,but the weaker the influence of environmental filtering and niche partitioning. A decrease of species competition but an increase of species facilitation in highlands,due to the adverse conditions imposed by low temperatures on the ecosystem functioning and survival capacity of plants(Coyle et al. 2014), could also promote the observed increase of β-deviation with elevation observed in our study.

One likely factor not assessed here that could have influenced the pattern of variation in local β-diversity is the expected biotic homogenization caused by forest disturbance(Karp et al.2012;Solar et al.2015).The high fragmentation and historical degradation of the tropical Andes(Armenteras et al. 2013), could have caused some of our sites to display a lower local β-diversity than under undisturbed conditions. In mountainous ecosystems, we expect the steep terrain at the highest mountain peaks to limit site access and act as a shield against human disturbances(Spracklen and Righelato 2014),thus generating a higher biotic homogenization in lowlands than in highlands.Indeed,the plots located in the smallest forest fragments (Carepa,Necoclí and El Bagre;see methods),were all located in lowlands. However, the systematic decline in the observed βdiversity (BDobs) does not support the hypothesis of biotic homogenization as a major cause of the observed pattern.For example, we did not find statistical differences (unpaired t-test) when comparing the β-deviation between the three sites located in the smallest forest fragments, which we assumed were exposed to higher disturbances, and the rest of the plots located in lowlands(＜1000 m asl).This result was a generalized outcome for any grain size for both large trees(50 m×50 m:p=0.79;20 m×20 m:p=0.82;10 m×10 m: p=0.42) and small trees (20 m×20 m: p=0.92;10 m×10 m:p=0.78;5 m×5 m:p=0.64).

Methodological remarks

First, for large trees, the LMMs selected species pool size(species richness) as a significant variable to explain the variation of the β-deviation with elevation (Table 2).This finding indicates that the applied null-model did not,in some cases, entirely and effectively remove the influence of the size of the species pool.Understanding the effect that changes in the shape of the species abundance distribution models have on determining the β-diversity along elevational gradients is still under debate (e.g. Qian et al. 2013). However, it could be seen as an alternative way to analyze the effect from changes in community size.Second,the absence of plots ≥1-ha that include small individuals in the Andean mountains prevents the use of sampling sizes along the elevational gradient which are large enough to properly assess the grain size and cut-off size hypotheses together in this complex ecosystem. Although our study is the first attempt in the Andean mountains to test the species pool hypothesis using plots ≥0.1 ha, our results were based on very few replicates of the largest grain sizes and need to be seen as preliminary evidence of an expected pattern rather than a conclusive view. To truly understand the pattern of β-diversity variation in mountainous tropical forests,it appears we need to transition towards a new generation of larger forest sampling schemes(e.g Garzon-Lopez et al.2014;Duque et al.2017;Sreekar et al.2018)that goes beyond the valuable heritage left by A.L.Gentry.Such a big challenge should be a priority in the tropical Andes,where the availability of information is much more scarce than in their Amazon lowland counterparts(Feeley 2015).

Conclusion

We determined that the effect of the grain size, species pool size and tree cut-off size, are paramount to identify the underlying processes that shape species assembly of tree communities. Our findings suggest that grain size enlargement and the inclusion of small size classes can help improve our ability to identify the extent to which the species assembly mechanisms shape the patterns of local β-diversity change along elevational gradients in tropical ecosystems.However,in future field campaigns that aim to assess tree local β-diversity along the elevational gradient in tropical forest inventories, we need to evaluate the limitation of the relatively small plot size employed so far.Overall, our study emphasizes the scale-dependent nature of β-diversity assessments. It showcases the advantage to decreasing the tree cut-off size and increasing the plot size in forest inventories (De Cáceres et al. 2012; Barton et al.2013; Sreekar et al. 2018) to improve our understanding about the likely response of tree diversity to global change in tropical mountain ecosystems.

Supplementary information

Supplementary information accompanies this paper at https://doi.org/10.1186/s40663-020-0214-y.

Additional file 1 Forest Ecosystems. Figure S1: Graphical representation of each one of the plots.Methods: Schematic description of the analytical procedure employed to extract the soils data set.Figure S2:Results post hoc ANCOVA analysis using “Tukey” test,comparisons among each slope in the linear mixed models. Figure S3: Mixed linear model validation for large trees using variograms with model residuals using Pearson method and geographical coordinates of the plots.Figure S4: Mixed linear model validation for small trees using variograms with model residuals using Pearson method and geographical coordinates of the plots. Figure S5:Mixed linear model validation for all trees using variograms with model residuals using Pearson method and geographical coordinates of the plots.Table S1 Analysis of covariance(ANCOVA). Comparison of slopes between grain size and elevation for the β-deviation and for all of the DBH cut-off sizes.

Abbreviations

BDdev: Deviation in Beta diversity; BDexp: Beta diversity expected by the null model; BDobs: Beta diversity observed

Acknowledgments

We like to thank Jonathan Myers for his reading and comments to the manuscript. Miquel De Cáceres provided code in R to evaluate the null models. Elysa Cameron kindly revised the language of the manuscript. We are indebted to two anonymous reviewers and the handling editor for all their wonderful comments that certainly helped to improve the quality oof the manuscript.

Authors’contributions

AD designed the study. JA M-V and SG analyzed the data. JA M-V and AD wrote the paper. All authors jointly discussed and agreed to the final version.

Authors’information

JA M-V is a doctoral student in Université du Quebec a Montreal, Canada.Her main research focuses on the ecological mechanisms that drive taxonomic and functional diversity along environmental gradients. S G-C is a doctoral student at the National University of Colombia-Medellín. His interests focus on the ecological and evolutionary drivers of tree communities.AD: is an Associated Professor at the National University of Colombia. His research focuses on tree community structure and function.

Funding

The project “Bosques Andinos”was funding by Helvetas Swiss development organization and developed by Medellín Botanical Garden “Joaquín Antonio Uribe”.

Ethics approval and consent to participate

The subject has no ethic risk.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interest.

Author details

1Departamento de Ciencias Forestales,Universidad Nacional de Colombia sede Medellín, Cra. 65 #59a-110, Medellín, Colombia.2Département des Sciences Biologiques, Université du Québec à Montréal, 141, Avenue du Président-Kennedy, Montreal, Canada.

Received: 7 May 2019 Accepted: 5 January 2020