Material and methods

Study area. We carried out the study at LTBS (18° 35' N, 95° 05' W), which is located in a volcanic complex of the Upper Tertiary and Middle Quaternary (< 7 Ma), that is composed of basaltic andesite and basalt (^{Verma et al. 1993}) and bordered by lava flows from the Holocene (^{Nelson & González-Caver 1992}). The landscape in LTBS is topographically heterogeneous, with slopes ranging from 2 to 45°, covering an area of 640 ha, where the elevation ranges between 150 and 700 m asl (^{Ibarra-Manríquez et al. 1997}). The annual average temperature is 24.5 °C and the mean annual precipitation (MAP) is 3,928 mm, with a dry season from March to May, where it rains less than 10 % of the total annual precipitation (^{SMN 2010}). The most abundant tree species are Damburneya ambigens (S.F. Blake) Trofimov, Guarea glabra Vahl, Omphalea oleifera Hemsl., and Pseudolmedia glabrata (Liebm.) C.C. Berg (^{Ibarra-Manríquez et al. 1997}). Additionally, 13 species of palms have been recorded; the most abundant include A. mexicanum, Chamaedorea alternans H. Wendl., C. pinnatifrons (Jacq.) Oerst., and C. tepejilote Liebm. Although there are few palm species, their high density makes them an important element of this TRF (^{Ibarra-Manríquez et al. 1997}, ^{Martínez-Ramos et al. 2016}). Despite being surrounded by cattle pasture, the protected status of the LTBS has protected the forest from fires, logging, and plant extractions for more than 50 years, such that it is representative of conserved TRF in the region (Figure 1).

Figure 1 (A) Study area and location of plots on the land units (High Altitude Slopes Strongly Inclined, HSSI; Low Altitude Slopes Strongly Inclined, LSSI; Lava Flows, LF; Piedmont, PD; Steeply Dissected Slopes, SDS). (B) Satellite image of the study area on Google Earth. Data: Google ^©2021, Maxar Technologies/TerraMetrics. 

Environmental heterogeneity characterization. We stratified the study area based on the five LUs delimited by ^{Navarrete-Segueda et al. (2021)}. These discrete units are based on the systematic integration of environmental components such as geology, relief, and soil (^{Zinck et al. 2016}). This approach has been used to stratify the environmental mosaic and the potential habitats in the study area (^{Navarrete-Segueda et al. 2021}) as well as in other Neotropical rainforests (^{Clark et al. 1995}, ^{Brown et al. 2013}, ^{Denslow et al. 2019}). Four of the five LUs were associated with cinder cones (late Tertiary and early Quaternary), that differ in inclination, elevation, microclimate, and their soils, which were formed from volcanic ash. The fifth landscape unit corresponds to a landscape associated with Holocene volcanic events and is characterized by shallow stony soils (Table 1). The selected LUs are representative nearly of 75 % of the landforms of the LTBS (640 ha). In each of the LUs, three plots of 20 × 50 m (0.1 ha) were established, covering a total sample area of 1.5 ha. Plots on steep slopes were oriented transverse to the slope to minimize the effect of topographic variation. Estimated inter‐site distances ranged between 400 and 4,550 m.

In each plot, we described soil profiles, based on ^{Schoeneberger et al. (2012)}, within the rooting depth (^{Siebe et al. 1996}). The available water holding capacity (AWHC) of each horizon was estimated and interpreted according to ^{Eckelmann et al. (2005)} and ^{Siebe et al. (1996)}. This variable integrates the soil texture, organic matter content, bulk density, percentage of coarse fragments, and the thickness of each horizon estimated in the field. In addition, texture, rooting depth, and percentage of coarse fragments were selected as indicators of nutrient storage and soil fertility (^{Bünemann et al. 2018}). We also integrated elevation and slope as relief variables since they have been shown to influence the spatial variation in vegetation structure and diversity (^{Baldeck et al. 2013}). To obtain the abiotic components related to relief, we parameterized the topographic data of slope and elevation, based on a digital elevation model (DEM) constructed using contour lines to 10 m resolution of ^{INEGI (2010)}, processed in ArcMap 10.1 (ESRI®). We extracted the average value and standard deviation of these variables for each plot.

We generated climatic data because in the Sierra de Los Tuxtlas there is a low density of meteorological stations (^{Gutiérrez-García & Ricker 2011}). To accurately interpolate and strengthen the climate surface, we included weather data from the surrounding area. We employed eight meteorological stations, with elevations ranging from 4 to 1,200 m. Yearly average precipitation (mm yr^-1) and temperature (°C) were generated with Inverse Distance Weighting (IDW) interpolation through ArcMap (ESRI^®). IDW interpolation is a mathematical (deterministic) function that assumes that closer values are more related than farther values (^{Chen & Liu 2012}).

In each 0.1 ha plot (20 × 50 m), we recorded tree individuals with a diameter at breast height (dbh) ≥ 10 cm and identified them to the species level. Using dbh data, we calculated the biomass per tree species (t dry mass 0.1 ha^-1) based on the allometric equation [1] applied to tropical regions with a MAP between 2,000- and 4,000-mm year^-1 (^{Brown 1997}, ^{Rügnitz et al. 2008}):

To estimate AGB per ha we used the equation [2]:

Where ∑at is the sum of the dry mass of all trees in the plot (AGB_t) of 0.1 ha and 1000 is the factor to convert kg into Mg.

The use of specific allometric equations based on forest type has proven to be an efficient method that captures 95 % of biomass variation (^{Brown 2002}, ^{Berenguer et al. 2015}). In each plot of 20 × 50 m (0.1 ha). We recorded and identified all the palms with a height greater than 1 m, except the adult individuals of Reinhardtia gracilis (H. Wendl.) Drude ex Dammer, which are smaller than 1 m. However, it has been shown that this species responds to soil conditions (^{Mendoza & Franco 1998}) and is an important element of the palm community in the study area (^{Ibarra-Manriquez & Sinaca-Colín 1995}). Therefore, individuals of R. gracilis were sampled only when flowers or fruits were observed.

Data analysis. We carried out generalized linear models (GLM) and multiple comparisons of means with post hoc Tukey tests to assess differences in abiotic and biotic components among LUs using ‘multcomp’ package in R (^{R Core Team 2022}). We apply error distribution of Poisson to count data and Gaussian to continuous data (after normality test), and log and identity link functions, respectively (^{Crawley 2007}). We applied a principal component analysis (PCA) to describe major gradients in the soil and topography to reduce the number and collinearity of variables (^{Laurance et al. 1999}). Therefore, we expected that PCA would allow us to detect and classify meaningful variables related to the palm community parameters. The analysis was performed with the "scale" argument to standardize the data, results were centered by mean and, variables were presented by their correlations. The PCA was performed using the dudi.pca function of package ˈade4ˈ (^{Dray & Dufour 2007}).

Data of individuals and richness were used to obtain mean (± S.E.) values of density (number of palms in 0.1 ha) and effective number of species (⁰D and ¹D) for each landscape unit, following ^{Jost (2006)}. We carried out GLM and multiple comparisons of means to assess differences in ⁰D, ¹D, and density of palms (0.1 ha) among LUs.

We constructed rank-abundance curves to assess changes in the dominant species of palms among the LUs, lumping the data of the three plots per unit. The curves were constructed following ^{Magurran (2004)}. Linear models were fit for rank- abundance curves and differences between units were tested using an analysis of covariance (ANCOVA) considering as an explanatory variable the rank (rank; as numeric) and landscape unit (Landscape unit; as categorical). To compare the slope of rank-abundance curves according to LUs, abundance of species was transformed to Log₁₀. Differences among slopes were evaluated between LUs (^{Izsák 2006}).

The effects of environmental heterogeneity on palm community were analyzed by a permutational multivariate analysis of variance (PERMANOVA) (^{Anderson & Braak 2003}). This analysis was performed with the adonis2 function of ˈveganˈ package of R. The environmental data were normalized with the scale function. PERMANOVA is especially useful for partitioning multivariate data in response to complex designs (e.g., LUs) (^{Anderson 2017}). Values of P were obtained using 9,999 permutations. Non-metric multidimensional scaling (NMDS) techniques were used to visualize the PERMANOVA results, which allowed us to examine the multivariate relationship of palm species (presence and abundance) with biophysical factors among LUs. The NMDS analysis were executed using the function metaMDS (ˈveganˈ package) in R (^{Oksanen et al. 2019}). We implemented NMDS with Bray-Curtis dissimilarities to density and with the Jaccard method to the presence-absence data. Since NMDS uses rank order information, this analysis represents a highly flexible method to explore the proximities and resemblances of the species composition and the structure of palm community in relation whit environmental heterogeneity (^{McCune & Grace 2002}).

Results

Environmental heterogeneity among the LUs. The GLM results showed that LUs are different in several of the analyzed environmental factors (Table 1). The slope presented a high variation. The units with low inclination were Piedmont (PD) and Lava Flows (LF) (< 11° each). These units differed significantly from Low Altitude Slopes Strongly Inclined (LSSI) and Steeply Dissected Slopes (SDS) (> 19° slope) (F_{4, 10} = 56.86, P < 0.000). The PD unit was found at lower elevations than LF (Table 1; Figure 1). High Altitude Slopes Strongly Inclined (HSSI) and SDS have a higher elevation, corresponding to the highest parts of the volcanoes. The LF unit, whose soils were formed by basaltic rock flows, had high surface stoniness (> 85 %), leading to soils with shallow rooting depth and low water storage (Table 1). In the other units, because the soils were formed from volcanic ash, the rooting depth and the soil water storage were greater, while the presence of coarse fragments was low. There were no significant differences in relief aspect and the density, biomass, and tree diversity among the LUs (F_4,10= 36844, P = 0.301; F_4,10 = 8.115, P = 0.0875; F_4,10 = 3330134039, P = 0.78; F_4,10 = 99.68, P = 0.133, respectively) (Table 1).

The first axis of the PCA (42 % of the variance) showed a gradient associated with the AWHC for plants and the coarse fragments in the soil (Table 2; Figure S1, Supplementary material). Coarse fragments in the soil are related to the LF unit, while, at the opposite extreme are the LUs recognized in soils formed from volcanic ash, with less coarse fragments content and larger stores of water in the soil, mainly SDS (Table 1; Table 2). On the second axis (26 % variance), the LUs are differentiated by a gradient determined by the elevation and the annual precipitation (Table 2). At one end of this axis are the PD unit and one plot corresponding to LSSI, located at low elevation and with high precipitation, separated from HSSI that had the opposite trend in these two environmental variables (Figure S1, Supplementary material).

Palm composition and density in the context of environmental heterogeneity. A total of 10 palm species and 2,761 individuals were recorded. No significant difference was found in species richness (⁰D) (F_{4, 10}= 0.81, P = 0.92) and the number of common species (¹D) (F_{4, 10} = 7.63, P = 0.153) of palms among the LUs (Table S1, Supplementary material). We found that palm density differed among LUs, except between LSSI and SDS (F_{4, 10} = 61.98, P < 0.000). Palm density was highest in the LF unit and lowest in HSSI (Table S1, Supplementary material), which was clearly due to a low density of A. mexicanum and C. alternans in the latter unit (Figure 2; Table S2, Supplementary material).

Figure 2 Species-rank curves based on species abundance for palm community present in the land units of Los Tuxtlas Tropical Biology Station reserve, SE Mexico. In the X-axis, species are ranked in decreasing order of their abundance, while in the Y-axis, the number of palms is indicated. HSSI (High Altitude Slopes Strongly Inclined), LSSI (Low Altitude Slopes Strongly Inclined), LF (Lava Flows), PD (Piedmont), SDS (Steeply Dissected Slopes). 

Rank-abundance curves showed that the most abundant species in all the LUs was A. mexicanum (Figure 2). The second most abundant species differed among LUs; C. alternans in LF and PD, C. ernesti-augusti H. Wendl SDS and LSSI, R. gracilis in HSSI. We found significant differences among rank-abundance curves of LUs with contrasting soil properties (Table S3, Supplementary material). The comparison showed that HSSI curve differs from LF, and PD curves (F₁= 0.455, P < 0.000 and F₁= 0.3143, P = 0.013, respectively). Simultaneously, we found difference between SDS and LSSI curves (F₁= 0.207, P = 0.036). We did not find significative differences among other LUs curves (Table S3, Supplementary material).

There were distinct distribution patterns in palm assemblages among the LUs. Results from the PERMANOVA showed significant effects of slope of relief (P = 0.009), elevation (P = 0.002), MAP (P = 0.04), soil coarse fragments (P = 0.002), and tree D¹ (P = 0.008; Table 3) on palm species presence. These results clustered the plots with low inclination (PD) and separated them from the LUs with abrupt relief and higher elevation, as visualized by NMDS (stress value 0.15) in Figure 3A. Species that responded to this gradient are Bactris mexicana Mart., which was registered mainly in low-inclination plots, while C. ernesti-augusti and C. elatior Mart. were found mostly in plots with abrupt relief (Table S2, Supplementary material). We found that 33 % of the variation in palm species presence was explained by topography since relief elevation and slope explain 21 and 12 %, respectively (Table 3). Coarse fragments in soil and tree diversity (¹D) accounted for 25 and 13 % of the variation, respectively.

Figure 3 Non-metric multidimensional scaling (NMDS) of (A) presence and (B) abundance of the palm community in the study landscape. Species presence and abundance in 20 × 50 m plots were used. Both species scores and plot scores for the first two ordination axes are plotted. Topography, soil, and tree community attributes in the plots were fitted to ordination to test the relationship between species composition and abundance, and environmental heterogeneity. MAP (Mean annual precipitation), HSSI (High Altitude Slopes Strongly Inclined), LSSI (Low Altitude Slopes Strongly Inclined), LF (Lava Flows), PD (Piedmont), SDS (Steeply Dissected Slopes). Initial letter of the genera: A (Astrocaryum), B (Bactris), C (Chamaedorea), and R (Reinhardtia). 

Results from the PERMANOVA showed significant effects of relief slope (P = 0.034), elevation (P = 0.03), soil coarse fragments (P = 0.003), and tree D⁰ (P = 0.028; Table 3) on palm species abundance. The NMDS in Figure 3B (stress value < 0.2) shows that, except for LSSI, the plots were grouped consistently with the LUs. Coarse fragments in the soil explained 19 % of the palm density variation (Table 3). C. tepejilote had the highest density in plots with high coarse fragments content (LF), while R. gracilis was most abundant in plots with deeper soils (and HSSI) (Table S2, Supplementary material). Similar to the palm species presence patterns, the density of B. mexicana, C. ernesti-augusti, and C. elatior responded to elevation and slope, which accounted for 24 % of variation in palm species density among the plots. Tree richness (D⁰) accounted for 13 % of density variation, with mainly, R. gracilis and C. alternans at the opposite extremes concerning tree richness and soil rooting depth interaction (Figure 3B).

Discussion

Environmental heterogeneity and the LUs. The delimited LUs based on geomorphology and soil captured the spatial variation of the abiotic component, detecting two main gradients: i) the coarse fragments and the rooting depth in the soil and ii) elevation and precipitation. The first gradient separates deep soils derived from volcanic ash (SDS, PD, LSSI, and HSSI) from shallow soils formed from basalt (LF) (Figure S1, Supplementary material). This gradient is an indicator of the water store and nutrients in the soil (^{Tetegan et al. 2015}, ^{Drobnik et al. 2018}). Variation in elevation and soil properties related to topographic position have been reported as predictors of palm density and diversity in other tropical forests (^{Clark et al. 1995}, ^{Cámara-Leret et al. 2017}). We found that elevation and its relationship with precipitation influenced the spatial distribution of the palm species. In this regard, ^{Perrigo et al. (2020)} and ^{Svenning (2001a)} indicate that in landscapes with significant elevation variation (e.g., > 100 m difference), as occurs in our study area, precipitation may vary even over small distances. The variation in this climatic variable was particularly relevant in LF, with lower soil water available for plants, since palm individuals could experience water stress during three months in which evapotranspiration exceeds rainfall. This finding may indicate that precipitation variations due to climate change can have the potential to affect the distribution of palm species at the landscape scale, especially those with lower soil AWHC.

Palm composition and density in the context of environmental heterogeneity. The ordination analysis showed that LUs are a good stratification approach to identify the effect of environmental heterogeneity on palm community composition in Los Tuxtlas reserve. Site distributions in the ordination reflected palm species distribution between units. This result is in line with other studies that proved the efficiency of sampling based on soil and topography to capture and explain spatial variation of tree, understory plants, and palm communities that are structured by species responses to environmental factors in volcanic Neotropical landscapes (^{Clark et al. 1995}, ^{Denslow et al. 2019}, ^{Navarrete Segueda et al. 2021}).

The altitude and the slope had a strong effect on the presence of B. mexicana and C. ernesti-augusti and was evident from the decrease in the density of A. mexicanum and C. alternans in the higher elevation units (SDS and HSSI; Figure 3B). This result coincides with those of ^{Cepeda-Cornejo & Dirzo (2010)}, for C. alternans and C. elatior that showed specific affinity with the LUs having lower precipitation and well-drained soils, either due to a low percentage of clay or high stoniness (LSSI and LF; Figure 3A). According to studies of ^{Clark et al. (1995)}, ^{Svenning (2001a)}, and ^{Muscarella et al. (2019)}, this result may be because soil water conditions are regulator of the composition of TRF palm species. Variables that differentiated PD from LF-elevation and coarse fragments in soil (Figure S1, Supplementary material) influence the distribution of C. ernesti-augusti and C. tepejilote, which are practically absent in PD and relatively abundant in LF. ^{Bacon & Bailey (2006)} and ^{Miceli-Méndez et al. (2013)} also reported higher density of C. tepejilote and C. ernesti-augusti in sites with shallow and stony soils, respectively. ^{Renninger et al. (2013)} suggest that palms are more efficient in resource allocation and distribution of aerial and root biomass in the face of constraints such as shallow soils with low AWHC. This strategy has been shown by C. tepejilote (^{Oyama & Dirzo 1988}), which may have an advantage over tree species, which have deeper roots. The allocation of biomass to the stem and roots has been reported as an important indicator of competitive interactions in the face of restrictive factors (^{Qi et al. 2019}).

Four species were recorded in all the LUs (A. mexicanum, C. alternans, C. pinnatifrons and R. gracilis) (Figure 2). The wide mesoscale distribution of C. pinnatifrons has been reported by ^{Svenning (2001b)}, who found that the density of this species is associated with deep, well-drained soils. This result contrasts with our findings, since soil depth was not a factor that influenced the density of this species. In the case of A. mexicanum, its presence of has been widely reported in this TRF (^{Popma et al. 1988}, ^{Ibarra-Manríquez et al. 1997}, ^{Hernández-Ruedas et al. 2018}) and is usually related to its high seed production capacity in the face of disturbances (^{Martínez-Ramos et al. 2016}). However, the hyperdominance of A. mexicanum decreased significantly in the highest elevation and rooting depth units (SDS and HSSI) (Figure 2), contrasting with the proposal of ^{Martínez-Ramos et al. (2016)}. This divergence in results is likely attributable to the fact that the plots of ^{Martínez-Ramos et al. (2016)} were all located on the PD unit, which excluded the effect of the elevation on this species. Our results support the argument of ^{Arroyo-Rodríguez & Melo (2016)} that it is important evaluate palm populations throughout the landscape in the Sierra de Los Tuxtlas, since the soil and topographic position have a clear relationship with A. mexicanum distribution throughout the analyzed forest.

The low density of A. mexicanum in SDS could also be a consequence of the increase in light in the lower strata of the forest. According to ^{Clark et al. (1996)} and ^{Alves et al. (2010)}, the lateral incidence of light at sites with abrupt relief can be equivalent to light levels in clearings or forest edges. This effect limits the density of palm species with low tolerance to prolonged light exposure, which includes A. mexicanum (^{Martínez-Ramos et al. 2016}, ^{Hernández-Ruedas et al. 2018}). Reduced exposure to light would also explain the increase of individuals of this species in the LUs with less sloping relief (LF and PD).

As expected, slight variations in the tree community, together with the abiotic factors, influenced the palm community (Figure 3). This effect is significant only when soil-topography variables are included in the analysis (Table 3). Interspecific competition with the tree community could explain this relationship. The interaction between soil and tree community had a positive effect on palm density, especially in landscape unit whit high stoniness, where palm species were favored due to their more efficient use of sources, as soil water, compared to tree species. This result lends support to the interpretations of ^{Martínez-Ramos et al. (2016)} by indicating that a higher density of C. alternans, C. tepejilote, and A. mexicanum (Figure 3B) alters the composition and structure of tree community. However, our results indicate that the composition of the palm community is driven by environmental heterogeneity.

The increase in the populations of C. elatior and C. ernesti-augusti in SDS could be related to the decrease in the density of A. mexicanum and C. pinnatifrons. Particularly, it has been reported that A. mexicanum in the LTBS is a highly competitive species that affects density of another palm and tree species (^{Noriega‐Piña et al. 2021}). ^{Cepeda-Cornejo & Dirzo (2010)} suggested that there is a microhabitat separation, with C. pinnatifrons distributed at low elevation sites, while C. ernesti-augusti typically occurs at slightly higher elevations. This relationship suggests that the spatial variation of the most abundant palm species, in the context of environmental heterogeneity, affects the density of other palms. This result is highly relevant because palm community is an important structural component in the understory layer (plants up to 10 m) (^{Ibarra-Manríquez 1988}). However, further studies are required to evaluate the density-dependent effects on palm populations, considering the density of trees-relief-light interaction.

Our results show that there is high edaphic and topographic variation that influences the distribution of species at the landscape scale. However, geopedological heterogeneity has not been considered to explain the distribution patterns and species diversity in the TRF studied. To integrate the complex variation of both abiotic and biotic ecosystem components, we proved that landscape stratification into spatially explicit LUs is highly effective. We found that palm species density differs among LUs, and our study supports the idea that abiotic and biotic components simultaneously regulate the palm composition and abundance. However, the response to each factor differs among species. Therefore, further studies are needed to fill the information gap on these species’ functional attributes and ecological niches in the context of environmental heterogeneity. While the composition of adult palms was related mainly to abiotic factors at the landscape level, there are other factors such as biotic interactions (e.g., herbivory, pollination, seed dispersal, seed predation) that influence palm composition, and operate at different scales and stages of the life cycle of palms. The systematic integration of landscape elements allows us to address the multifactorial effect of the distribution and density of palm species. Therefore, under this approach, the analysis of other plant groups together (e.g., trees or lianas) can provide useful information to understand species coexistence in plant communities and critical factors for their conservation, as well as to understand the importance of environmental filters in the face of environmental changes and their effect on the distribution of species.

Supplementary material

Supplemental data for this article can be accessed here: https://doi.org/10.17129/botsci.3204

Supplemental material