Small-Scale Environmental Drivers of Plant Community Structure and Diversity in Neotropical Montane Cloud Forests Harboring Threatened Magnolia dealbata in Southern Mexico

: Gradient analysis was used to determine factors driving small-scale variation of cloud forest communities harboring Magnolia dealbata , a threatened species and bioculturally relevant tree for the Chinantecan, Mazatecan, Nahuan, and Zapotecan ethnicities in southern Mexico. Particularly, we aimed to: (a) determine factors explaining major community gradients at different heterogeneity scales along a small-scale elevational gradient, (b) test the Decreasing and the Continuum hypotheses along elevation, and (c) classify vegetation to assist in identifying conservation priorities. We used a stratified random sampling scheme for 21 woody stands along a small-scale (352 m) elevational transect. Four main data matrices were used (presence-absence, density, basal area, and guild data). Through Non-metric Multidimensional Scaling (NMS), Principal Coordinates Analysis (PCoA), and distance-based Redundancy Analysis (db-RDA), we found that major community variation was explained by soil pH, displaying an outstanding vegetation discontinuity, separating the species-rich relic Oreomunnea-Ticodendron -stands from stands with higher importance values for M. dealbata . The high species richness observed was explained by a combination of the windward effect of dry-seasonal maximum cloud condensation gain and habitat differentiation-specialization, a phenomenon that may also explain the mid-peak hypothesis and ensure the survival of relic species. Sampling-truncation and conservation status also played a role in this. Our results do not support the Decreasing and Continuum hypotheses along elevation. Investigation: R.D.-Y., J.A.V.-G., M. Á E.S.-P., C.R.-P. and S.I.G.-Y.; Methodology: R.D.-Y., J.A.V.-G., M. Á C.R.-P. and S.I.G.-Y.; Project administration: J.A.V.-G.; Resources: J.A.V.-G., M. Á and E.S.-P.; Software: R.D.-Y. and J.A.V.-G.; Supervision: J.A.V.-G. and M. Á Validation: J.A.V.-G., M. Á and G.H.-V.; Visualization: R.D.-Y., J.A.V.-G. and M. Á Writing—original draft: R.D.-Y. and J.A.V.G.; Writing—review and editing: R.D.-Y., J.A.V.-G., M. Á G.H.-V., E.S.-P., C.R.-P. and S.I.G.-Y. All authors read and


Introduction
Neotropical montane cloud forests (CF) are known to display high gamma, beta, and alpha species diversity, environmental heterogeneity, endemism, and a high number of species of conservation concern. They form a continental pattern of island-like fragmented habitats resembling archipelagos [1][2][3][4][5][6][7][8]. Despite their biological relevance, these complex ecosystems are being replaced or in the central orographic axis of the physiographic sub-province Sierra Madre de Oaxaca [12], particularly characterized by steep topography (15-100%) [76]. Major rocks in Juquila Oaxaca include igneous intrusive: monzonite (12.52%), sedimentary: lutita-sandstone (2.73%), metamorphic slate (44.67%), and shale (40.08%) [73]. Soils include luvisol (85.35%) and cambisol (14.65%) and are mostly acidic [77][78][79]. Juquila Vijanos includes the Juquila and Río Blanco rivers, the two drain into the Cajonos River, and finally reach the Papaloapan River. Climate is semi-warm humid with mean annual temperatures ranging from 16 to 22 • C, with abundant mean annual rainfall of 2600 mm (San Juan Yae, 1900 m a.s.l.), with 69.34% of rainfall occurring during the summer [73]. Due to the abundance of temperate elements and diagnostic families such as Actinidiaceae, Clethraceae, Chloranthaceae, Hamamelidaceae, Symplocaceae, Theaceae, and Winteraceae [76], the original vegetation is considered an upper Tropical Montane (TM)-CF [80]. Cloud forests in the Sierra Norte de Oaxaca occur between 1300 and 2300 m a.s.l., mostly on the Atlantic windward side of the mountains [14], on hilly sides with pine-oak forest and "acahuales" (secondary successional vegetation) of different ages resulting from the abandonment of agricultural lands [81]. However, much of the original vegetation has been converted to coffee plantations, agriculture (including polycultures), and some pasturelands [78,81]. Oaxaca is an important center of domestication of Mesoamerican plants, including currently complex agroforestry systems like those in Sierra de Juárez (Chinantla), where useful species are tolerated within field crops or rustic coffee plantations involving up to 34 native tree species, including threatened species of Magnolia, which serve as nurse plants. This type of forest management with high beta diversity has important implications for biodiversity conservation and sustainable development [13]. In the area prevail allogenic disturbances such as logging, slash and burn agriculture, firewood, and non-timber product extraction; meanwhile, autogenic disturbances, such as landslides, fires, and natural tree falls, occur to a lesser extent [78,81,82]. In contrast with other successional processes in humid forest areas, soil becomes progressively more acidic along succession, despite a decrease in pine dominance [78,82].

Community Sampling
A stratified random sampling including 21 0.1 ha forest stands (sites) within the habitat of Magnolia dealbata was carried out from December 2016 to August 2017, in San Juan Juquila Vijanos, on areas with minimum disturbance and excluded from cattle grazing. A 60 × 48 m rectangular universe was defined, divided into a 4 × 5 grid (each cell 12 × 12 m). At the center of each cell, a circle was delimited to include an area of 100 m 2 per quadrat (subplot). The centers of the possible circular quadrats (11.28 m in diameter) were separated by 12 m, to avoid any overlap between the adjacent quadrats. Ten circular quadrats per site were chosen in a stratified random sampling design [56]. In each circle, the diameters of all woody species ≥2.5 cm in diameter at breast height (dbh, at 1.3 m) were measured and recorded. A triplicate collection of voucher plant specimens was obtained: they were mostly identified by the authors, and difficult taxa were determined by specialists. For accepted names of species, we used Tropicos.org [83], and for plant distributions, we consulted Plants of the World Online (POWO) [84]: their nomenclature and author abbreviations follow the International Plant Name Index [85]. The main voucher collection was deposited in the herbarium of the Instituto de Botánica de la Universidad de Guadalajara (IBUG) herbarium at the University of Guadalajara and sets of duplicates were deposited at Centro Interdisciplinario de Investigación para el Desarrollo Integral Regional (CIDIR) and the herbarium Sociedad para el Estudio de los Recursos Bióticos de Oaxaca, A. C. (SERO).

Community Sampling
A stratified random sampling including 21 0.1 ha forest stands (sites) within the habitat of Magnolia dealbata was carried out from December 2016 to August 2017, in San Juan Juquila Vijanos, on areas with minimum disturbance and excluded from cattle grazing. A 60 × 48 m rectangular universe was defined, divided into a 4 × 5 grid (each cell 12 × 12 m). At the center of each cell, a circle was delimited to include an area of 100 m² per quadrat (subplot). The centers of the possible circular quadrats (11.28 m in diameter) were separated by 12 m, to avoid any overlap between the adjacent

Environmental Sampling
At each site, a total of 31 environmental variables were recorded in five subsets: Climatic (5), topographic (4), edaphic (14), canopy structure (5), and disturbance (5) ( Table 1). Highly correlated variables (r > 0.8) of each subset were excluded from the analysis to avoid collinearity. Climatic data were extracted from the WorldClim database [86] and an index of the heat load was estimated [87]. Topographic variables were directly measured in the field and elevation was recorded by a Global Positioning System GPS Garmin 60CSx. Edaphic variables were obtained from soil analyses. One mixed Diversity 2020, 12, 444 6 of 29 soil sample (500 gr) per site, was taken from five quadrats of 30 × 30 cm each, and 30 cm deep (one quadrat in each site's corner and one in the center) [88]. The samples were analyzed in the Agroecology Laboratory of the Centro Universitario de Ciencias Biológicas y Agropecuarias, University of Guadalajara. Soil density was determined through the test tube method, and sand, silt, and clay texture were determined by Bouyoucus. Usable water was determined by the difference between the soil field's capacity and permanent wilt point. Organic matter was estimated through the Walkley and Black [89] (1934) method; for interchangeable cations (Ca), we used the volumetric method. Mg was inferred from the volumetric method, Na and K were assessed by flammometry, and soil pH was measured by a potentiometer; for electrical conductivity, we used a conductometer, and for soil macronutrients, phosphorus (P), and percentage of nitrogen (N), we used the Kjeldahl method [90]. Canopy structure variables were obtained from hemispheric photographs at each site, using a Nikon D7100 digital camera and a hemispherical (fisheye) Canon lens. To reduce the effects due to the slope, in each site, the camera was mounted on a tripod at a height of 1.20 m [91] and a level device was used to ensure horizontal alignment of the camera. Each photograph was taken at the center of each site and oriented to match the top of the photograph with the north [91,92]. Photos were taken on cloudy days, to prevent the sun's rays from directly affecting the camera lens. Each canopy photograph was taken from a leveled camera [93]. Images obtained were analyzed through the Gap Light Analyzer software (GAP free version 2.0, 1999), to correct the distortion inherent in this type of photographs. The disturbance was assessed, assigning a score from 0 to 10 to each variable, depending on the intensity of the disturbance: a score of 0 was considered as insignificant, 1 as low, 5 as intermediate, and 10 as high. A disturbance index was calculated following Mir et al. [94].

Ordination
We used Non-metric multidimensional scaling (NMS), a distance-based ordination considered superior to currently available model-based methods [95], to map the location of sample units from the high-dimensional space of the distance matrices, to positions along community dissimilarity gradients in a low-dimensional space. In NMS, species similarity between samples represent their environmental similarities but in a highly coalescent configuration, including similarities in intra-(density-dependent mortality, resource competition) and inter-specific interactions (herbivory, dispersal, predation, pollination) and biogeographical events [42]. Four main raw matrices (P/A, density, basal area, and guild data) were used, each consisting of 21 sites and 31 environmental variables. We used the Sørensen (Bray-Curtis) distance measure, having a monotonic relationship to environmental distance, and considered as a robust measure of ecological distance [42,43,96]. The secondary matrix consisted of 21 sites × 31 variables (Table 2). NMS is considered among the most effective methods for the ordination of ecological data because it does not require the assumption of a linear relationship among the variables and because it determines the dimensionality of the data [89]. Preliminary runs were performed to determine the appropriate dimensionality using the autopilot mode, choosing the slow and thorough option, the Sørensen distance, an instability criterion of 0.00001, with 50 runs with real data, and 500 iterations to evaluate instability. We used Pearson's correlation coefficients to evaluate the relationship between the identified axes of the ordination and the environmental variables. Principal Coordinates Analysis (PCoA) was performed only for guild data, since these data showed the lowest gradient length (low heterogeneity), making it appropriate for this metric multidimensional technique. We used a Sørensen distance matrix of 21 sites × 60 guilds for this ordination method, in which the distances among sites in the ordination diagram are maximally correlated with the ecological distances [97]. Additionally, we used distance-based Redundancy Analysis (db-RDA), a reliable constrained ordination technique, to assess the influence of environmental variables [98,99]. The db-RDA allows for the flexible selection of many resemblance functions; again, we used Sørensen dissimilarity in each of the four matrices: P/A, density, basal area, and guild data. The variables most associated (with highest r values) with the NMS axes were selected as input for a standard RDA. NMS, PCoA, and db-RDA were performed using the program PC-ORD 7 [97].

Classification
A hierarchical cluster analysis, was performed using PC-ORD 7 [97], with the flexible beta algorithm, was performed to determine community groups, using a Sørensen distance matrix of 21 sites × 77 species. We used a restriction criterion (β = −0.25) as a linkage method to obtain a dendrogram [95]. Compositional differences among groups, generated by cluster analysis, were assessed using a multi-response permutation procedure (MRPP), a technique not requiring distributional assumptions (multivariate normality and homogeneity of variance) [95]. This technique provides a statistical test (T) describing the separation among groups, a p-value evaluating the likelihood that an observed difference is due to chance, and the agreement statistic A, describing within-group homogeneity compared to that expected by chance. When all the items within groups are identical, the weighted mean within-group distance (δ) = 0 and A = 1, the highest possible value for A. A neutral value (A = 0) is obtained when heterogeneity matches that expected by chance. If there is less agreement within groups than expected by chance, then A < 0 [95]. We run MRPP available on the PC ORD 7 software [97] using the Sørensen distance measure, the same one used in both the ordination and classification. Indicator Species Analysis (ISA) was used to determine the most representative and exclusive tree species in each of the groups detected in the cluster analysis [89], using the relative frequency and relative abundance of the species in each selected group. Indicator values were tested through 1000 Monte Carlo randomizations [95]. ISA was also used to objectively define the optimum number of groups, selecting the hierarchical cluster step with both the smallest average p-value and the highest number of significant indicator species [95,100]. ISA was performed for P/A data, using the Tichý and Chytrý [101] method; for abundance data (density, structure, and guilds), we used Dufrêne and Legendre [100], with the Sørensen distance measure. We used the software package PC-ORD 7 for all multivariate analyses [95].

Diversity and Structure
A regression analysis was carried out using species richness (α-diversity) as the dependent variable and elevation as the independent variable. The resulting pattern for the 21 sites of this study was compared against data of alpha diversity for 146 sites from previous studies using the same sampling scheme with circular quadrats [47]. The relative species-richness of families was examined along the elevational gradient: the 21 sites were grouped and averaged by elevation into four classes: 1600−1699 m a.s.l. (sites 15 −17, 18, 20), 1700−1799 m a.s.l. (sites 5, 10−11, 13−4), 1800−1899 (sites 1,4,6,9,12), and 1900−1999 m a.s.l. (sites 2-3, 7−8, 19,21). Species turnover (β-diversity) was assessed through detrended correspondence analysis (DCA), performed with PC ORD 7 [97] to quantify the heterogeneity, as gradient length, for each of the four data matrices (P/A, density, basal area, and guilds). This measure is ecologically meaningful, and it was also used to assess the magnitude of discontinuities along the primary gradient in ordination and among clusters in vegetation classification. The total amount of the species found (γ-diversity) was obtained from the 21 sites along the 350 m elevational gradient. Structure data were obtained from summary statistics of individual tree species (frequency, density, and dominance) in the 21 sites, using the PC-ORD software package version 7.0.

Endemism and Conservation
We used Tropicos.org and POWO to assess the geographic distribution and endemism of each species. The risk status and endemism of all species were determined using various criteria: The International Union for Conservation of Nature and Natural Resources Red List [102], the Mexican endangered species act NOM−059 [22], and the Red List of Mexican Cloud Forest Trees [20].

Ordination
NMS suggested three dimensions for P/A data, and two and one for density and basal area, respectively. Axis 1 accounted for a greater proportion of variance for basal area data, followed by density and P/A data (Table 3), axis 1 was consistently explained by pH (Table 4) in all three datasets (P/A, density, and basal area). Axis 2 was explained by elevation and MAT for P/A data and by Na for density data, and none of the measured variables explained axis 3 ( Table 4). The lowest "stress" value was achieved with P/A data and with fewer iterations, followed by basal area and density data (Table 3). Groups overlays (convex hulls), showed no overlap in all three datasets (Figure 2A-C), except for P/A data, where the Siparuna-CF overlapped with the Vismia-CF in axis 1-2. Oreomunnea-Ticodendron-CF were consistently separated in all three datasets along axis 1 and explained by soil acidity. For P/A data, the Vismia-CF was separated from the Siparuna-CF along axis 2 and explained by elevation. For density data, the Zinowewia-CF was separated from the Oreomunnea-Ticodendron-CF along axis 2 and explained by Na and K. For basal area data, all three groups where separated along axis 1 by pH. For guild data, PcoA showed high eigenvalues and accounted for a high total variance, inertia (76.22%) (68.01%, corrected), suggesting the existence of strong environmental gradients in the Juquila river watershed. Axis 1 explained 69.74% and axis 2 explained 6.48% ( Figure 3A). Groups overlays (convex hulls) showed no overlap among all four groups, with three of them separated along elevation in lower (lw-CF), mid-elevation (me-CF), and upper CF, with the latter subdivided and separated by soil acidity into low acidity (up-la-CF) and high acidity (O-T-CF) ( Figure 3B).  The db-RDA showed that the canonical axis 1 was explained by pH in three different datasets (P/A, density, and basal area), except for guild data, which was explained by MAT. Canonical axis 2 was explained by MAT in two datasets (P/A and basal area), by elevation for P/A data, by Na and K for density data, and by pH for guild data. Guild data accounted for a greater proportion of both cumulative explained variance from total variation (Table 5, Figure 4A-D). In summary, the db-RDA confirmed that vegetation change was best explained by pH, MAT, and elevation, and the guild dataset was the most efficient. The db-RDA showed that the canonical axis 1 was explained by pH in three different datasets (P/A, density, and basal area), except for guild data, which was explained by MAT. Canonical axis 2 was explained by MAT in two datasets (P/A and basal area), by elevation for P/A data, by Na and K for density data, and by pH for guild data. Guild data accounted for a greater proportion of both

Classification
ISA, applied to every hierarchical step of the flexible-clustering, determined where to prune each dendrogram, suggesting four groups for P/A data (Table 6, Figure 2D), 3 groups for density data ( Table 6, Figure 2E), four groups for the basal area (Table 6, Figure 2F), and 4 groups for guild data ( Table 6 and Figure 3B). P/A and density data were the only ones that resulted in having species with

Classification
ISA, applied to every hierarchical step of the flexible-clustering, determined where to prune each dendrogram, suggesting four groups for P/A data (Table 6, Figure 2D), 3 groups for density data ( Table 6, Figure 2E), four groups for the basal area (Table 6, Figure 2F), and 4 groups for guild data ( Table 6 and Figure 3B). P/A and density data were the only ones that resulted in having species with significant indicator values, while for basal area and guild data, the analysis failed to highlight species with significant indicator values (Tables 6 and 7). The Oreomunnea-Ticodendron-CF had an outstanding number of species with significant IV's. The MRPP showed homogeneity within groups for P/A, density, basal area, and guild data ( Table 6).

Diversity
Species richness increased with increasing elevation (r = 0.39, degrees of freedom (d.f.) GL 19, p < 0.05) ( Figure 5A), inversely to the expected pattern observed for 167 Mexican CF sites (r = −0.38, GL 165, p < 0.005) ( Figure 5A), with the same sampling scheme. The elevation gradient of CF sites in Mexico with transect sampling of 0.1 ha with the Gentry method was 750 to 2750 m a.s.l., and the richness gradient was 15 to 75 species per 0.1 ha (average of 40 species). Transect richness, in general, for 30 Mexican CF sites tended to decrease with increasing elevation (r = −0.7242, d.f.L 28; p < 0.005) ( Figure 5B). Species turnover or gradient length for ecological guild data for the first axis, as determined by DCA, had an eigenvalue of 0.01, and the shortest gradient length of 0.39 of standard deviation SD, representing a low beta diversity. For the P/A data, the eigenvalue of axis 1 was 0.3 and a gradient length of 2.2 SD, using importance value data the eigenvalue of axis 1 was 0.5, and the gradient length was 2.8 SD, using density data the eigenvalue of axis 1 was 0.6, and the gradient length was 2.9 SD, and using basal area data the eigenvalue of axis 1 was 0.7, and had the largest gradient length was 3.9 SD, representing a high beta diversity or high species turnover with almost no species in common between the endpoints of axis 1.  and Myrsine juergensenii (Mez) Ricketson & Pipoly. Asteraceae was the most species-rich family in most elevation ranges, except at 1800−1899 m a.s.l., where it was matched by Lauraceae, Ericaceae, and Fagaceae. Families showed greater equitability at this elevational range. Ericaceae unexpectedly increased toward lower elevations and Lauraceae increased toward higher elevations (1900−1999 m a.s.l.) ( Figure 6). The family richness in terms of species/genera were Asteraceae (9/6), Ericaceae and Rubiaceae (4/4), Lauraceae and Pentaphylacaceae (4/3).  Mexican CF sites, using circular quadrats, with a fixed sampling universe [56]. Pacific slopes (dots in light blue), Atlantic slopes (dots in yellow), this study (dots and central trend line in red). (B) Thirty Mexican CF sites, using transects: with undefined sample universe (19 sites) [1] and within a defined sampling universe (11 sites, brown circles) [17].

Environmental Community Drivers
Community composition and structure of CF with Magnolia dealbata in the Sierra Norte de Oaxaca may be explained by habitat specialization to soil gradients and altitude. On short elevational gradients, edaphic variables are often more relevant than climatic variables in explaining the strongest community gradients [51,103]. For instance, moisture and topography can change within small scales, impacting species composition, while soil nutrients across topography gradients may result in fertility gradients affecting species dominance. Spatial environmental variation in temperature, light, soil pH, moisture, nitrogen, and altitude at any scale may drive to the occurrence of habitat-specialized plant species [104,105], which leads to changes in species composition even on small scales.
Soil acidity (pH) along the short (352 m) elevational gradient was more relevant than any other measured climatic variables in explaining the variance in each of the three strongest community gradients recovered by NMS (compositional, structural, and dominance). In contrast to other studies in Fagus-CF, out of all the analyzed microenvironmental factors, pH was the only one not correlated with community variation [106]. Studies in the study area report soils with a pH below 5, so they are considered acidic and easily leached [77,78]. Despite the fact that conifers tend to increase the acidity of the soil [77,[107][108][109] and that conifer needles provide lower amounts of Ca, K, and Mg than deciduous species [110], the dominance of Pinus chiapensis was inversely correlated to the acidity gradient. The acidity gradient was instead correlated with precipitation of the driest month (r = 0.59, p < 0.01) and inversely correlated with N (r = −0.48, p < 0.05). It is well known that acidity is also attributed to soil leaching [56,111], which is enhanced under high rainfall and steep slopes, like those reported for our sites, generating a higher concentration of aluminum and displacement of cations [78,112]. A significant number of species that prefer acidic soils grow in high-altitude areas because soil pH is negatively correlated with elevation (higher tendency to acidification at higher altitudes), since at higher altitude, the decomposition of organic matter is slower, and the acidification process is more intense due to higher precipitation [113].
However, elevation was not correlated with the strongest gradient (axis 1) on any of the three ordinations. Soil acidity is common in CF [114]; in Mexico, it has been recorded in both the Atlantic slopes of Veracruz and Oaxaca [58,78,115], and the Pacific slopes of Guerrero and Jalisco [116,117] from early to the late-successional stage and under different conservation and extraction regimes, as shown for Juquila Vijanos in the Sierra Norte de Oaxaca [78,79,110].
However, this is the first time that an acidity gradient explains the strongest CF gradients along a short elevational range. The outstanding discontinuity along the acidity gradient (axis 1), separating two major groups of CF does not support the Continuum hypothesis [118]. The great majority of CF stands (19 sites) occur in acidic soils (pH 4.3-5.15), while the other group includes only two sites occurring in strongly acidic soils (pH 4.4-4.0); site 19 (1951 m a.s.l.), mono-dominated by Oreomunnea mexicana (relative dominance, 18%), and the site 21 (1960 m a.s.l.), co-dominated by Persea liebmannii Mez (14%), Ticodendrum incognitum and Dendropanax populifolius (11% each), and Oreomunnea mexicana (10%). However, the acidity level is not as extreme as shown in other Oreomunnea populations in Oaxaca, with pH values from 3.7 to 3.5 [119]. Monodominance of Oreomunnea mexicana has also been registered from La Esperanza (1600 and 1800 m a.s.l.), Santiago Comaltepec, Oaxaca [14], other nearby areas in the same region [120], and in Santa Cruz Tepetotutla, south of San Felipe de Usila, Oaxaca at 1840 m a.s.l. [121]. It has also been inferred from density data, in one single site east of Cofre de Perote, Veracruz [58]. Codominance of this species has been reported in La Esperanza, in coexistence with different taxa (mostly Quercus spp.), between 1800 and 2050 m a.s.l. [14]. Meave et al. [121] recorded for the first time a forest dominated by O. mexicana and T. incognitum in Santa Cruz Tepetotutla at 1840 m a.s.l. After twenty-one years, the present study is the second record for the Sierra Norte of Oaxaca, where the two biogeographic relicts, O. mexicana and T. incognitum, are reported as codominant in a 0.1 ha site, among the other few dominant tree species. Ticodendron incognitum has been found at high densities (39 trees/0.1 ha) in La Esperanza, Oaxaca (1750 m), Mexico [17], while in Santa Cruz Tepetotutla, the highest density (10 trees/0.1 ha) was recorded at 1840 m in elevation [121]. In San Miguel Tiltepec (418 trees/0.1 ha), it was recorded at 1640 m [15], and in lower densities at the Barva volcano, Costa Rica (14−27 trees ha), at 1750−2000 m [54].
Elevation was of secondary importance in shaping vegetation structure, and composition also changes with altitude. Environmental variables such as temperature, precipitation, potential evapotranspiration, and radiation can play a vital role in determining the distribution of species along the altitudinal gradient [65,[122][123][124]. One of the first patterns observed in nature was the change of species richness in an altitudinal gradient, which was widely accepted to decrease as elevation increased [70,125]. In elevation gradients, depending on the type of vegetation, species richness may be associated with different climatic variables [65]. The fact that elevation explained a secondary gradient and not the first may have to do with the range of the elevation of the gradient. The longer the gradient elevation, the more expected to be the overriding variable, suggesting that for each elevational mountain gradient, there must be a gradient range threshold below which edaphic or topographic variables become more important than elevation to explain community structure [51].

Cloud Forest Types and Relationships
In general, the CF studied in the present work correspond environmentally to the warm-temperate lower-montane moist forests [126]. The high woody species turnover and heterogeneity found in Juquila Vijanos has been documented frequently for CF transects in different biogeographic provinces with contrasting floristic composition, for instance, the three cloud forest communities reported from the calcareous Cerro Grande Manantlán Massif in the Sierra de Manantlán  [127]. We found no records of CF inhabiting soils with low pH values such as those observed in O-T-CF in Juquila Vijanos, Oaxaca. Concerning ISA, indicator species for O-T-CF are associated with low acidity, while indicator species for A-O-CF are associated with high elevation. The absence of indicator species in V-CF suggests that their species have a wide distribution range and could be shared with other neighboring CF.

Alpha Diversity
The effect of elevation on species richness is not universal but rather scale-and context-dependent. Some studies report that species richness depends on the life form considered [56,128]. Four major patterns of species diversity along elevation are described among different organisms [70]: (1) Midpeak is the most common pattern in plants [34,37,101,102], (2) decreasing richness with increasing elevation is the most common for woody plants in tropical montane forest [1,23,[28][29][30]32,34,43,57,110,111], and less common patterns are: (3) low plateau and (4) low plateau with a mid-peak [70]. Our study, showing increasing richness with increasing elevation in CF transects, does not conform to any of the major trends, however, to some extent, our result is related to studies showing the increasing pattern from the lowlands to some mid-elevation peak [129][130][131][132]. Thus, the observed increasing pattern for CF could be considered a sampling artifact. Sampling truncation at the upper portion of an assumed unimodal CF species richness pattern could result in a short gradient displaying a positive relationship. Alternatively, it could be explained in part by a combination of (a) habitat specialization to edaphic conditions along a gradient of high soil heterogeneity, where forests develop on soils ranging from low to high acidity, as evidenced in this study (see below), (b) a higher rainfall and fog from eastern trade winds out of the Gulf of México toward higher elevations, which is expected to influence species richness [1,70], (c) biogeographic refugia, including paleoendemic relics, which may also contribute to the high diversity. For instance, among the sites with to the richest and highest elevation group: the O-T-CF, which are considered Oligocene-Miocene relics [14], and (d) better conservation status on steep slopes that are farther away from towns, e.g., disturbance or successional gradients [77,78,82], which in turn may be correlated to altitude [70]. The resulting diversity values for the Juquila Vijanos gradient increased two-fold  along the elevational gradient, and every value fell within the expected species richness ranges for Mexican CF: 8-42 spp. ha −1 using circular quadrats ( Figure 5A) and 15-75 spp. ha −1 using transects ( Figure 5B). Only two sites (19 and 21) showed unusually high values (32 and 35 spp., respectively) for their elevational zone, being greater than that reported for Neotropical CF [1,17,27]. These two sites correspond to the richest and highest elevation group, the O-T-CF.

Family Shifts Along Elevation
The top five species-richest families at the Juquila, Oaxaca gradient, were Asteraceae, Ericaceae, Fagaceae, Lauraceae, and Rubiaceae. The latter two are also among the top five families at similar elevation range in the Andes, in Antioquia, Colombia, however, in the latter, Melastomataceae is also among the top families [133]. The family Asteraceae dominated at most elevations in congruence with studies in western México: the karstic Cerro Grande massif, Sierra de Manantlán [56], the volcanic and plutonic Sierra de Cacoma, Talpa de Allende, Jalisco, from 1600 to 2200 m a.s.l., with its greatest contribution being from 1750-2000 m a.s.l. [61]. Geographic isolation and optimal growth conditions led to the evolution of a particularly large number of Asteraceae in the montane regions of Guatemala and Mexico [134]. High species richness of Asteraceae along most of the Juquila gradient could be attributed to the high dispersal efficiency of the family through their anemochorous and ectozoochorous mechanisms, but also to their antiherbivore defenses through the production of secondary metabolites [135]. However, this family showed the lowest values at the middle portion of the CF gradient (1800-1899 m a.s.l.), this may be explained by the fact that Asteraceae do better in dryer habitats and this elevation range is associated to the highest moisture regime observed between 1835 and 1944 m a.s.l., possibly corresponding to a unimodal distribution of precipitation and or an area of high condensation of moisture from trade winds [136] out of the Gulf of Mexico, as evidenced by the 102% gain from fog through water collectors during the dry season, peaking at 1898 m a.s.l. in Altotonga, Veracruz, also in the windward side of the Gulf of México, in the neighboring Sierra Madre Oriental [137]. Other CF with low diversity of Asteraceae include the mesic CF in the Cerro Grande Sierra de Manantlán, in western Mexico, where Orchidaceae and Aspleniaceae out-number the Asteraceae only from 2300 to 2400 m a.s.l. [56]. The observed decrease of Asteraceae at 1800−1900 m a.s.l. in the Juquila gradient could be explained by high moisture. For instance, along an elevational gradient in Cerro Grande, Manantlán [56], Asteraceae shows its lowest richness value at 2300 m a.s.l. in moist CF with higher fog incidence, where umbrophilic ferns (Aspleniaceae), orchids, and Lamiaceae outcompete heliophytes like Asteraceae due to greater efficiency in using understory diffuse light [56,138]. In the same Manantlán gradient, at 2400 m, again, orchids out-number Asteraceae as the leading family. Similarly, but in a contrasting ecosystem, wet grasslands in Uganda have less diversity of Asteraceae than dry grasslands [139]. Higher moisture or rainfall could also make pappus and bristles wet, interfering with wind-dispersed species, except for those with ectozoochory.

β-Diversity
Gradient length estimates, except for guild data (0.39), match the expected outcome, ranging from 1.36 to 11.98 for field studies [71,72], and the length observed on structural vs. P/A data tends to be higher, representing their high structural heterogeneity, complexity, and turnover. Comparison of distance estimates between datasets should, however, be restricted to situations with the same abundance scale, range of the abundance scale, and species deletion level [71].

γ-Diversity
Gamma diversity values are not directly comparable among CF studies, due to different sampling schemes, elevational range, and area size covered, among many other factors. For instance, a 1300 m elevational CF gradient in Cofre de Perote, eastern Mexico, had 128 tree species ≥5 cm dbh [58], while for Juquila Vijanos, after excluding species with dbh <5 cm and ≥2.5 cm, we reported 64 woody species ≥5 cm dbh. Thus, gamma species richness of the Juquila gradient is half of that of the more northern Cofre de Perote (Veracruz), however, the altitudinal range of the Juquila gradient is only 26% of that of Cofre de Perote. This suggests that the gamma diversity in the southern Juquila gradient is proportionally greater than that in the Cofre de Perote gradient, in agreement with the globally known pattern of increasing species richness with the decrease in latitude [140].

Floristic Affinities
Floristic affinities of these CF, in terms of species composition, are similar to that of other CF in the Sierra Madre de Oaxaca and the Chimalapas, including La Chinantla, Sierra Mazateca, and Cerro Salomón, Chimalapas [6]. It is worth noting that the generic composition of the Juquila gradient is highly similar (ca. 50%) to that of the Miocene flora of northern Chiapas [64]. This Sierra de Juárez gradient also has affinities to the Cretaceous formations: Tarahumara (Sonora) and Cerro del Pueblo (Coahuila), which include Alnus, Betulaceae; Hedyosmum, Chloranthaceae; Quercus, Fagaceae; Juglandaceae (possibly), Magnolia, Magnoliaceae; Myrtaceae; Pinus, Pinaceae [141], while floristic affinities with Early Eocene floras of the formations La Carroza, La Trinidad (Chiapas), and Jackson (Tamaulipas) include Alnus; Ilex, Aquifoliaceae; Liquidambar, Altingiaceae; Oreomunnea, Juglandaceae; Myrica, Myricaceae; Eugenia, Myrtaceae, and the families Asteraceae, Caryophyllaceae, Fabaceae, and Poaceae [139]. The presence of additional boreal elements such as Pinus, Quercus, and Ilex represent the effect of decreasing temperatures before the Pleistocene, which allowed colonization to southern latitudes, as evidenced in Veracruz (Paraje Solo) of Mid-Pliocene [141,142].

Endemism and Conservation
Temperate semi-humid zones have higher endemism at the species level than warm humid zones [143]. The high endemism of these CF is explained by various factors: (a) the long permanence of the topographic relief of the Sierra de Juárez since the Miocene ca. 16 Ma, (b) the constant volcanic activity in eastern Oaxaca during the Miocene [144] contributed to the modification of environments, climates, and niches, as well as the formation of geographical barriers, inducing vicariance and allopatric speciation, (c) the effective ecological isolation of these CF archipelagos [8] that have resulted from the decrease of the space matrix with the increase in elevation, and (d) the high physiographic, geological, and edaphic heterogeneity [143]. However, some taxa managed to overcome barriers, and some constitute major disjunctions: Eastern Mexican and eastern North American disjunction, separated by northern Mexico and the southern Texas gap [145], Transisthmic disjunctions, separated by the Tehuantepec Isthmus, and northern vs. southern Mesoamerican disjunctions, separated by de Nicaraguan depression [146].
The proportion of Critically Endangered species in Juquila, Oaxaca (4%), were similar to that reported for the ecological reserve El Cielo, in Tamaulipas (4%) [67], and in Cofre de Perote National Park, Veracruz (3%) [58] and Sierra de Cacoma, Jalisco (3%) [57,61]. Endangered species in Juquila, Oaxaca (9%), represented a lower fraction than in Cofre de Perote, Veracruz (16%) and El Cielo, Tamaulipas (11%), and it was similar to Sierra de Cacoma, Jalisco (8%). Near Threatened species had a higher proportion in Juquila, Oaxaca (12%), than in Sierra de Cacoma, Jalisco (1%), Cofre de Perote, Veracruz (9%), and El Cielo, Tamaulipas (7%). The proportion of Vulnerable species in Juquila, Oaxaca (18%), was lower than in Cofre de Perote, Veracruz (21%), and higher than in El Cielo, Tamaulipas (11%) and Sierra de Cacoma, Jalisco (9%). In Juquila Vijanos, the high proportion of species under a risk category (31%) [20] suggests that this area may have been functioning has a biodiversity refuge [14,64]. This result confirms that the Sierra Madre de Oaxaca is one of the major biodiversity hotspots in Mexico [147] and the area with the largest number of threatened taxa reported for the Sierra Norte of Oaxaca (Ixtlán, Villa Alta y Mixe districts) [148]. However, the conservation status of many tree species is poorly documented, and many species are yet to be included in the Red List of Mexican cloud forest trees [20].
A considerable number of floristic elements along this gradient indicate disturbance, for instance, the dominant shade-intolerant woody Asteraceae and several other gap-phase regeneration species such as Hedyosmum mexicanum, Phyllonoma laticuspis, Magnolia dealbata, Clethra mexicana, Liquidambar styraciflua, Palicourea padifolia, Pinus chiapensis, Miconia glaberrima, and Vaccinium leucanthum. This disturbance aspect is acknowledged as a limitation in this research.

Conclusions
Major small-scale CF community variation was explained by soil pH, while elevation, K, Na, and firewood extraction were environmental variables of secondary importance. This supports the habitat specialization hypothesis as the major driver of small-scale community organization of CF stands harboring M. dealbata in Juquila, Oaxaca.
The high discontinuity separating the Oreomunea-Ticodendron forest stands do not support the Continuum hypothesis. This relic Oreomunnea-Ticodendron-CF resulted isolated at the high acidity end of the gradient, and there, species richness reached its maximum, where Magnolia dealbata had its lowest density and basal area. In contrast, Magnolia dealbata presented the highest density values in the Vismia-CF and the highest basal area in the Alnus-Ocotea-CF and Siparuna-CF.
The increasing species richness with increasing elevation does not support the Decreasing hypothesis, but this may be due to a sampling artifact, sampling truncation at the top of an assumed unimodal (mid peak) pattern, or may be explained in part by a combination of habitat specialization, contributions from biogeographic refuges, and from sites with good conservation status. The high species richness observed was associated with the sea windward dry-seasonal maximum cloud condensation gain and habitat differentiation-specialization. We hypothesize that a combination of these factors may also explain the hypothesis of mid-peak species richness and may have allowed long-term survival of relict species such as Oreomunea Mexicana and Ticodendron incognitum. Future studies may consider testing the unimodal (mid-peak) hypothesis for species-rich and mist-loving vascular families such as Aspleniaceae, Grammittidaceae, Orchidaceae, and Bromeliaceae, on longer altitudinal transects, a wider range of edaphic and topographic variables, using different scales, and incorporating land-use change, among additional factors.
The species-richest Oreomunnea-Ticodendrum-CF, considered an Oligocene-Miocene relic, must be legally protected urgently, and since it is isolated and specialized to the acidic end of the gradient, it requires specific maintenance of soil conditions. The generated or tested hypotheses about major environmental variables influencing community organization as well as the relationships among various CF groups should be used to guide conservation and sustainable development strategies for these communities harboring the bioculturally relevant Magnolia dealbata.