Assemblage of forest communities in subtropical montane forests of western Mexico

Abstract

Functional diversity indices have been used to differentiate the relative contribution of stochastic and deterministic processes that modulate the assemblage of communities; however, knowledge regarding the relative contribution of assemblage mechanisms in forest communities is scarce. We analysed the assembly mechanisms driving forest assemblages along a topographic gradient at two spatial scales (1000 m² and 3000 m²) for three different forest types from subtropical mountain forests (pine-oak, mixed pine-broadleaf and broadleaf forest) in western Mexico, using null models of multi-trait indices. The forest structure differed along the topographic gradient. Upper slopes were dominated by Pinus douglasiana with an importance value index (IVI) of 0.8, while 10 and 13 tree species were required in the middle and lower slopes, respectively, to reach the same IVI. The results support the idea that the subtropical montane forest of western Mexico is a mosaic of communities, when analysed at a scale of 1000 m² the forest assembly was mainly explained by stochastic processes, while analysis at the scale of 3000 m² showed that functional convergence of species were the main mechanisms of the assemblage of the pine-oak forest communities due to an abiotic stressful environment.

Introduction

Biological communities are characterised by species composition and relative abundance, which are produced by historical and random processes such as the Pleistocene glaciations, that occurred over wide geographic (thousands of kilometres) and temporal (millions of years) scales and shaped the group of species in a given region (HilleRisLambers et al. 2012). On a smaller local scale (tens of kilometres), community species composition has been explained by the effect of stochastic processes, such as ecological drift and the species’ seed dispersal capacity (Hubbell, 2001; Norden, 2014), and by deterministic processes such as environmental filters, these last can be abiotic, defined by the physicochemical properties of sites, or biotic, such as the competition for resources. Environmental filters limit the establishment of species and define their dominance within each community, increasing their functional convergence or niche differentiation (MacArthur and Levins, 1967; HilleRisLambers et al. 2012; Laughlin, 2014). Discerning the relative contribution of stochastic and deterministic mechanisms to structuring community has been a focal point in the field of community ecology, and functional diversity (FD) studies can provide quantitative evidence for achieving this goal (Cornwell and Ackerly, 2009; Spasojevic and Suding, 2012). FD describes the distribution of species in a functional space structured by axes, which represent the functional traits in a community (Rosenfeld, 2002; Mason et al. 2005; Suárez-Castro et al. 2022). The functional traits are the anatomical, phenological, morphological and physiological characteristics that affect the performance of an organism (i.e., fitness), and the sensitivity of these traits to environmental variation reveals the effect of environmental and evolutionary pressures on populations (Violle et al. 2007; Díaz et al. 2016; Maynard et al. 2022). Patterns of species abundance and functional traits reflect the degree of functional differentiation among species (i.e., biotic interactions) at a given locality as well as the effect of the environment on their functioning (i.e., abiotic filters) (Mason et al. 2007; Villeger et al. 2008).

Theoretically, when the set of species’ functional traits exhibits relatively low diversity and is associated with high resistance to environmental stress (functional underdispersion), the community assemblage is considered to be modulated by abiotic (deterministic) filters (HilleRisLambers et al. 2012). For example, tree species that establish in dry, fire-prone environments tend to develop thick bark and small leaves (Laughlin, 2014). On the other hand, when the combination of species’ functional traits comprising a community is highly diverse (functional overdispersion), this indicates a high degree of ecological niche differentiation, promoted by interspecific interactions, within the community assemblage (MacArthur and Levins, 1967; Mason et al. 2007; Villeger et al. 2008; HilleRisLambers et al. 2012). For example, along gradients of light and water availability in forests, interspecific differences in leaf functional traits related to shade and drought tolerance have shown to be associated with a niche differentiation (Engelbrecht et al. 2007; Poorter, 2009). To determine the importance of environmental filters (biotic and abiotic) as modulators of community assemblage, multidimensional FD indices have been proposed, including functional richness (FRic), functional evenness (FEve) and functional divergence (FDiv) (Mouchet et al. 2010). Although there is a correlation between these indices and species richness, this is removed by calculating the standardised effect size (SES) through null models (Swenson, 2014). So, the combined of these three SES-FD indices are considered an appropriate analytical tool to discriminate community assemblage rules, because of their multi-trait nature, they are complementary and independent, and uninfluenced by species richness (Mason et al. 2005; Villeger et al. 2008; Mouchet et al. 2010; Mason et al. 2013; Götzenberger et al. 2016).

Few studies in forest ecosystems have employed FD indices to empirically explain the mechanisms underlying community assemblage (Zakharova et al. 2019), surely because their conclusions about the drivers of community assemblage patterns are often difficult to compare, probably because they do not always consider the influence of spatial scale, they employ indices that do not capture different diversity components, or the functional traits employed do not capture the niche overlap among species (Münkemüller et al. 2020). For instance, there are studies that employed the same multi-trait indices, but each one analysed FD at different spatial scales: 200 m² plots (Sanaphre-Villanueva et al. 2016) and 1000 m² plots (Bhaskar et al. 2014), however, showed completely different assembly patterns for the same ecosystem (tropical dry forest), and this probably is an indication that overlooking emergent properties of scale in communities (Münkemüller et al. 2020). On the other hand, there are studies that analysed FD patterns at different spatial scales (accumulated), but using a single-trait indices (Kraft and Ackerly, 2010; Zhang et al. 2018), both studies showed a sign of functional convergence at the small scale and loss of same at the large scale, but this pattern may be an effect of sensitive to the rise of species richness and not necessarily a proxies for niche overlap (Mouchet et al. 2010; Münkemüller et al. 2020).

Las Joyas Research Station (LJRS) is located in west-central Mexico, within the Sierra de Manantlán Biosphere Reserve (SMBR). The vegetation within the LJRS comprises a mosaic of pine-oak forests, associated with convex landforms and high slopes where fire has been a frequent disturbance, broadleaf forests, occurring in hollows and humid ravines that function as shelters from fire, and intermediate mixed pine-broadleaf forests (mixed pine-broadleaf forests) (Jardel-Peláez et al. 2004b). Those forests occurring at elevations between 1600 and 2200 m asl, support the coexistence of tree species of both Holarctic and Neotropical affinity, and are referred to as subtropical mountain forests (Jardel-Peláez et al. 2004b; Saldaña-Acosta et al. 2008). Also, these forests have a protective function of the headwaters of watersheds (Corlett and Hughes, 2015; Galicia et al. 2015). So, these forests’ conservation is especially critical.

A previous study reported the convergence of functional traits in diverse tree lineages of different biogeographic origins (Saldaña-Acosta et al. 2008). Other studies have related the distribution pattern of different forest communities to slope topography and soil properties (water content and nutrient concentrations). The high and convex parts of the slopes, with soils that are relatively drier and with lower nutrient concentrations, are dominated by pine-oak species, while broadleaf forest species are scarce or absent. This latter community dominates the low and concave parts of the slopes in soils with higher fertility and greater moisture content, thereby supporting wider diversity of broadleaf species that interact and compete for resources. This suggests that the assemblage of tree species in the pine-oak forest is the result of abiotic filters, while that of the broadleaf forest is due to niche differentiation (Saldaña-Acosta, 2001; Castillo, 2019; Quintero-Gradilla et al. 2020). However, the spatial separation between the two forest types is not discrete since transition zones are present, in which pines and broadleaf species constitute a mixed forest (Jardel-Peláez et al. 2004a, 2004b; Saldaña-Acosta et al. 2008). Analysis of the functional traits of the trees in these forest communities allows us to elucidate the underlying mechanisms that contribute to species assemblages along a topographic gradient in the subtropical mountain forest of the Sierra de Manantlán.

Our aim was to discern the relative contribution of deterministic and stochastic processes in forest community assemblage along a topographic gradient, which probably is driving three communities’ composition: pine-oak forest, broadleaf forest and mixed pine-broadleaf forest of the SMBR. We measured leave functional traits related to structural (i.e., petiole stem-specific density and petiole length), biophysical (i.e., stomata density, specific leaf area and leaf dry matter content) and hydrological constraints (i.e., leaf water potential and stomatal conductance) (Silvertown et al. 2015), whose correlation is used to help detect morpho-physiological trade-offs in the utilisation of resources (i.e., water, light and nutrients) (Wright et al. 2004; Klein, 2014; Diaz et al. 2016). Three FD indices were estimated that are considered the most accurate in terms of identifying the mechanisms driving community assemblage at different spatial scales: FRic, FEve and FDiv (Mason et al. 2005; Villeger et al. 2008; Mouchet et al. 2010; Mason et al. 2013; Carmona et al. 2016). Our hypothesis is that if the mosaic of tree communities (pine-oak forest, broadleaf forest and mixed forest) along a topographic gradient is due to environmental filters as suggested by Castillo (2019), Quintero-Gradilla et al. (2020) and Saldaña-Acosta (2001), abiotic filters would be expected to explain the forest assembly of pine-oak, biotic interactions explain the assembly in the broadleaf forest, and stochastic processes (seed dispersion) explain the mixed pine and broadleaf forest assembly (Fig. 1). Furthermore, if FD patterns vary at the scale at which they are analysed, this indicates that scale has an influence on community assembly mechanisms (Münkemüller et al. 2020), which would support the community mosaic hypothesis open connected by dispersal (i.e., meta-communities) (Leibold et al. 2004).

Figure 1. Conceptual framework of the study. (a) Observation and description of forest community structure patterns along a topographic gradients, (b) research question, (c) hypothesis and (d) predictions of the relative contribution of species assemblage mechanisms in the communities, derived from functional diversity patterns (Liu and Wang, 2018; Ortega-Martínez et al. 2020). Bars show the standardised effect size value (SES-Value), estimated from the observed functional indices and null models. Community 1 (com 1) represents an assemblage of species that is mainly explained by abiotic factors (functional underdispersion), communities 2 and 3 (com 2, com 3) represent an assemblage of species that is explained by seed dispersal, and community 4 (com 4) represents an assemblage of species with high functional differentiation that responds to biotic factors (functional overdispersion). Flow diagram adapted from Münkemüler et al. (2020).

Materials and methods

Study area and sampling design

The study was conducted at the LJRS of the University of Guadalajara. The LJRS has an area of 1257 ha and is located at coordinates 19°14’49’’–19°37’30’’ N and 104°14’49’’–104°18’16 “W, and between 1550 and 2250 m asl, within one of the core zones of the SMBR.

The climate is temperate-warm humid with a summer rainfall regime. According to the Köppen system corresponds to Cw subtropical with a summer rainfall regime, or to subtropical lower montane humid forest of the Holdridge Life Zone System, characterised by hot and rainy summers and moderate winters in which frosts occur but not snowfall. The annual ambient temperature was 15.8 ± 0.1°C, the annual precipitation was 1992 ± 57 mm and the water balance (potential evapotranspiration to precipitation ratio) was 0.5 (SE; 1990–2022, LJRS weather station records, 1950 m asl). January was the coldest month (13.2 °C) and May the warmest (17.8°C); 63% of the precipitation occurs between July and September, and the driest period (relative air humidity between 46 and 60%) occurs between March and May.

The vegetation is mixed pine-oak and pine-broadleaf forests, moist broadleaf forest (mesófilo de montaña, sensu Rzedowski, 1978) and secondary scrub established on old abandoned agricultural fields (Jardel-Peláez et al. 2004b). Most of the pine-oak and pine-broadleaf forests were logged for commercial timber production between 1946 and 1980, and fires have been recorded with a return interval (Weibull median probability) of 5–9 years; logging, cattle grazing and forest fires were excluded after the establishment of the SMBR in 1987 (Jardel-Peláez, 2018).

Within the LJRS, three independent slopes were selected, referred to here by their local names: 1) Escobedo (Esc), 2) Triguito (Tri) and 3) Charco de los Perros (CP). These sites were located at least 1 km apart and shared similar topography, with a slope of between 30% in the upper part and 15% in the lower part, and no recorded fire events in the last 30 years. All had a northeastern aspect, conferring a drier environment compared to the slopes with a southwestern aspect, facing the Pacific Ocean. The higher part of each slope is dominated by pine-oak forest, the lower part by broadleaf forest and the middle part by mixed forest, consisting of a combination of pine and broadleaf (Appendix 1).

Nine 1000 m² (20 × 50 m) plots were established following the method of Jardel-Peláez et al. (2004b). With one plot in each forest type (pine-oak forest, mixed pine-broadleaf forest and broadleaf forest) on each of the three slopes. The distance between the upper (∼2100 m asl) and lower (∼1800 m asl) part of the slope was ∼300 m. Within each plot, all trees of diameter at breast height (DBH, 130 cm above ground level) ≥5 cm were censused. Each tree was labelled and identified to species level, and its height was recorded with a Haga clinometer (Haga GmbH & Co. KG., Germany) and DBH measured with a diameter tape.

A climatic microstation (H21-USB Onset HOBO, Bourne, MA, USA) with sensors was established in the centre of each plot, at a soil depth of 20 cm, to measure soil water content (S-SMD-M005 Onset HOBO) and soil temperature (S-TMB-M002 Onset HOBO). Systematic hourly recording of soil water content and soil temperature was conducted in May 2021. During the sampling month, the environment reached an average vapour pressure deficit of 1.24 ± 0.01 SE kPa, with a single rainfall event of 19.1 mm.

Functional traits

The dry season is the most suitable time to make functional comparisons between species since interspecific differences are more clearly manifested in that period. In May 2021, between 10:00 and 14:00 h, the branch with the leaves most exposed to the sun was collected from three individuals of each tree species in each plot. Telescopic scissors were used to collect branches from trees of up to 10 m in height. For trees taller than 10 m, a pneumatic pistol was used. Three healthy, fully expanded and hardened leaves were selected from each branch along with a portion of the branch stem of ∼10 cm in length and ∼2 cm in diameter. The following parameters were estimated directly in the field during sampling: (1) mid-day water potential (ψ, MPa) using a Scholander Chamber (Model 600, PMS); (2) mid-day stomatal conductance (gs, mmol m⁻² s⁻¹) using a Porometer (AP4, DELTA-T DEVICES). The material collected was transported to the LJRS laboratory following the protocols of Cornelissen et al. (2003) and Perez-Harguindeguy et al. (2013), and were analysed: (3) leaf dry matter content (LDMC, mg g⁻¹); (4) specific leaf area (SLA, cm² g⁻¹); (5) stem specific density (SSD, mg mm⁻³); (6) petiole length (L-Pet., cm), and (7) stomatal density (SD, mm²). To calculate the FD, the seven functional traits were used.

Forest structure along a topographic gradient of slopes

To describe the community structure along a topographic gradient of slopes, the importance value index (IVI) of each species was calculated from the average of the relative values of density, basal area and frequency (McCune and Grace, 2002) for the pine-oak forest, broadleaf forest and mixed pine-broadleaf forest.

Statistical analysis

FD was measured with the three multidimensional indices proposed by Villeger et al. (2008): (1) FRic, which represents the volume in N-dimensional space (as many dimensions as traits), (2) FEve, which provides the degree of functional similarity between rare and abundant species, and (3) FDiv, which generates the degree of functional differentiation among abundant species. The calculation of FD is based on a principal coordinate analysis, which reduces dimensionality and generates new axes with which the functional space defined by species attributes can be calculated (Swenson, 2014). To eliminate the influence of trait units on FD estimation, functional attribute values were standardised with a mean of 0 and variance of 1 (Villeger et al. 2008). For FD estimation, the FD package (Laliberté et al. 2014) was used. The three indices were estimated at the spatial scales of 1000 m² (plot) and 3000 m² (cumulative of each community).

Null models were employed to assess whether the observed values differ from the null expectations (Swenson, 2014; Götzenberger et al. 2016; Sanaphre-Villanueva et al. 2016; Ortega-Martínez et al. 2020). To conduct this analysis, the SES of each index was computed, representing the deviation of the observed value from the mean of the null distribution divided by the standard deviation of the null distribution (Swenson, 2014). The SES values were estimated through constructed null models with 999 iterations following Swenson (2014), which randomises the names of species in the functional traits matrix while preserving the fixed richness of species and abundances in the matrix of communities. The randomisation process was applied to the entire set of pine-oak, mixed and broadleaf forest types, assuming equivalence in their respective species regional pools. Positive SES values indicate that the observed result is greater than that expected, while negative values indicate that the observed result is less than that expected (Swenson, 2014). The biological interpretation of the values is as follows: (1) if the SES-Value < −1.96, abiotic filters drive the community assemblage, (2) if the SES-Value >1.96, biotic interactions modulate the assemblage, and (3) if the SES-Value is between −1. 96 and 1.96, community assemblage is the result of stochastic processes (Götzenberger et al. 2016; Liu and Wang, 2018; Ortega-Martínez et al. 2020) (see Fig. 1).

The null model was constructed using the package Picante (Kembel et al. 2010). A significant difference with a p < 0.05 is considered when SES values fall outside the range −1.96 to 1.96, assuming a normal distribution of deviation (Swenson, 2014; Liu and Wang, 2018; Ortega-Martínez et al. 2020).

To quantify the differences in soil properties along the topographic gradient, environmental data recorded at mid-day in May 2021, which corresponds to the driest month of the year, were used. Soil water content and soil temperature data did not fulfil the assumptions of normality, so comparisons among the three communities were performed with a Kruskal–Wallis test, for which the package FSA was used (Ogle and Ogle, 2017).

All the analyses described were performed using R software (R Core Team, 2021), and the packages used are in R software.

Results

Assemblage of communities in the plots (scale of 1000 m²)

Just in plot 1 (pine-oak forest) the SES-FEve was 2.01, and it was also the plot with the driest soil and highest temperature (Appendix 3), whilst in the other eight plots, SES values of the three FD indices did not go beyond the range of values of −1.96 and 1.96 (Fig. 2). Thus, in general, stochastic processes primarily determined the assemblage of the communities at the plot scale.

Figure 2. Standardised effect size (SES) of functional richness (FRic), functional evenness (FEve) and functional divergence (FDiv) of three plots per forest type: pine-oak forest, mixed pine-broadleaf forest and broadleaf forest (1000 m² scale).

Assemblage of communities by forest type (scale of 3000 m²)

The FRic in the pine-oak forest was lower than expected at −2.18, while the FD of the broadleaf forest and mixed pine-broadleaf did not go beyond the range of values of −1.96 and 1.96 (Fig. 3a). The value of SES-FRic suggests that the assemblages of the pine-oak community was due to deterministic mechanisms.

Figure 3. Community assembly of pine-oak forest, mixed pine-broadleaf forest, and broadleaf forest communities at a scale of 3000 m² by forest type, and the soil parameters in these forest types. (a) Standardised effect size (SES) of functional richness (FRic), functional divergence (FDiv) and functional evenness (FEve) of the species assemblage of trees with DBH ≥5 cm, along a topographic gradient of three slopes in a subtropical zone in western Mexico, the horizontal red lines show the value 1.96 and −1.96. (b) Importance value index for the tree community in the upper (pine-oak forest), middle (mixed pine-broadleaf forest) and lower (broadleaf forest) parts of the slope. (c) Average value of soil parameters recorded in May 2021 (±standard deviation).

In addition, the pine-oak forest recorded the highest soil temperature (H = 159.39, df = 2, p < 0.001) and the lowest moisture content (H = 77.04, df = 2, p < 0.001), compared to the broadleaf forest, which was cooler and wetter (Fig. 3c).

Structure of the communities

On calculating the IVI in each community, in the pine-oak forest, Pinus douglasiana presented an IVI of 0.8. In the mixed pine-broadleaf forest, it was necessary to sum the IVI values of Pinus douglasiana, Zinowiewia concinna, Styrax ramirezii, Magnolia iltisiana, Tilia caroliniana subsp. occidentalis, Clethra fragrans, Cornus disciflora, Carpinus caroliniana, Myrsine juergensenii and Symplococarpon purpusii to reach an equivalent value. On the other hand, in the broadleaf forest, 13 species were required to reach the IVI value of 0.8. Some of the species with the highest contribution were Carpinus caroliniana, Symplococarpon purpusii, Magnolia iltisiana, Fraxinus uhdei and Symplocos citrea (Fig. 3b). Moreover, 10 tree species were observed to be shared among the communities (Aiouea glossophylla, Citharexylum mocinnoi, Clethra fragrans, Cornus disciflora, Fraxinus uhdei, Magnolia iltisiana, Myrsine juerguensenii, Persea hintonii, Symplocos citrea and Zinowiewia concinna) (Appendix 2).

Discussion

The results confirmed that, along the topographic gradient of the subtropical montane forest of western Mexico, forest community structures and soil properties (water content and temperature) vary. In the upper and convex parts of the slopes, soil water content was the lowest and soil temperature the highest, and these areas are mainly dominated by Pinus douglasiana, which reached an importance value of 0.8. In the lower and concaved parts of the slopes, where soil water content was the highest and soil temperature the lowest, 13 broadleaf tree species were required to reach an importance value of 0.8. Finally, in the middle part of the slope, where soil temperature was intermediate and soil water content similar to that in the upper parts, 10 tree species between broadleaf and P. douglasiana were required to reach an importance value of 0.8. This is a pattern already reported in the study area by several authors (Saldaña-Acosta, 2001; Castillo, 2019; Quintero-Gradilla et al. 2020).

On analysing the contribution of deterministic and stochastic mechanisms to the tree species assemblage of three community types, from a multi-trait perspective of the species and considering two spatial scales: 1000 m² scale (plot) and 3000 m² (accumulated plots), community assemblage patterns from the FD indices were sensitive to the spatial scale at which they were analysed, particularly in the pine-oak forest, as reported in previous studies (Cavender-Bares et al. 2006; Carboni et al. 2013; Trisos et al. 2014; Zhang et al. 2018). However, in contrast to that reported for a subtropical forest by Zhang et al. (2018), who used single-trait indices (RaoQ), the present study showed that, at a scale of 1000 m², stochastic processes prevail over deterministic ones, while at a scale of 3000 m², deterministic processes become more significant to elucidating the composition of pine-oak communities. We suggest that the differences in ecological patterns between the two studies could be related to the fact that the single-trait indices (RaoQ) lose power to detect community assembly rules when species richness increases (Mouchet et al. 2010). And when the sampling area is expanded, the species richness increases.

Stochastic mechanisms best explained assemblage at the plot scale (1000 m²)

Surprisingly, just plot 1 in the pine-oak forest registered a SES-FEve value of 2.01, we believed that this result is due to the fact that the indices estimate the functional similarity of the rare species with the abundant ones (Villéger et al. 2008), and in the case of plot 1, only three tree species were recorded: Pinus douglasiana with 22 individuals, Pinus oocarpa with 2 individuals and Quercus obtusata with 1 individual, and functionally trade-off of the 3 species could be similar in trait value stem-specific density (∼0.6 mg mm⁻³), leaf-dry matter content (∼400 mg g⁻¹), specific leaf area (∼5 cm² g⁻¹), stomatal conductance (∼55 mmol m⁻² s⁻¹), water potential (∼−2 MPa) and long petiole (∼25 mm) (Appendix 2), as well, plot 1 was the one that recorded the most stressful soil conditions, with the lowest soil water content and highest soil temperature (Appendix 3). Some important aspects of this result are that the functional trait similitude between species on plot 1 could be related to the stressful soil conditions that would drive hydrological niche segregation (Silvertown et al. 2015), and probably one way to test this hypothesis is through an analysis of the relationship between the plant functional traits and soil properties (Anderegg, 2023). On the other hand, communities with high FEve may be more resilient to disturbance processes that alter the availability of resources (Oliver et al. 2015) since the loss of a rare species can be functionally replaced by another species, thus maintaining ecosystem functioning. Finally, the SES-FD difference between communities (plots) suggests functional heterogeneity in the topographic gradient of the ecosystem, highlighting the relevance of the protection of the mosaic of communities and their connectivity for the conservation of functions delivered by ecosystems.

On the other hand, the values of richness, evenness and FDiv, at the 1000 m² plot scale, did not exceed the −1.96 to 1.96 range in eight plots of the three forest types, suggesting that stochastic processes contributed more to the species assemblage of the communities at this scale.

This result differed from that reported in studies that employed single-trait indices in tropical rainforest at a scale of 10,000 m² (Kraft and Ackerly, 2010), subtropical forest at a scale of 400 m² (Zhang et al. 2018), as well as studies with multi-trait indices in tropical dry forest in plots of 200 m² (Sanaphre-Villanueva et al. 2016), but were similar to that reported for tropical dry forests by Bhaskar et al. (2014), who also conducted their analysis in 1000 m² plots. However, by not detecting environmental filters, these authors suggested that forest community assemblage is being influenced by forest management history, and they barely mention the stochastic mechanisms and the scale effect. These conclusions may be affected by what Münkemüller et al. (2020) call ‘common pitfalls’, since on the one hand unidimensional indices only consider a single component of diversity, and by not analysing the relationship between functional traits, they fail to capture the niche overlap among species and, on the other hand, community assemblage processes can occur at different spatial and temporal scales.

We can conclude from our results that, with exception of plot 1, at the scale of 1000 m² (plot) the species assemblage in subtropical montane forest communities was dictated by stochastic mechanisms such as seed dispersal, species demographic randomness or species ecological equivalence. Although it was not possible to identify the contribution of potential stochastic mechanisms, it was recorded that seed dispersal is mainly modulated by dispersers such as bats and birds in the study area (Rost et al. 2015). Moreover, some of the species with high importance values, including Pinus douglasiana, Clethra fragrans, Carpinus caroliniana and Fraxinus uhdei, produce anemochorous seeds, which allows them to benefit from wide dispersal (Cuevas et al. 2021). Thus, the two main mechanisms of seed dispersal in the study area may be promoters of randomisation in the colonisation of the entire diversity that comprises the region (Beckman et al. 2020). It is also known that the floors of the pine-oak forest and mixed pine-broadleaf forest harbour a seed bank that requires fire for germination (Zuloaga-Aguilar et al. 2016) and that species regeneration and succession are mainly controlled by their fire tolerance (Quintero-Gradilla et al. 2020). Tree functional traits related to water exchange between the plant (branch and leaf), soil and atmosphere were quantified in this study. It can therefore be considered that the ecological equivalence detected with the FD indices focused on the functional compensations of the leaf to capture resources (i.e., the leaf economic spectrum). Although it is also possible that the ecological equivalence of species was an artefact of the plot size in which the community assemblage was analysed (Hubbell, 2001), since at the 1000 m² scale it was not possible to detect niche differentiation or functional overlap among adult tree species.

Deterministic mechanisms drive the assemblage at the accumulated plot scale (3000 m²)

At the 3000 m² scale, the hypothesis that forest community assemblage is modulated by abiotic filters in the pine-oak forest and by stochastic processes (seed dispersal) in the mixed pine-broadleaf forest was supported but was not supported by biotic filters in the broadleaf forest.

The value of FRic of the pine-oak forest at the scale of 3000 m² was lower than expected (SES-FRic = −2.74). This indicates that the functional volume is small and formed by a few strategies with low dispersion in the functional space, suggesting functional convergence. Convergence of functional attributes among various tree species in the study area has been proposed previously (Saldaña-Acosta et al. 2008) but not tested appropriately. The results of this study identified that functional convergence occurred on the upper parts of the slopes, where species potentially tolerant of the drier environment of these areas along the topographic gradient (see Fig. 3c), with consequent effects on the assemblage of the pine-oak forest community. As at the plot scale, the community assemblage of the mixed forest and broadleaf forest analysed at the 3000 m² scale is mainly explained by stochastic processes, suggesting that the richness, divergence and functional similarity in the two vegetation type is mainly explained by the randomness of seed dispersal or ecological equivalence.

Differences in FD patterns in pine-oak forest between the 1000 m² and 3000 m² scales suggest that tree communities in subtropical montane forests are comprised of a mosaic of communities (Leibold et al. 2004), the niche differentiation or functional convergence of which is detected when the scale is enlarged. Interaction with the environment may have effects at the 3000 m² scale, while seed dispersal may have effects at the 1000 m² scale, regulating the assemblage of communities. Nonetheless, to test the meta-community hypothesis, detailed information on the demographics of the populations is needed.

On the other hand, the FD analysis of the accumulated three pine-oak plots (3000 m²) may be affected by the environmental variability of each plot. An example of this is the coefficient of variation of water content in the soil of the pine-oak forest (see Fig. 3c: standard deviation/mean × 100), which reached 37%, indicating a notable variation between plots. This coincides with what was reported by Zhang et al. 2018 who point out that the greater the environmental variability, there is greater the sensitivity to identifying the functional convergence of tree communities. At this point, we consider that not only the spatial scale can affect assembly patterns but also the gradients of environmental variation, since in essence it is the environment that modulates the functional response of the species.

Finally, our findings highlight the value of multi-trait indices as a good proxy for niche overlap (the more functional traits, the more precision in the description of the niche of the species), detect community assembly rules and are independent of species richness. This last is important due to the community assembly rules could change through different spatial scales and with the rise of area sampling the richness increase. We do not rule out that the functional patterns in the broadleaf forest are due to the scale at which the analysis was made, and it is probably necessary to test the community assemblage hypothesis over a broader spatial and environmental gradients. Moreover, the effects of the most common disturbance processes in the study area, such as fire, on the assemblage of communities should be evaluated, and functional attributes specifically related to fire tolerance should be quantified.