Environmental and spatial processes shaping Quercus‐dominated forest communities in the Neotropics

Abstract

Quercus‐dominated forests in the Neotropics compose one of the broadest distributed ecosystems of mountainous zones whose species distribution has been explained by climate change adaptation over an abrupt physiography. This study aims to comprehend the contribution of environmental and spatial processes influential in species composition of Quercus‐dominated forests. A database of 86 plots, randomly located in a highly diverse spot of the Mexican Neotropics, was examined; it involved vegetation, environment, and spatial data. Local spatial variables were derived from distance‐based Moran's eigenvector maps, while regional spatial variables were obtained from geographical coordinates (latitude/longitude). Redundancy analysis and variance partitioning were run to reveal the relative importance of environmental and spatial processes explaining plant species composition of Quercus‐dominated forests and ascertain the environmental variables controlling the species composition at the different spatial scales (local and regional). Results showed that both local and regional spatial variables significantly control species composition in Quercus‐dominated forests, with the regional scale accounting for a higher explained variance than the local scale. At the regional scale, elevation, aspect, litter, maturity, and Mg were significant in explaining species community composition; elevation variance was twofold higher than variance for aspect and litter. At the local scale, elevation accounted for the highest variation followed by maturity, P, N, and Mg. The analysis also tested the variation of species community composition over environmental gradients mainly determined by elevation. This study indicated the importance of elevation at both local and regional spatial scales as an important driver in Quercus‐dominated forests distribution, probably related to Quercus species adaptation to temperature and precipitation. Finally, evidence from this research emphasizes the significance and necessity of a multiscale analysis to discern the spatial structured environment responsible for the high species diversity and broad distribution of forests dominated by Quercus genus in the Neotropics.

INTRODUCTION

Quercus (oaks) is a genus globally represented by 435 species, mostly distributed in North America with 240 species, followed by China with 100 species (Denk et al., 2017; Huang et al., 2013; Nixon, 1997). Species of this genus thrive well in various forest ecosystems, from temperate to tropical, from dry to wet, and from deciduous to evergreen forests (Bölöni et al., 2021; Manos & Hipp, 2021). Given their wide distribution, Quercus species are important components of high biodiversity regions, providing habitat for wildlife and ecosystem services (Gil-Pelegrín et al., 2017; Johnson et al., 2019; Mölder et al., 2019; Nixon, 2006).

The Mexican Neotropics is a region imprinted with one of the highest biodiversity worldwide (Figueroa-Rangel et al., 2020); their mountains are key areas for speciation due to their environmental heterogeneity and topographic complexity. This region harbors a high number of forest communities such as Quercus forests, Pinus–Quercus forests, and cloud forests, with a characteristic in common: the presence of Quercus species, either dominating in density or diameter and stem-high sizes (Alfaro-Reyna et al., 2020; Encina-Domínguez et al., 2018; Hipp et al., 2017).

Quercus diversification has been attributed to past climate change, ring porosity, seed dispersion, wind pollination, and deciduousness emergence (de la Riva et al., 2019; Hipp et al., 2017; Manos & Stanford, 2001).

In Mexico, the main environmental variables responsible for high Quercus diversity and broad distribution are climate (temperature and precipitation) and physiography (altitude and topography) (Morales Pacheco et al., 2018; de la Riva et al., 2019); the abrupt topography of its mountainous regions is characterized by a steep elevation gradient over short distances, changes in slopes percentages, and slope orientation producing a multiplicity of microclimate conditions where temperature and precipitation fluctuate at small spatial scales (Krasilnikov et al., 2013). At local scale, soil variables such as chemical composition, humidity, litter, and organic matter content (OMC) are associated with Quercus distribution (Aguilar-Romero et al., 2016); at regional scale, climate allied to physiography is important driver determining Quercus diversity (Martin et al., 2021; Uribe-Salas et al., 2018).

Spatially structured environmental variation plays an important role in plant distribution (Dray et al., 2012); consequently, the discrimination of spatial ecological structures in species distribution, composition, and diversity is a pivotal question in community ecology (Legendre, 1993). However, the spatial scale can also constitute a confounding variable; thus, it should be included as a predictor or a covariable in statistical models (Dray et al., 2006). In general, the local scale is considered the spatial autocorrelation created by community dynamics, while the regional scale corresponds to the scale of the environmental drivers (Dray et al., 2012).

To discern the multiscale framework determining Quercus diversity and its distribution, it is essential to ascertain the scales at which the environmental processes are shaping Quercus-dominated communities. In this context, the present study aimed to understand the contribution of environmental and spatial processes shaping species community composition in Quercus-dominated forests by addressing two questions: (1) What is the contribution of regional and local spatial processes elucidating variation in species composition of Quercus-dominated communities? (2) Which environmental variables are controlling species composition at the different spatial scales (local and regional)? We hypothesized that the environment at the local spatial scale is determining the highest variation in species composition.

METHODS Study area

The study was conducted in the Sierra de Manantlan in west-central Mexico (130,000 ha), between the coordinates: 19°24′32″–19°31′02″ N and 103°57′44″–104°01′09″ W (Figure 1). A total of 13 vegetation types according to Rzedowski (2006) have been described in the Reserve encompassing an elevation range from 430 to 2450 m asl. The main vegetation types at high and medium elevations consist of Quercus and Pinus–Quercus forests (Olvera-Vargas et al., 2010), cloud forest (Sánchez-Rodríguez et al., 2003; Santiago-Pérez et al., 2009), Pinus forest (Cerano-Paredes et al., 2015; Hernández Vargas et al., 2000), and Abies forests (Figueroa-Rangel & Moreno-Gomez, 1993); at low-elevations, tropical dry forests dominate the landscape (Olvera-Vargas, 2006). Temperature ranges from 12 to 18°C, mean annual precipitation around 1300 mm; precipitation peaks between July and September, while dry season occurs between October and June (Olvera-Vargas & Figueroa-Rangel, 2018). Physiography is one of the most outstanding features of the study area, portrayed by disjointed sandy ridges with a wide variety of slopes (10%–60%) and abrupt ravines of nearly 90%; there are also numerous sinkholes and limestone karst of variable depths (Olvera-Vargas & Figueroa-Rangel, 2012).

FIGURE 1. Map of the study area with the location of the 86 plots in Quercus-dominated forest communities in Sierra de Manantlan, Biosphere Reserve, west-central Mexico.

Data collection

The database consisted of 86 circular permanent plots, 500 m² each, established at random in Quercus-dominated forest at the eastern portion of the Sierra de Manantlan. In each plot, all individual trees with diameter at breast height (DBH) ≥5 cm and ≥1.30 m tall were tagged, counted, their DBH recorded, and identified to species level.

In each plot, the following environmental variables were recorded according to the protocol by Olvera-Vargas et al. (2014): elevation, aspect, slope, and catena in seven classes (flat terrain, ridge slope, upper slope, middle slope, lower slope, base slope, and gully).

Soil chemistry and texture edaphic information was derived from composite soil samples (c. 30 cm deep excluding the litter layer) collected from three points within each plot. The following variables were derived from these samples: pH, cation exchangeable capacity, OMC, and soil nutrient elements (Ca, K, Mg, Mn, N, and P).

Litter was measured in the field considering the depth of vegetable material on the forest floor measured in centimeters. The percentage of rocks and stones on the forest floor were visually assessed. Tree component included canopy openness and maturity; maturity was classified in four categories: immature, young reproductive, mature, and old-growth (Olvera-Vargas et al., 1996). Canopy openness entails the percentage of sky not covered by vegetation (Mitchell & Whitmore, 1993); its estimation was based on hemispherical photographs (8 in total), taken at every cardinal point around the center of the 500-m² plot. A Cannon 7.5/5.6-mm fish-eye lens mounted on a Cannon AE-1 Program body was used to take the photographs at 1.50 m above forest floor. The photos were further analyzed for canopy openness with HemiView Canopy Analysis Software v. 2.1. Geographical coordinates (latitude and longitude) were also documented for each plot.

Data analysis

Tree species importance was estimated by the Importance Value Index (IVI) (Etherington et al., 1987). This index is composed by the sum of relative density (number of individuals per species/total number of individuals all species), relative frequency (number of plots in which the species is present/total number of plots), and relative basal area (total basal area per species/total basal area all species).

Generalized dissimilarity models (GDM) were run to model species composition variation as a function of environment and geographical distance (latitude and longitude) using dissimilarity matrices (Bray–Curtis on species abundance data and Euclidean on environmental and geographical variables). This is a method extended from matrix regression analysis; it uses adaptable splines and generalized linear modeling. It is appropriated when two nonlinearity elements typical of ecological data sets are encountered: the curvilinear association between increasing ecological distances, and observed compositional dissimilarity between plots and variation in beta diversity along environmental gradients (Ferrier et al., 2007). The model produced splines, indicating how species composition change along geographical distance and with each environmental variable (Fitzpatrick et al., 2021).

To evaluate the influence of space and based on the sampling procedure, two classes of spatial variables were used to represent the spatial variation at regional and local scale. The geographical coordinates (X-longitude and Y-latitude) of each plot were used to represent the global trend at regional scale; these coordinates were used to identify the linear trend at a regional scale. Spatial variables at local scale were generated through distance-based Moran's eigenvector maps (dbMEM); these variables can be used in regression or canonical analysis as explanatory variables. The dbMEM analysis produces a matrix of pairwise geographical distance among plots; then, the matrix is truncated to preserve only those distances among close neighbor plots (Borcard et al., 2011; Dray et al., 2006), to keep all plots connected in a network. Subsequently, a principal coordinate analysis (PCoA) was ran on the truncated matrix and its eigenvalues with positive signs were considered as spatial variables; the sizes of the eigenvectors were used to represent spatial structures (Borcard et al., 2011). As a result, 14 eigenvectors (dbMEM) with positive values (positive spatial autocorrelation) were produced by the PCoA, which were used as spatial predictors in RDA. The dbMEM analysis was developed in R v4.05 package with function adespatial (R Core Team, 2021).

Redundancy analysis was performed to determine the variables explaining community species composition in Quercus-dominated forests. The response data consisted of a matrix with species abundance by plot, which was Hellinger transformed to reduce the weight of abundant species following Legendre and Gallagher (2001). As explanatory variables, we used environmental and spatial variables (X and Y coordinates and dbMEM).

Previous to RDA, the forward selection criterion through the adjusted R² (R²_adj) and the alpha significant level (95%) were applied to reduce the explanatory variables (dbMEM's and environmental variables), keeping those that better describe the most variation in the species matrix (Borcard et al., 2011; Dray et al., 2006).

For environmental variables, forward selection was implemented on the undetrended response variables. For dbMEM variables, it was run on the response variables detrended by the X and Y coordinates (once the trend was significant) (Borcard et al., 2011).

After environmental variables selection, collinearity was evaluated by the variance inflation factor (VIF), reserving those variables with VIF ≥ 10. Finally, the explanatory variables consisted of (1) the selected environmental variables, (2) the selected dbMEM variables representing the local scale, and (3) the X and Y coordinates for the regional scale. Variance partitioning and RDA were run on the three classes of the explanatory variables to calculate their contribution to variation in community species composition (Peres-Neto et al., 2006); their significance was test by permutation with 999 runs.

To determine how species composition varied along the environmental gradient, RDA was conducted only on environmental variables for community species composition and, the first two axes of RDA of Quercus-dominated forest community composition were plotted together with environmental variables.

To further evaluate the most important environmental variables at the two different spatial scales, we perform RDA of community species composition on spatial variables at regional and local scale; subsequently, for each of the two spatial scales, we run RDA using the fitted canonical values of the community species composition on environmental variables. Finally, ANOVAs were run to test the RDAs to inspect the significance of every environmental variable by the two different scales (regional and local).

All statistical analysis was performed in R v4.05 using the vegan, adespatial, and gdm packages (R Core Team, 2021).

RESULTS Community composition of Quercus-dominated forest

Quercus-dominated forest was composed of 38 species from 28 genera; the most important genera according to IVI was Quercus with 9 species. At species level, Quercus crassipes was the most important with the highest IVI (58.64%), followed by Ternstroemia lineata and Quercus laurina (42.73% and 26.21%, respectively); in terms of relative density, T. lineata was the most abundant (26.5%) followed by Q. crassipes (23.7%) and Q. laurina (8.8%), while the rest of the species contributed with low densities. The following five Quercus species were the most dominant using basal area as an indicator: Q. crassipes (26.4%), Q. laurina (11.5%), Quercus rugosa (8.9%), Quercus candicans (7.8%), and Quercus castanea (7.1%); the most frequent species were T. lineata (10%), Q. crassipes (8.3%), Styrax ramirezii (6.9%), and Alnus jorullensis (6.4%) (Figure 2a).

FIGURE 2. Importance Value Index (IVI), relative density, dominance, and frequency for all species recorded in the Quercus-dominated forest community (a). Boxplot for number of individuals (b) and mean diameter at breast height (DBH) (c) for Quercus species (acronyms for species in parenthesis).

The analysis of Quercus species revealed a high heterogeneity in number of individuals and mean DBH per plot (Figure 2b,c). Quercus crassipes was the most variable in terms of number of individuals with a median of 25 individuals by plot, followed by Q. laurina (median ~ 10 individuals). The highest DBH median was recorded for Q. candicans (~38 cm) although with a high variation in values; however, Q. rugosa was the species with the highest DBH variation even showing one individual with 121.3 cm.

The contribution of regional and local spatial processes shaping variation in species composition

The forward selection procedure retained the following environmental variables: elevation, slope, aspect, litter, OMC, pH, Ca, K, Mg, N, P, canopy openness, catena, and maturity.

For spatial variables, the selection criteria preserve six local spatial variables (dbMEM's): 3, 4, 5, 6, 7, and 8. Although complex pattern was revealed between variation in species composition with the environment and with space (Appendix S1: Figure S1), together, environment and spatial variables explained ~44% of total variation. After controlling for the other explanatory variables, partial RDA presented significant partial contribution of environmental and spatial variables (p ≤ 0.001) to the community species composition (Table 1). Variance partitioning results revealed that environmental and spatial variables contributed with ~20% and ~45%, respectively. Spatial contribution regarding geographical coordinates explained 26.75%, while dbMEM account for 18.17%. This result indicates that both local and regional spatial variables are robustly regulating community composition in Quercus-dominated forests (Table 1 and Figure 3).

TABLE 1 Redundancy analysis and variance partitioning of community species composition based on environmental explanatory variables (Envi) and two sets of spatial variables: (1) geographic coordinates (XY) and (2) dbMEM variables.

Variables	R2adj	df	F	p
Total	0.44287	22	17.488	0.001***
XY	0.2675	2	16.526	0.001***
XY|Envi + dbMEM	0.0640	2	4.7336	0.001***
dbMEM	0.1817	6	4.1459	0.001***
dbMEM|Envi + XY	0.1229	6	3.5383	0.001***
Envi	0.2000	14	2.5185	0.001***
Envi|XY + dbMEM	0.0554	14	1.5477	0.007**

Abbreviation: dbMEM, distance-based Moran's eigenvector maps. *equals significance at 0.05; **at 0.001; ***at 0.001.

FIGURE 3. Variance partitioning for species community composition in Quercus-dominated forest in west-central Mexico, constructed on environmental variables, regional (geographical coordinates: X-longitude and Y-latitude), and local distance-based Moran's eigenvector maps (dbMEM) spatial variables.

Environment alone (not considering spatial variables) accounted for ~33% (proportion = 0.3318). The most important and significant environmental variables were elevation, litter, aspect, and maturity (Table 2).

TABLE 2 Environmental variables contribution to explain species community composition by regional and local spatial scale in Quercus-dominated forest.

Scale	Variables	df	Variance	F	Pr(<F)
Regional	Elevation	1	0.0035884	32.488	0.001***
	Aspect	1	0.0019065	17.2606	0.001***
	Litter	1	0.0018502	16.7513	0.001***
	Mg	1	0.0013144	11.9005	0.001***
	Maturity	1	0.0015845	14.3456	0.002**
Local	Elevation	1	0.0019655	19.2735	0.001***
	Maturity	1	0.0008483	8.3184	0.007**
	P	1	0.0006903	6.7688	0.013*
	Mg	1	0.0004269	4.1858	0.041*
	N	1	0.0004467	4.3801	0.042*
No space	Elevation	1	0.06654	10.2810	0.001***
	Litter	1	0.04851	7.4955	0.001***
	Aspect	1	0.02387	3.6884	0.007**
	Maturity	1	0.01701	2.6281	0.029*

*equals significance at 0.05; **at 0.001; ***at 0.001.

Environment locally structured accounted for ~42% (proportion = 0.4198), while environment regionally structured accounted for ~59% (proportion = 0.58915).

At the regional scale, two physiographic variables (elevation and aspect), two variables related to soil (litter and Mg), and one variable related to vegetation (maturity) were significant in the explanation of species community composition; elevation variance was almost two times higher than aspect and litter, signifying its highest explanation on community composition (Table 2). At local scale, elevation accounted for the highest variance followed by maturity and soil variables P, N, and Mg (Table 2).

As maturity, the only variable related to vegetation, was significant at both spatial scales, a multivariate homogeneity of group dispersions through a permutation of beta diversity analysis was run to further understand how species composition varies across plots with different maturity (Anderson, 2006; Anderson et al., 2006). The beta diversity permutation demonstrated species composition varied across plots with different maturity (F_3,82 = 9.54, p = 0.000017).

Results of the RDA using the fitted canonical values of the species community composition on environmental variables at local scale were only significant at axis 1 (p = 0.002). The pattern observed revealed that positive scores were located at the center of the study area, with the highest sizes between 19.44 and 19.46 in latitude, corresponding to plots 1–5, 40–42, and 36–38, which are located in intermediate elevations and mostly dominated by Q. castanea, Q. candicans, and Q. rugosa (Figures 4 and 5). Negative scores were located between 19.43 and 19.44 in latitude, with the highest sizes for plot 39 and plots 46–48 found in the lower elevations of the study plots, with Q. castanea dominating the canopy; an exception was plot 28 with negative scores and big size but located in the highest elevation with Q. crassipes dominating the canopy, low values in Mg and classified in a young reproductive stage (Figures 4 and 5). Positive scores with the lowest sizes were between 19.47 and 19.48 in latitude; these plots (52–56 and 79–81) were located in high-elevations, mostly at mature stages and high values in Mg, but low in P and intermediate and low values of N; species composition corresponded to a combination of Q. laurina, Quercus gentryii, and Quercus scytophylla. Negative scores with the lowest sizes corresponded to plots 6–11 and 26 in high-elevations (over 2300 m asl) with heterogeneous (high, intermediate, and low) values in Mg, N, and P. Finally, negative scores with intermediate sizes were evident for plots 73–78, 82, and 83 at the top left of the figure between 19.48 and 19.49 in latitude; these plots occurred in elevations from 2250 to 2300 m asl, classified as old-growth plots with high values in Mg but low in P and N; they were mainly dominated by Q. laurina and Q. aff. excelsa (Figures 4 and 5).

FIGURE 4. Plot location (longitude and latitude) produced by the distance-based Moran's eigenvector maps (dbMEM) of the redundancy analysis fitted plot scores for axis 1; color and sizes of the squares represent the value of the eigenvalues modeling spatial correlation; it also denotes the spatial scale regarding Quercus-dominated community composition; black squares are positive, and white squares are negative with sizes proportional to the absolute values; high sizes are associated with the regional scale, while small sizes are related to local scale. Distance resolution = 0.01° latitude and longitude (d = 0.01).

FIGURE 5. Environmental variables significant in explaining Quercus-dominated forest communities at local and regional spatial scales grouped by Quercus species with the plot number where they were recorded: in parts per million (ppm) and in meters above sea level (m asl).

Environmental variables alone explained ~33% of the Quercus-dominated forest community composition according to the RDA results (Figure 6). Quercus species were clearly separated in the ordination space with Q. crassipes on the left side of the ordination diagram; Q. laurina to the right top of the ordination diagram, together with T. lineata and S. ramirezii; Q. castanea, and Q. candicans positioned at the bottom on the right side, opposite to the elevation vector; Q. rugosa was located at the center together with Q. scytophylla, Q. gentryii, Quercus obtusata, and Q. aff. excelsa (Figure 6).

FIGURE 6. The first two axes of redundancy analysis (RDA) of Quercus-dominated forest community composition on environmental variables. Red rectangle on the right contains the species acronyms for those species clustered at the center of the diagram (species names for acronyms as in Figure 2). Ca, calcium; CAOP, canopy openness; K, potassium; Mg, magnesium; N, nitrogen; OMC, organic matter content; P, phosphorus

Results of GDM (percent deviance explained = 14.358) revealed geographical distance, using latitude and longitude, was important explaining changes in species composition in oak-dominated forests. For the environment, the most important variables were elevation, maturity, litter, phosphorus, and nitrogen. The highest change in species composition occurred from 0.04 to 0.06° (latitude/longitude); it increased at elevations above 2200 m asl, above 6 cm in litter and in old-growth plots regarding maturity class; and change in species composition stabilized at 15 ppm in phosphorus and nitrogen (Appendix S1: Figure S2).

DISCUSSION

Species composition of the Quercus-dominated forest in the Sierra de Manantlan Biosphere Reserve in the Neotropics revealed a diverse community where Quercus species are co-occurring at different mixtures and levels of abundance and sizes. Considering the total number of species in the study area, Quercus species represented around 24% of the entire community. This is a common feature of mountain forests in Mexico, where Quercus species dominate numerous vegetation types and environments at elevations from 50 to 3000 m (Arenas-Navarro et al., 2020; Luna-Vega et al., 2004).

The nine Quercus species reported in this study were also very different in terms of structural variables and, although only four Quercus species depicted an IVI greater than 15%, the most important species in the community was Q. crassipes with high IVI (58.64%); this Quercus species surpassed T. lineata subsp. lineata (IVI = 42.73%), the most abundant and frequent species but in the second place in terms of IVI due to its individuals were represented by smaller DBH sizes and consequently lower basal areas. Quercus laurina, Q. castanea, Q. rugosa, and Q. candicans followed in order of importance; this result was mainly due to their higher basal areas. Species with intermediate IVI were Q. aff. excelsa and Q. scytophylla, while species with low IVI were Q. obtusata and Q. gentryii. All Quercus species belong to the Mexico province, which include Mexico, Central America, Arizona, and New Mexico (Edler et al., 2017); considering their taxonomy, seven Quercus species belong to section Lobatae (red Quercus): Q. candicans, Q. castanea, Q. crassipes, Q. gentryii, Q. laurina, Q. scytophylla, and Q. aff. excelsa, while only two (Q. obtusata and Q. rugosa) are in section Quercus (white Quercus) (Hipp et al., 2017; Uribe-Salas et al., 2018; Valencia-A & Gual-Díaz, 2014). In leaf habit, only Q. obtusata is deciduous, Q. laurina and Q. castanea are brevideciduous, and the rest of the species are evergreen; in terms of evolutionary history, a closer relationship is reported among Q. crassipes, Q. gentryii, Q. candicans, Q. castanea, and Q. scytophylla (Oyama et al., 2018; Vega-Ortega, 2021).

Ternstroemia lineata, the second important species in the study area, has been reported in cloudy habitats, with high humidity, over elevations from 1400 to 3140 m (Alcantara-Ayala et al., 2020). The most common vegetation types where T. lineata develops in Mexico are the cloud forest, Quercus forest, and Quercus–Pinus forest (Benz et al., 1996, Luna-Vega et al., 2004); this is an species considered as near threatened according to the Red List of Mexican cloud forest (González-Espinosa et al., 2011); it also represents a target species in studies of niche modeling used to predict species distribution, in conservation evaluations, and in responses to climate change (Alcantara-Ayala et al., 2020; Luna-Vega et al., 2012). In Quercus-dominated forest, it is one of the most abundant successful species regenerating in the understory (Olvera-Vargas & Figueroa-Rangel, 2018). Studies on population ecology revealed T. lineata as a very dynamic species with moderate levels of recruitment (1.4%–2.25%) and mortality (0.54%–1.15%) rates (Olvera-Vargas et al., 2015).

Quercus-dominated forest in the study area is shaped significantly by the contributions of both environmental and spatial variables. Both sets of variables contribute almost equally to explain variation in species composition in Quercus-dominated forests, implying that both neutral (sensu Hubbell, 2001) and niche-based processes are essential for these communities (Gravel et al., 2006). The regional scale, determined by the geographical coordinates, exceeded in around 8.5% to the local scale with elevation as the main constraint at both regional and local scales; therefore, the hypothesis, posed in this study, that the environment at the local spatial scale was determining the highest variation in species composition, was rejected. Elevation has been one of the most quoted environmental variables explaining community species composition, specifically in Quercus, Pinus, and Pinus–Quercus forests of the Mexican Neotropics (Encina-Domínguez et al., 2018; Sandoval-García et al., 2020). Mainly for Quercus forests, a number of authors found an important relationship between elevation and Quercus diversity, mostly driven by an adaptation to temperature and precipitation (Sabás-Rosales et al., 2015; Vega-Ortega, 2021). Aguilar-Romero et al. (2016) reported Q. candicans, Q. crassipes, and Q. rugosa frequently located in environments with high rainfall and low temperature. In this study, Q. crassipes, Q. laurina, Q. gentryii, and Q. scytophylla were located over the highest elevations associated to low temperatures and high rainfall.

At the regional scale, aspect was also important; this environmental variable is an influential factor for species distribution (Kabrick et al., 2004) and, together with slope position and soil depth, they regulate soil water availability, which in turn reallocated Mg available in soils along time (Kabrick et al., 2011). Magnesium was significant at both scales; this is one of the most important soil nutrients constituting up 10%–30% of the exchangeable cations; it is considered as a fundamental plant nutrient, contributing in cell processes such as protein synthesis and pH regulation; consequently, its scarcity conduces to a decrease in plant growth (Oloffson, 2016). Mean exchangeable Mg concentrations in forest soils have been reported to range from 0.1 to 0.7 g/kg, values, which are within the values, recorded in the present study (0.2–0.6). Quercus species were more abundant in low Mg concentrations in the Ozark Highland forest in United States (Kabrick et al., 2011), similarly to the findings in the present study, where low Mg values were recorded in plots with the highest number of Quercus species.

Nitrogen and phosphorous were both significant at the local scale; these two elements are fundamental macronutrients regulating plant growth and productivity in terrestrial ecosystems with a crucial function in plant metabolism (Crous et al., 2019); while N is associated with plant carbon assimilation, P takes part in energy metabolism forming macromolecules (e.g., nucleic acids and ribosomes) (Rao & Terry, 1995; Reich et al., 2008). Nitrogen in ecosystems derives from the biological fixation of N₂ and from the atmosphere; it stores along the development of the ecosystems, strongly tied to the decomposition of soil organic matter, predominantly in mature plots (Brödlin et al., 2019). Phosphorous originates from rock weathering and, as pedogenesis evolves, its content decreases giving place to inorganic P, which is tied to secondary minerals (some derived from Mg such as saponites) (Turner et al., 2007; Walker & Syers, 1976); therefore, in mature plots such as the Quercus-dominated forests in our study area, plant nutrition depends on the mineralization of N and P (Bünemann et al., 2016; Lang et al., 2017). Low quantities of N have been linked to increased diversity in plant communities; the accumulation of this element has even considered as central factor driving changes in species composition along a gradient of different ecosystems by means of competitive interactions; this process has also been related to the amount of P in soil (Bobbink et al., 1998). In this study, low values in N were linked to plots with species composition resembling cloud forests whose biodiversity is one of the highest among the mountain forest in Mexico (Gual-Díaz & Rendón-Correa, 2014). Moreover, the elevation gradient in this study is related to low temperatures, which is an important factor in the regulation of N and P (Zheng et al., 2017); steep slopes also allow deep soil accumulation, which increases bedrock weathering to preserve mineral-derived nutrients, such as Ca, Mg, K, and P. High values of precipitation and humidity are also critical factors affecting P distribution in soil, even more limiting that N (which is deposited largely from the atmosphere); therefore, its mineralization in organic form accounts for its supply to plants (Shiau et al., 2018). The process explained above must be taking place in the study area as litter was a significant variable explaining species composition in the Quercus-dominated forest.

Maturity was also a significant variable explaining species composition in Quercus-dominated forests at both local and regional scales; plots dominated by Q. crassipes were mostly classified in the young reproductive stage and few plots in the old-growth stage; on the contrary, almost the totality of the plots dominated by Q. laurina were classified in the old-growth stage. Quercus crassipes plots mostly correspond to dry habitats with open environments, while Q. laurina is mostly growing in the wettest habitats of the study area (Olvera-Vargas, 2006), particularly in plots dominated by cloud forest, which is mainly composed by mature and old-growth plots. In the opposite, Q. crassipes is mostly associated with Pinus leiophylla and A. jorullensis; therefore, it is likely that these plots were logged in the past to extract Pinus species; such activity could have opened the canopy allowing forest regeneration, which may explain the presence of a young reproductive stage.

Finally, plots dominated by Q. castanea, Q. candicans, and Q. rugosa were mostly linked to regional scales, while Q. laurina, Q. gentryii, Q. aff. excelsa, and Q. scytophylla were related to local scales with elevation, maturity, and soil-related variables (P, N, and Mg) significantly explaining species composition of Quercus-dominated forests, mainly those located in cloud forests plots; therefore, these results are consistent with the location of cloud forests in restricted and fragmented habitats in many regions of Mexico where Quercus species play a key ecological role (Valencia-A & Gual-Díaz, 2014).

CONCLUSIONS

The environment regionally structured accounted for the highest explained variation in comparison with the local scale with elevation as the main constraint at both scales; Q. crassipes, Q. laurina, Q. gentryii, and Q. scytophylla were the species located at the highest elevations. Magnesium was the only soil element significant at both scales; nitrogen and phosphorous were only significant at the local scale. Maturity was a significant variable at both spatial scales discerning differences among plots dominated by Q. laurina (mostly classified in the old-growth stage) and plots dominated by Q. crassipes (mostly classified in the young reproductive stage).

Finally, the discernment of the local and regional scales structuring the environment ultimately defined variation in species composition in Quercus-dominated forest; the evidence presented in this study provides an insight into spatial local ecological processes that have vastly molded the high plant species diversity of the Mexican Neotropics.

ACKNOWLEDGMENTS

We are greatly indebted to J. Martín Vázquez-López, Oscar Sánchez Jiménez, and Abel Ceja Gutiérrez for his appreciated assistance during the fieldwork. Universidad de Guadalajara provided funds for the fieldwork. We thank two anonymous reviewers for their thoughtful feedback.

CONFLICT OF INTEREST

The authors declare no conflict of interest.

DATA AVAILABILITY STATEMENT

Data (Olvera Vargas & Figueroa-Rangel, 2022) are available from Mendeley: https://data.mendeley.com/datasets/szdzkbzzxd/1.