Diversity and niche differentiation of a mixed pine–oak forest in the Sierra Norte, Oaxaca, Mexico

Abstract

This study examined tree species diversity, distribution, and community differentiation patterns along an elevation gradient in pine–oak forest in the Sierra Norte mountains in Oaxaca, Mexico. Pine and oak are ecologically and economically valuable both locally and globally, but their dynamics are poorly understood in seasonally dry montane forests. This is a biome that is both widespread, with high human use and importance, and widely understudied. The community‐managed forest we studied contained high levels of tree diversity (32 total species), with especially high levels of oak (10 species) and pine (eight species) diversity compared to other pine–oak forests in Mexico. Tree communities in the study area demonstrated high levels of species turnover across sites, especially at mid and low elevations, as well as high levels of oak species coexistence within communities, with a mean of three oak species per 1000 m². We identified three distinct tree vegetation types using multivariate ordination and cluster analyses and found that both tree distributions and community assemblages are primarily differentiated by elevation, but also by soil type, topography, and likely successional disturbance from historical land use. Oak communities in the study area followed patterns of phylogenetic overdispersion with species from different sections (red and white) co‐occurring more frequently than species from the same section, and demonstrated differences in reported reproductive phenology, with coexisting species alternately fruiting in rainy and dry seasons. This differentiation in both oak species' environmental associations and in fruiting phenology has important management and conservation implications as Mexico becomes more arid with climate change. This study also provided key information for local management as different forest types should have different silvicultural management regimes, as well as essential baseline data useful for a broader theoretical understanding of how closely related species coexist in communities.

Introduction

Pine and oak species are extremely important both ecologically and economically around the world. They are harvested for timber, fuelwood, resins, and fodder and provide essential food sources, habitat for wildlife, and key ecosystem services associated with forest cover in a wide range of environmental conditions. In the continental United States and Mexico, pine (Pinus) and oak (Quercus) are the top two genera in terms of both species diversity and biomass in natural forests, and also in timber production, where oak and pine comprise between 20–29% and 18–25%, respectively, of the total aboveground biomass (Cavender‐Bares 2016). Given the importance of these floristic groups, increasing understanding of pine–oak ecological dynamics is key for improved conservation and silvicultural management of these forest types.

However, the mechanisms driving pine–oak forest composition and ecosystem functioning remain poorly understood, especially in seasonally dry and subtropical systems in Mesoamerica and Asia where climate change and anthropogenic disturbance have strong impacts. The lack of research is especially striking in Mexico, which is the center of diversity for both the Quercus and Pinus genera. Mexico contains an estimated 49 of the 120 pine species worldwide, 22 of which are endemic to the country (Farjon and Styles 1997, Gernandt and Pérez de la Rosa 2014), and 172 of the world's 435 oak species, 85 of which are endemic to the country (Zavala‐Chávez 1998, Nixon 2006, Hipp et al. 2018). In southern Mexico, the state of Oaxaca is home to some of the highest diversity and endemism of conifers in the country (Contreras‐Medina and Luna‐Vega 2007), and the highest richness of oak in the country, with 48 Quercus species found in the state alone (Valencia 2004). Mexican pine–oak forests are also unique systems in that they contain a high diversity of tree species but a relatively low diversity of tree genera, making them particularly interesting because the overarching processes governing community structure in forests dominated by closely related species remain poorly understood.

Montane pine–oak forests are some of the most vulnerable ecosystems to the impacts of climate change at both the ecosystem (Villers‐Ruiz and Trejo‐Vázquez 1997) and species levels in Mexico (Gómez‐Mendoza and Arriaga 2007, Galicia et al. 2015). This is in part due to the likely increase in intensity of seasonal drought as the country becomes more arid (Seager et al. 2009, IPCC 2013), as well as the higher elevation at which these forests are found (Lenoir et al. 2008, Anderegg and HilleRisLambers 2019). These shifts in drought may be especially impactful in oak forests, as trade‐offs between drought tolerance and growth rates are a mechanism that drives oak species coexistence in seasonally dry forests (Poulos et al. 2007, Aguilar‐Romero et al. 2017, Fallon and Cavender‐Bares 2018). Unfortunately, Mexican montane pine–oak forests have received little attention in research or reviews of forest resources, forest ecology, or global climate change models (Gómez‐Mendoza and Arriaga 2007, Richardson et al. 2007, Galicia et al. 2015).

For forest management and conservation efforts to be effective in these pine–oak forests, there is a need for increased understanding of the floristic and species composition dynamics across altitudinal and environmental gradients. Studies of ecological niche differentiation sensu Grinnell (1917) examine species' different requirements for limiting resources and the resulting differences in associations between species and environmental variables, life history, and growth strategies. Such studies are especially important in the face of global climate change, as such knowledge can help predict how stand composition and dynamics may change as precipitation and temperatures shift, and can provide information about where and how to target management efforts or conservation actions to aid the establishment and growth of particular species and their associations (Guisan et al. 2013). The few studies that have examined niche differentiation in mixed pine–oak forests have relied on species presence data from herbarium records of oak and pine species on a broad national or regional scale; examinations of actual species distributions within a forest are rare (Alba‐López et al. 2003, Gómez‐Mendoza and Arriaga 2007, Aguirre‐Gutiérrez et al. 2014).

This study examined pine and oak species distributions in a diverse community‐managed forest in the Sierra Norte region of Oaxaca, Mexico, in order to assess patterns of niche differentiation and species coexistence among these two genera. The few other studies of pine–oak forest in Oaxaca have examined tree community compositions, but have delineated vegetation types solely based on the dominant pine species in a stand without examining the broader floristic assemblages or examining the environmental factors driving community differentiation across the landscape (Castellanos‐Bolaños et al. 2010, Ríos‐Altamirano et al. 2016). While such an emphasis on pine is certainly understandable given its high economic importance, research that examines entire tree communities is key for more holistic forest management that includes the management of subsistence and smaller‐scale products such as oak charcoal and firewood, and is critical for providing baseline information for conservation efforts in the face of climate change. As such, this is the first study in this region to examine entire floristic groups for this forest type.

We addressed three main objectives in this study: (1) to characterize the diversity and floristic composition in this forest type; (2) to characterize patterns of niche differentiation among pine–oak communities, and examine the factors that influence where these assemblages are found across the landscape; and (3) to characterize the niche differentiation and distribution of individual pine and oak species.

Finally, although we focused on a particular forest in Oaxaca in this study, the patterns we describe in pine and oak species dynamics provide valuable insight to similar understudied forest types especially in Central and subtropical Asia (i.e., the Himalayan region) and the central mountains of northern Central America, where forests are also heavily relied on by local communities for pine timber and resin and oak firewood, leaf fodder, and charcoal.

Methods

Study location

We conducted this research in the communal forest of Nuevo Zoquiapam in Ixtlán Province in the Sierra Norte of Oaxaca, Mexico (located at 17°17′19″ N, 96°37′15″ W; Fig. 1). The Sierra Norte region in Oaxaca is one of the most diverse mountain systems in the country and has been considered a priority region for conservation in Mexico (Arriaga et al. 2000). Nuevo Zoquiapam is a predominantly Zapotec community composed of two small towns—San Matias Zoquiapam and Nuevo Zoquiapam—consisting of roughly 320 and 1300 inhabitants, respectively (INAFED 2010). Currently, the community holds 9554 ha of land, of which about 7340 ha are in forest. Within the forested area, 3620 ha are managed for timber production, 1085 ha are managed for pine resin production, 1380 ha are designated as areas for restoration due to low timber production, and 1250 ha are designated as protected areas. The forest ranges in elevation from 1800 m.a.s.l. to over 3200 m.a.s.l. and transitions from pure pine stands at highest elevations to mixtures of pine and oak to pure oak stands at the lowest elevations.

1 Fig.. Map of study site with inset showing location relative to the state of Oaxaca, and the location of Oaxaca within Mexico. The three shades of gray indicate the three elevation zones used in the stratified sampling design, and the stars show the locations of the two villages. The black dots indicate the 30 randomly located sampling sites; three plots are systematically located at each site, for a total of 90 plots.

Mean annual rainfall in the region ranges from 800 to 1200 mm and falls predominantly in May–October. Temperatures during the coldest month (usually December or January) range from −3°C to 18°C and during the warmest month (usually May) range from 6.5°C to 30°C. Soils in the region are derived from basalt material, have relatively high clay components, and are classified in the FAO system as Humic Acrisol (e.g., Humults in the Ultisol soil order; USDA Soil Conservation Service 1975, Saynes et al. 2012).

Field methods

We sampled tree species composition and stand structure across an elevation gradient divided into three zones of low (2300–2600 m asl), middle (2600–2900 m.a.s.l.), and high elevation (2900–3200 m.a.s.l.). In each zone, we randomly generated 10 site locations using ArcGIS 10.3, and at each of the 10 sites, we located three plots for a total of 30 plots per elevation zone (i.e., 90 plots across the entire forest; Environmental Systems Research Institute [ESRI] 2014). We used the random GPS point as the center point for each site and set each of the three plots 50 m from the center. The direction for the first plot was always set in the direction heading up the slope, and the remaining two were each set using a compass to be 120° apart from the first. If a plot fell on a trail, the distance between the center of the site and the plot location was extended in the same compass direction by an additional 20 m. Each plot was a 20 × 20 m square (400 m²) for a total of 1200 m² per site and 1.2 ha per elevation zone. All plots were adjusted to correct for slope, as measured by a clinometer and assigned to a 10° slope class. We generated alternate site points as well, and in two cases these were used when the initial point fell in small P. patula plantations.

We identified the species and measured the diameter for all trees with dbh >10 cm within each plot. We initially used morphotypes to identify tree species using the common names both in Spanish and in Zapotec with the help of a local field assistant. We later identified all morphotypes to species with the assistance of the herbarium at the Instituto Tecnológico del Valle de Oaxaca (ITVO) and stored voucher specimens there.

At each plot, we recorded altitude (m), aspect, and slope in degrees. In the first of the three plots at each site, we also took soil samples immediately below the litter and O horizon from the first 10 cm of mineral soil and measured depth of the organic layer. Mineral soil samples were analyzed for organic material content, pH, carbon, and nitrogen contents by the Laboratorio de Analysis de Suelo FYPA using Micro‐Kjeldahl protocols for nitrogen and Walkley and Black protocols for organic matter analyses (Walkley and Black 1934, Ma and Zuazaga 1942).

Data analysis

We calculated summary diversity statistics, including total species richness (i.e., across the entire sampling area), mean species richness per plot, and beta‐diversity metrics (Magurran 2004). We partitioned beta‐diversity into separate components of species turnover (measured as Simpson's dissimilarity) and species nestedness (measured as the nestedness‐resultant fraction of Sorenson's dissimilarity) using the betapart package in R (Baselga 2009, 2012, Baselga and Orme 2012). Species turnover represents the addition of new species across sites and can overemphasize the role of rare species, while nestedness looks at species loss (i.e., how often species‐poor sites are a subset of assemblages in more species‐rich areas; Wright and Reeves 1992). We also generated species accumulation curves and species rank curves using the vegan package for biodiversity (Oksanen et al. 2019). We performed these analyses on all samples together and on each of the three elevation zones separately. For rank abundance curves based on species' importance values, we calculated importance value as the sum of the relative frequency, relative density, and relative basal area of each species in each elevation zone. We performed all statistical analyses using the program R (R Core Team 2020).

To characterize the number of distinct vegetation communities present in the forest, we used two hierarchical cluster analyses with Ward's distance using the hclust function in the cluster package and the pvclust function in the pvclust package, which assigns P‐values to cluster partitions (Suzuki and Shimodaira 2015, Maechler et al. 2018). In the first cluster analysis, we clustered species based on square‐root‐transformed tree densities in all plots. In the second analysis, we clustered plots based on species composition and used these plot clusters as groups in the ordination analysis described below. For both cluster analyses, we used the NbClust function to determine the best number of clusters (Charrad et al. 2014). We also used the function multipatt from the indicspecies package in R to assess which species could be considered significant indicators for the clusters generated from the hierarchical clustering of plots (De Cáceres and Legendre 2009).

To visually explore community structure and niche differentiation with environmental factors, we used non‐metric multidimensional scaling analysis (NMDS) with a Bray–Curtis distance matrix on count data for all identified species by plot using the vegan package in R (Oksanen et al. 2019). We found a three‐dimensional solution was the better fit based on stress values and ran up to 20 Monte Carlo permutation tests, stopping once a convergent solution was reached. We then fit linear vectors of the environmental variables with the envfit function from the vegan package in R to examine which factors were significantly correlated with the spread of sites and species (Oksanen et al. 2019). We color‐coded the ordination diagram using the cluster solutions from the hierarchical cluster analysis for visual assessment of the cohesiveness of the clusters.

To more directly examine the relationship between environmental factors and species composition, we performed multivariate regression analysis using the environmental variables of altitude, aspect, slope, and soil characteristics (pH, depth of organic layer, percent organic matter, percent N, and percent C). For this analysis, we transformed the circular variable of aspect (in radians) into the continuous variables of northness and eastness using sine and cosine (Beers et al. 1966). For the response variables in the multivariate regression model, we used square‐root‐transformed tree counts from plots for the top 16 tree species (i.e., all species that were found in 10 or more plots) with the MANOVA function in R. We examined the direction of univariate correlations using linear regression coefficients on untransformed tree counts.

Finally, we examined the relationship in distribution and co‐occurrence between the two oak clades (red, Lobatae; and white, Quercus) and between individual species. To assess the relationship in distributions of red and white oaks, we used a log‐likelihood (G test) of distributions with the GTest function in the DescTools package (Signorell et al. 2020). The G test is a goodness‐of‐fit test of frequency distributions and has been recommended over the chi‐square test for this type of data on theoretical grounds (Hope 1968, Sokal and Rohlf 2012). To assess the co‐occurrence rates and patterns between species, we created a species co‐occurrence matrix based on the presence–absence data of all species found in at least five plots using the co‐occur package in R (Griffith et al. 2016).

Results

Tree species richness and diversity

We measured 2146 trees which were represented by 32 different tree species, 13 genera, and 12 families (we were unable to identify four understory tree morphospecies; see Appendix S1: Table S1 for a full list). We found 10 Quercus species (six in the red oak section Lobatae and four in the white oak section Quercus) and eight Pinus species (seven subgenus diploxylon and one subgenus haploxylon). Of the oak species, seven are endemic to Mexico (with one endemic to this region of southwestern Mexico) and the remaining three are endemic to Mexico and Central America. Of the pine species, five are endemic to Mexico, two are found in Mexico and Central America, and one is found in Mexico and parts of the southwestern United States. While species accumulation curves for the entire study area and for the high and the middle elevation zones indicated that sampling likely captured most species present, the curve for the low elevation appeared to be still increasing, implying that increased sampling might have uncovered more tree species in this zone (Fig. 2).

2 Fig.. Species accumulation curves by elevation zone (with gray open shapes) and for the entire study area (in filled black circles). Points show expected (mean) species richness with vertical lines representing 95% confidence intervals based on 100 permutations.

Mean species richness and Shannon's diversity index were lowest in the high elevation zone and highest in the middle elevation zone, although the highest total number of species was found in the low elevation (Table 1). The maximum richness per plot (400 m²) was 10 species (found in one plot in the middle elevation zone and one plot in the low elevation zone), and the minimum richness was one species (found at two plots in the high elevation zone and one plot in the low elevation zone). Mean richness per plot ranged from 3.4 to 5.6 species across the three zones (Table 1). Mean richness per plot of pine was 1.7 species and per site was 2.7 species with a maximum of four species per plot and five species per site. Mean richness per plot of oak was 2.3 species and per site was 3.4 species, with a maximum of six species per plot and seven species per site.

1 TableSummary of tree species diversity statistics by elevation zone (alpha diversity) and for the total forest (gamma diversity).

Zone	N	Species richness	Mean richness	Shannon's index	Species turnover	Species nestedness
High	30	11	3.4	1.83	0.830	0.066
Middle	30	22	6.5	2.47	0.841	0.047
Low	30	24	5.2	2.44	0.871	0.043
Total	90	32	5.0	2.85	0.958	0.014

Notes:

Species richness is the total number of tree species encountered in each zone (or in total), while mean richness is the mean number of species per plot (400 m²). Species turnover and species nestedness are two components of beta‐diversity. Species turnover is measured as Simpson's dissimilarity, and species nestedness is measured as the nestedness‐resultant fraction of Sorensen's dissimilarity.

The two metrics of beta‐diversity showed increase beta‐diversity at lower elevation zones (Table 1). Species turnover increased with decreasing elevation across the three zones, while species nestedness (which is negatively correlated with beta‐diversity) decreased with elevation such that high elevation zone had substantially higher levels of nestedness across plots compared to both the middle and low elevation zones.

For all zones, the dominant tree genera were pines and oaks, although the dominant oak and pine species shifted between elevation zones. Across all zones, the oak species were present at higher densities but pine trees represented larger proportions of basal area. Based on importance values, the most dominant species in the high and middle elevation zones were pines (P. hartwegii in the high and P. pseudostrobus in the middle), followed by oak species (Fig. 3). In the low elevation zone, the dominant species was an oak (Q. scytophylla) followed by several pine species (P. oaxacana and P. douglasiana; Fig. 3). In the high zone, the top five species represented 84% of all stems, and the top 10 species represented 99.8%. In the middle zone, the top five species represented 70.2% of all stems, and the top 10 species represented 85.6%. In the low zone, the top five species represented 66.2% of all stems, and the top 10 species represented 89.4% of all stems. Tree species showed distinct distributions along the elevation gradient, with some species more widely distributed than others (Fig. 4). For a summary of the soil characteristics by elevation zone, see Appendix S1: Table S2.

3 Fig.. Species rank curves based on species' importance value for each elevation zone. Importance values represent a combination of relative frequency, relative density, and relative basal area for each species in each elevation zone.

4 Fig.. Altitudinal distribution of the most common tree species shown as density per hectare by elevation. Oaks are highlighted in light gray and pines in dark gray.

Differentiation of tree communities

Hierarchical cluster analysis of species found three distinct clusters based on NbClust and P value determinations (P < 0.06, NbClust best cluster solution = 3; Fig. 5). The first cluster contained species predominantly found at high elevations (P. hartwegii, P. pseudostrobus, Q. crassifolia, Q. rugosa, Q. obtusata, and Q. laurina). The second cluster consisted of predominantly hardwood species found at middle and lower elevations (along with P. leiophylla and P. ayacahuite, which were both infrequently found in sampling). The third cluster contained species found predominantly at low elevations (P. oaxacana, P. douglasiana, P. patula, P. teocote, Q. conzattii, Q. scytophylla, and Arbutus xalapensis).

5 Fig.. Hierarchical cluster analysis of species using an Euclidean distance matrix with Ward's method on square‐root‐transformed tree count data. The au values in black on the left represent unbiased P values supporting each cluster, while bp values on the right represent a bootstrap probability value. The three outlined clusters are based on an alpha > 0.9 cutoff point, as well as on the results of NbClust function for determining the best number of clusters.

Hierarchical cluster analysis of plots also generated a three cluster solution based on the NbClust determinations. Of the 16 species selected by the multipatt algorithm for assessing indicator species, 13 spp. were associated with just one cluster and three spp. were associated with two clusters. The first cluster had only one significant indicator species (P. hartwegii) and no other associated species (Multipatt, stat = 0.942, P = 0.005). The second cluster had four significant indicator species (Q. obtusata, P. pseudostrobus, Litsea glaucescens, and Q. glabrescens; Multipatt, stat > 0.394, P < 0.01), and six other associated species that were not significant indicators (understory specialists of Frangula capreifolia, Myrsine juergensii, Buddleja cordata, Cercocarpus macrophyllus, Cleyera integrifolia, and Oreopanax xalapensis). The third cluster had eight significant indicator species (P. douglasiana, Q. scytophylla, P. teocote, Q. conzattii, P. oaxacana, P. patula, Q. laeta, and Q. affinis; Multipatt, stat > 0.340, P < 0.025) and one associated species (P. leiophylla). Four species were associated with both the first and second clusters: Q. laurina, Q. rugosa, and Q. craossifolia as significant indicators (P = 0.005) and P. ayacahuite (at the P = 0.09 significance level). Three species were associated with both the second and third clusters, although none were significant indicators (Alnus acuminata, Q. candicans, and Prunus serotina).

Niche‐differentiation of tree species

Non‐metric multidimensional scaling analysis showed strong differences in species composition across the altitudinal gradient and soil gradients (non‐metric multidimensional scaling with Bray–Curtis distance, stress value = 0.167, non‐metric fit R² = 0.97; Fig. 6). Tree species distributions significantly correlated with gradients of altitude and soil pH (P < 0.001), slope and soil nitrogen content (P < 0.01), and the percent content of organic material and percent carbon content (NMDS, R² > 0.09, P < 0.5). Based on a visual assessment, the three clusters of plots generated by the hierarchical cluster analysis grouped well together (Fig. 6a). The first cluster largely separated out along the elevation axis and was located at highest altitudes in the NMDS, while the second cluster contained that largest number of plots and was associated with mid‐elevations and higher soil nutrient contents, and the third cluster was associated with lower elevations, steep slopes, and higher soil pH.

6 Fig.. Two‐dimensional diagrams showing the three‐dimensional NMDS using Bray–Curtis distance of tree species densities by plot. Plots are represented by a circle and are color‐coded by cluster based on the hierarchical cluster analysis. (a), (b), and (c) show different possible axis combinations to show the separation of plots and species from different angles. Species names are located where their densities are highest, and arrows show the statistically significant fitted environmental variables (altitude, slope, soil pH, nitrogen, and carbon content). Arrow length and direction signify the strength and value of the correlation. Carbon content and organic material were highly correlated, so only percent carbon is shown on the figure for clarity. In (d), lines connect members of the three clusters to show their grouping in a three‐dimensional rendering.

On a species level, P. hartwegii was correlated with highest elevations compared to other species, and among high elevation oaks, Q. crassifolia was located at higher elevations compared to Q. rugosa or Q. laurina (Fig. 6b). At lower elevations, P. teocote and Q. conzattii were associated with a higher pH (less acidic soils) but with lower soil nutrient content of both carbon and nitrogen compared to P. oaxacana and Q. scytophylla. Unlike in the hierarchical cluster analysis of species, P. leiophylla was found associated with other low elevation pine species of P. douglasiana, P. patula, and P. teocote. Similarly, A. xalapensis associated with both high and low elevation communities in the ordination, whereas it clustered with the low elevation group in the hierarchical cluster analysis. Finally, high organic matter, and carbon and nitrogen content were associated with a community of hardwood species, including Q. glabrescens, L. glaucescens, P. serotina, F. capreifolia, Cleyera integrefolia, and B. cordata.

The multivariate model of the top tree species found that altitude, slope, the depth of the organic layer, the percent organic material, and the N and C contents were significant factors for species composition (Table 2; MANOVA, df = 16 and 66, Pillai's trace > 0.31, F > 1.86, P < 0.05). Univariate responses for each species in the model found that most species were significantly correlated with altitude, with the exception of a few species with broader distributions or distributions at middle elevations (i.e., the understory species A. xalapensis and L. glaucescens, the oak species Q. obtusata and Q. laurina, P. ayacahuite, and A. acuminata). Although pH and aspect were not significant factors for the multivariate species composition, several species were significantly correlated with pH at the P < 0.1 level (P. ayacahuite, P. douglasiana, and Q. obtusata) and two species (P. patula and P. teocote) were significantly correlated with aspect in the univariate models (Table 2; ANOVA, P < 0.05).

2 TableResults from multivariate regression analysis of top tree species density by altitude, slope, aspect (northness), and soil variables pH, depth of organic layer, percent organic material, and percent N and C content.

Species	Altitude	Slope	Northness	pH	Depth organic layer (cm)	Organic material (%)	Nitrogen (%)	Carbon (%)
Model results (Pillai's trace)	0.87 (***)	0.33 (*)	0.21	0.2	0.51 (***)	0.30 (†)	0.33 (*)	0.37 (**)
Univariate responses (linear coefficients)
A. acuminata	2.42 × 10−4	0.01 (*)	0.11	0.13	0.06	1.93 (†)	−30.81	17.87 (*)
A. xalapensis	8.44 × 10−5	0.00	−0.60	−0.35	−0.11	1.66	45.89	−26.60
L. glaucescens	1.25 × 10−3	0.01 (†)	0.08	0.43	−0.06	4.42	34.73	−20.14
P. ayacahuite	2.15 × 10−3	0.01	0.21	0.42 (†)	−0.28 (***)	6.17	88.07	−51.11 (*)
P. douglasiana	−6.50 × 10−4 (***)	−0.05	−0.33	2.55 (†)	−0.23	−57.29 (**)	10.81	−6.15
P. oaxacana	−6.21 × 10−3 (†)	−0.01	0.42	−1.12	0.04	19.78 (**)	34.24	−19.89
P. patula	−3.27 × 10−4 (**)	−0.03	0.28 (*)	0.97	−0.10	−30.49	−42.27 (***)	24.66
P. pseudostrobus	2.90 × 10−4 (*)	0.01	0.22	0.07	−0.15 (*)	−8.51	−190.30	110.4 (**)
P. hartwegii	1.45 × 10−3 (***)	−0.01	0.41	0.31	0.62 (**)	4.77	119.40	−69.31
P. teocote	−2.62 × 10−4 (†)	−0.01	−0.86 (*)	0.48	−0.13 (*)	4.71	56.09	−32.55
Q. conzattii	−4.87 × 10−3 (***)	−0.04 (**)	0.05	0.66	−0.11	−15.47 (*)	20.38	−11.79
Q. crassifolia	5.41 × 10−3 (*)	0.03 (†)	0.13	0.42	0.05	−12.38	21.02	−12.17
Q. laurina	1.52 × 10−4	0.06 (*)	1.17	0.68	0.48 (*)	−8.17	−127.00	73.69
Q. obtusata	3.99 × 10−4	0.03 (***)	−1.09	−0.86 (†)	−0.04	−8.16	−84.37 (**)	49.09
Q. rugosa	5.42 × 10−3 (***)	0.01	0.40	−0.58	−0.12	−14.35	−83.94	48.73
Q. scytophylla	−6.57 × 10−5 (*)	0.04	−0.55	−0.18	−0.30 (*)	17.81	100.1 (†)	−58.17

Notes:

Results for the multivariate regression are shown first as Pillai's traces, and univariate results for each species are presented with linear coefficient estimates. For the MANOVA, tree counts were square‐root‐transformed, but linear coefficients are based on untransformed tree densities.

*** P < 0.001, ** P < 0.01, * P < 0.05, † P < 0.1.

Tree species co‐occurrence

The log‐likelihood test of independence (G test) of distributions of the count of red section oaks and white section oaks in each plot found that there was a significant relationship between these clades (log‐likelihood ratio G test, G = 951.33, df = 89, P < 2.2 × 10⁻¹⁶). Overall, 60% of the plots contained at least one oak from both the red and the white sections, while only 32% of plots contained oaks from just one or the other clades, and 7% of plots had no oaks present.

Discussion

Diversity patterns were striking in this forest, in that over half of the 32 tree species measured were either pine (eight species) or oak (10 species), leading to relatively low diversity of genera (13 in total) but high diversity of species within just two genera. We found that species composition was significantly correlated with altitude, topography, and soil characteristics (particularly soil nutrients and organic matter content) and that species diversity patterns shifted along elevation zones with increased richness and diversity at lower elevations. Tree communities were delineated into three broad vegetation types that aligned with a high elevation, a mid‐elevation, and a low elevation community with significant unique pine and oak indicator species in each. Finally, we found a notable pattern of phylogenetic overdispersion in oak distribution, with a high mean richness of oak species in each plot and a significant correlation between the presence of species of red and white clades. These results have important implications for future forest management of these species, as well as for potential impacts of climate change‐induced drought on these species distributions and vegetation communities.

Tree species diversity

The species richness encountered in the study area is high compared to reported rates for other montane forests in Mexico, especially considering that the diameter limit for sampling (10 cm) was higher than that used in many other studies. For example, studies have found similar overall species numbers but smaller numbers of oak species in the neighboring community of Ixtlán de Juarez (Castellanos‐Bolaños et al. 2010), in La Chinantla, a hyper‐humid forest area in Oaxaca (Meave et al. 2006), in pine–oak forest in Chiapas (González‐Espinosa et al. 1991), and even in the Sierra de Manantlan Biosphere Reserve in western Mexico (Olvera‐Vargas et al. 2010). Considering this last study took place in a national park and used a smaller diameter limit, it is impressive that the community‐managed forest of our study area is maintaining comparable levels of tree diversity while also harvesting timber and fuelwood products from all eight pine species and nine of the 10 oak species.

Across the three elevation zones sampled, we found distinct patterns in species richness and diversity that correspond with the way climate and environmental conditions shift along the altitudinal gradient. We found that both overall species richness and beta‐diversity measures increased at lower elevations such that the low zone was both more species rich and more spatially heterogeneous. Both the low and middle zones showed higher rates of turnover as well as lower rates of nestedness, which aligns with the overall higher number of species (twice as many as in the high zone) and higher richness per plot in these zones. In contrast, species nestedness was markedly higher in the high elevation zone, indicating that differences in composition between plots here were largely due to a loss of species at poorer sites rather than additions of new or rare species. Declining beta‐diversity with altitude has been documented in other studies in tropical dry forest in Oaxaca (Gallardo‐Cruz et al. 2009), in Jalisco, Mexico (Vázquez and Givnish 1998), and in other regions (Wang et al. 2002), with varying explanations. As Vázquez and Givnish (1998) note, there is often a higher incidence of drought, fire, and human disturbance at lower elevations in Mexico, which may lead to higher heterogeneity. The low zone also seemed to contain more varied topography, with riparian slopes around streams and dry ridges.

While overall richness was highest in the low zone due to the higher species turnover rates there, mean species richness per plot and species diversity were highest at mid‐elevations, following the humped pattern seen in many other regions (Rahbek 1995, Wang et al. 2002, Grytnes 2003, Kessler et al. 2011, Sanders and Rahbek 2012). The middle zone receives high levels of orographic precipitation and therefore seems able to support a more stratified canopy with more understory specialists per plot (Gentry 1982, Condit et al. 1996, Vázquez and Givnish 1998, Givnish 1999), while the low water availability at low elevations and the low temperatures and low nutrient cycling at high elevations are likely limiting in those zones (Marrs et al. 1988, Vázquez and Givnish 1998, Kessler et al. 2011).

Species coexistence and phylogenetic dispersion

The oak diversity in this study area was noteworthy, with multiple oak species per plot (a mean of 2.3 species and a maximum of six species recorded per 400 m²) and at least one oak from both the red and the white sections present in the majority of plots (the presence of the two clades were significantly correlated). Other research on oak communities has also demonstrated this pattern of phylogenetic dispersion where co‐occurring oaks are more likely to be from different sections of the genus (i.e., red and white oaks together) than from the same section (i.e., only red or only white oaks; Whittaker 1969, Mohler 1990, Cavender‐Bares et al. 2004, Cavender‐Bares 2018).

The mechanisms that maintain species diversity have long been a focus of ecological research and remain unclear—especially among closely related species such as oaks (Cavender‐Bares 2018). Such a diversity pattern may be partially explained by small‐scale site heterogeneity leading to niche‐partitioning (MacArthur and Levins 1967, Chesson 2000) through sympatric parallel radiation in drought adaptations across the two oak clades (Cavender‐Bares et al. 2004, Poulos et al. 2007, Olvera‐Vargas et al. 2010, Cavender‐Bares 2018, Fallon and Cavender‐Bares 2018). Significant differences in physiological and leaf structural characteristics have been found in Central America in oak species (Poulos et al. 2007, 2008, Aguilar‐Romero et al. 2017), and across populations (Ramírez‐Valiente and Cavender‐Bares 2017), and arguably demonstrate a functional mechanism driving oak differentiation and diversification in the region. In addition, oaks can display differences related to shade tolerance and successional roles in Mexico (Poulos et al. 2008, Rivas‐Rivas et al. 2017), which may also contribute to niche‐partitioning.

Another mechanism that can contribute to the coexistence of similar species is the storage effect (Chesson and Warner 1981, Usinowicz et al. 2012). The storage effect represents resource partitioning through time, where species with long‐lived populations which show temporal fluctuations in reproduction can drive high diversity patterns (Chesson 2003, Angert et al. 2009). Mast seeding (i.e., variable production of seeds, sometimes synchronized, at intervals of two or more years) has been documented for pine species in the United States and Mexico (White 1985, Savage et al. 1996, Mast et al. 1999, Barton et al. 2001, Barton 2002, Owen et al. 2017) and oak species across North America in temperate areas (Sork et al. 1993, Koenig et al. 1994), tropical zones (Sork 1993), and in the southwestern United States and Mexico (Olvera‐Vargas et al. 1997, Alfonso‐Corrado et al. 2007, Parmenter et al. 2018). Modeling indicates that temporal variability and asynchrony in seed and seedling production can translate into large reductions in interspecific competition and therefore into increased likelihood of coexistence (KelIy and Bowler 2002, Usinowicz et al. 2012), which may contribute to the high diversity of species in both of these genera in the study area. Additionally, many of these coexisting oak species show alternate flowering and fruiting in either the dry or rainy season, which may further contribute to diminished temporal competition and a resulting storage effect.

Tree species distribution and community differentiation

Despite the ecological and economic importance of pine and oak in Mexico, relatively little is known about the floristics of most species. We found that species composition was significantly correlated with altitude, topography, and soil characteristics (particularly soil nutrients and organic matter content), similar to results of other studies in Mexico where altitude is the most significant explanatory factor for Quercus (Meave et al. 2006, Olvera‐Vargas et al. 2010) and Pinus species distributions (Yeaton 1982, Alba‐López et al. 2003, Ríos‐Altamirano et al. 2016). In general, altitude‐related moisture availability is a well‐established factor driving differentiation of pine and oak in montane regions around the globe—including in Costa Rica (Kappelle et al. 1995), the Himalaya region in India (Chawla et al. 2008, Singh et al. 2009, Saha et al. 2016), and in temperate United States (Abrams 1990, Kubiske and Abrams 1992, Ashton and Berlyn 1994, Fallon and Cavender‐Bares 2018). Although a study of species composition in the Sierra Madre Oriental in northern Mexico found that aspect (and the resulting moisture differences between slopes) was an important predictor of species richness and floristics (González‐Tagle et al. 2008), we found that northness was only a significant predictor for two individual pine species (P. patula and P. teocote) and not for broader species composition or vegetation clusters. This may be due to the more southern location of this research site, as fewer differences in solar radiation between northern‐ and southern‐facing slopes likely occur at this latitude.

Overall, we found three main vegetation types in this study area that corresponded with high, middle, and low elevations as well as differences in soil characteristics. There were also differences in species groupings between the two cluster methods (i.e., by species or by plot) which indicate that subdivisions may exist within these broader categories. For example, both clustering methods were consistent in grouping low elevation pines and oaks together (P. douglasiana, P. teocote, P. oaxacana, P. patula, Q. conzattii, Q. scytophylla, Q. laeta, and Q. candicans; Fig. 7b, c), but both the hierarchical clustering of species and the NMDS showed separation between P. oaxacana and the other low elevation pine species. Observationally, many P. oaxacana stands in the area appeared to be structurally distinct and sometimes managed quite differently (i.e., for pine resin rather than for timber and fuelwood). The stands dominated by P. teocote and P. douglasiana, on the other hand, seemed to be often found in areas that local community members reported to have been previously used for agriculture, indicating that this subgroup of species may the result of old‐field succession in this region.

7 Fig.. Examples of each vegetation type and subgroup differentiated in the cluster analyses, including (a) the high vegetation type dominated by Pinus pseudostrobus with high elevation oak components, (b) the low vegetation dominated by Pinus oaxacana with oak components, (c) the low tree group dominated by Quercus scytophylla with Pinus teocote and Pinus douglasiana components, (d) the highest vegetation community dominated by Pinus hartwegii, and (e) the mid‐elevation hardwood vegetation type dominated by oaks and understory specialists. Photo credits: a–d, M. Martin; e, H. Haas.

In the high and mid‐elevations, there were also a few differences in species groupings between clustering methods. In the clustering of species, P. pseudostrobus and Q. obstusata were grouped with P. hartwegii and other high elevation oaks, but in the clustering of plots with indicator species method, P. pseudostrobus and Q. obtusata were grouped with the mid‐elevation hardwoods and understory specialists and P. hartwegii formed its own distinct community type at highest elevations (Fig. 7a, d). Likewise, the mid‐elevation cluster of hardwood and understory specialists seemed to form a clearly separate sub‐vegetation type in the NMDS and the clustering of species that was primarily located in more mesic areas at mid‐elevations, with occasional pine elements from P. pseudostrobus or P. ayacahuite (Fig. 7e).

Species niche differences, drought tolerance, and climate change

Differential drought sensitivity shapes tree species' distributions at both local and regional levels, and across wet and dry topographic sites (Engelbrecht et al. 2007). One of the few studies to examine drought response of two Mexican oak species argues that there is a strong relationship between a species' adaptation to water stress and its altitudinal distribution, with a lower elevation oak found at drier sites showing greater drought tolerance (Poulos et al. 2007). In tropical forests, species found on wetter topographic sites tended to be more sensitive to drought compared to those distributed on drier sites (Engelbrecht et al. 2007), so species here that are found in more mesic conditions (e.g., P. patula, P. ayacahuite, Q. laurina, Q. glabrescens, and associated hardwoods) may actually be more sensitive to drought, and therefore to climate change, compared to those found in more xeric areas (e.g., P. teocote, P. douglasiana, Q. crassifolia, Q. scytophylla, Q. conzattii, and Q. laeta).

High elevation species

Often species found at higher elevations are considered more vulnerable to shifts in climate (Chen et al. 2011, Laurance et al. 2011, Freeman et al. 2018). The two pines in the high elevation community, P. hartwegii and P. pseudostrobus, are both considered among the most vulnerable of Mexican pine species to climate change based on range estimates from herbarium records due to their restricted distributions and high altitudinal ranges (Gómez‐Mendoza and Arriaga 2007, Galicia et al. 2015). Of these, P. hartwegii showed a more restricted altitudinal distribution and therefore may be more vulnerable to climate shifts as it would not have the space to migrate upslope as temperatures increase.

The oak species found here at higher elevations (Q. rugosa, Q. crassifolia, and Q. laurina) are considered moderately sensitive to climate change according to range estimates from herbarium collections records (Gómez‐Mendoza and Arriaga 2007), but these climate sensitivity studies were conducted from a large‐scale perspective and may be missing some of the differentiation between species that we see here. For example, of these three oaks, Q. crassifolia was more common at slightly higher elevations and on steeper slopes based on our MANOVA and NMDS results and therefore may be more drought‐tolerant compared to the other high elevation oak species (Engelbrecht et al. 2007, Poulos et al. 2007, 2008), but may also grow in areas that are more susceptible to increased drought as climate shifts.

Middle elevation species

In the mid‐elevation group, P. ayacahuite (a soft pine) is considered a moderately sensitive species (Gómez‐Mendoza and Arriaga 2007, Galicia et al. 2015) and is often found at more humid sites compared to P. pseudostrobus (a hard pine), which is associated with drier areas, warmer annual temperatures, and poor soil quality (Galindo‐Jaimes et al. 2002, Alba‐López et al. 2003). In our analyses, P. ayacahuite was associated with hardwood understory species in both clustering analyses and in the NMDS, likely in more mesic areas, and appears to be both less drought‐sensitive and less resistant to fire (Rodríguez‐Trejo and Fulé 2003, Poulos et al. 2018).

For the oak species that dominated the more mesic areas at middle elevations, herbarium record estimates argue that Q. obtusata should be considered a tolerant species (Gómez‐Mendoza and Arriaga 2007). This aligns well with our assessment, where Q. obtusata showed a wide range in vegetation associations, implying that it may be a more generalist species; and indeed, it has one of the widest geographic distributions and one of the widest altitudinal ranges among oaks in Mexico (Valencia 2004). Quercus obtusata (locally called roble) is not a preferred species for fuelwood—community members report that the wood is too hard to burn well—so a lack of harvesting may also allow this species to spread into a broader range of habitats compared to other, more preferred oaks. The impact of climate change on these species will be complicated by these differences in management and harvesting, and as seen with Q. obtusata, disentangling the effects of these two drivers can be challenging.

Low elevation species

While two of the oaks found here at lower elevations (Q. scytophylla and Q. laeta) are considered tolerant to climate change (Gómez‐Mendoza and Arriaga 2007), four of the oak species found in the study (i.e., Quercus affinis, Q. candicans, Q. conzattii, and Q. glabrescens) have not been studied from a climate modeling perspective. Of the lower elevation pine species found in this cluster, P. patula, P. teocote, and P. leiophylla are considered moderately sensitive to climate change, while P. douglasiana is considered tolerant (Gómez‐Mendoza and Arriaga 2007, Galicia et al. 2015). Of these pines, P. patula had the broadest altitudinal range in our study, a pattern which was also seen in another study of Oaxacan pines (Ríos‐Altamirano et al. 2016). P. patula is associated with more mesic sites in another studies (Rodríguez‐Trejo and Fulé 2003), and in our model was one of the few species to be significantly correlated with northern aspects, which are likely to be slightly more mesic compared to southern‐facing slopes at latitudes above the equator. In contrast, P. teocote distribution was correlated with more droughty sites on southern aspects and with a thinner organic layer, and therefore may be more drought‐tolerant.

Finally, it would be remiss to discuss pine and oak species distributions without mentioning the role of fire, especially as fire, drought, and climate change are so closely entwined (Liu et al. 2010, Abatzoglou and Williams 2016, Marín et al. 2018). Almost all of these pine and oak species in this study show multiple adaptations to fire, including serotinous cones, thick bark, self‐pruning, and resprouting abilities (Rodríguez‐Trejo and Fulé 2003, Rodríguez‐Trejo and Myers 2010, Galicia et al. 2015). As noted by Poulos et al. (2018), many of these pine species demonstrate multiple fire adaptation traits that are associated with different fire regimes, possibly as an evolutionary bet‐hedging strategy favored after exposure to mixed regimes over time. For example, P. hartwegii, P. pseudostrobus, P. teocote, and P. leiophylla have traits associated with fast, hot, and non‐flammable regimes and P. patula can be associated with fast and hot‐flammable, while only P. ayacahuite was associated with non‐flammable alone (Poulos et al. 2018). Observationally, fire appeared to have been more commonly used in management and agricultural practices at lower elevations, but without knowledge of fire histories of specific locations, it is difficult to interpret how these adaptations align with fire regimes and community differentiation.

Conclusions

This community‐managed forest contains high levels of tree diversity, with especially high levels of diversity of oak compared to other pine–oak forests so far studied in Mexico. This diversity is noteworthy given that well over half of the species encountered are actively harvested by the local community—all of the pine species are harvested for timber, eight of the oak species are preferred for fuelwood. This study found three main community assemblages, which primarily differentiated across the landscape along the altitudinal gradient, with subgroupings of species potentially differentiated by soil characteristics and differences in land‐use history and management. Oak in the study area demonstrates high levels of beta‐diversity across these community types, as well as high levels of species coexistence within communities, which followed a pattern of phylogenetic overdispersion with species from different sections (red and white) co‐occurring more frequently than species from the same section. This diversity pattern could be caused by small‐scale niche‐partitioning in drought tolerance or shade tolerance traits, or through temporal niche‐partitioning (a storage effect) due to masting dynamics or to oak species alternately fruiting across dry and rainy seasons.

Unfortunately, little research has examined either physiological and anatomical leaf and wood drought adaptations or differences in reproductive phenology and masting for these particular oak species in Oaxaca. Such dynamics also have important implications in the face of climate change, as understanding of species' drought adaptations will be key to future management, as will knowledge of species fruiting phenology. In particular, oaks fruiting during the dry season may be subject to increased drought pressures. Mexico becomes more arid and may suffer reproductive failure and population declines relative to species fruiting during the wet season, as dry season seedling mortality appears to be a key driver of tree species distributions (Comita and Engelbrecht 2009). New research suggests that co‐occurring oaks may benefit each other more than previously understood through mechanisms such as the promotion of ectomycorrhizal symbionts (Desai et al. 2016), by enhancing soil conditions (Chávez‐Vergara et al. 2016), or other forms of facilitation (Cavender‐Bares 2018), underscoring the importance of understanding how these diverse communities are assembled and maintained as a basis for conserving and managing these forests going forward.

Finally, further assessment of how land‐use history has impacted stand dynamics in these forest types is essential. Several of these pine and oak species showed potential impacts of differing harvest preferences or links to agricultural abandonment, and so understanding the ways in which these land‐use practices may interact with future climate shifts in altering species ranges will be key to sustainable management and conservation going forward.

Acknowledgments

Thank you to the funders of this research: the Lewis B. Cullman Fellowship, the Yale Institute of Biospheric Studies, the Tropical Resources Institute, the Garden Club of America, and IdeaWild. We are grateful to our local collaborators at IXETO, Heriberto Aguirre Diaz and Gliserio Marin Garcia and to the community of Nuevo Zoquiapam for their permission and support to work on their lands. Thank you also to Gerardo Rodriguez Ortiz and the Instituto Tecnologico del Valle de Oaxaca for logistical and herbarium support and to field assistants Yazmín Pérez Alavez, Iván José Aguilar Pinacho, Humberto Caceres, Wenceslao Robles, and Israel Betata. We also thank Danica Doroski for her helpful comments on an early draft of this manuscript and the two reviewers for their invaluable feedback.