Study Site

El Cielo Biosphere Reserve (CBR) comprises 144,000 ha in northern Mexico (Figure 1). It is located in the Sierra Madre Oriental, a cordillera with great ecological and historical complexity (Halffter, 1962). Human activities, such as deforestation and agriculture, and natural events, such as multiple treefall, caused serious disturbances in the CBR before 1984, when the CBR was declared as a protected natural area (Arriaga, 2000). Currently, low-impact agriculture is still practiced in areas adjacent to semideciduous forest (SDF), causing some disturbance (Vargas-Vázquez et al., 2019). Most cloud forest in Mexico is found in the Sierra Madre Oriental, and this region, specifically the CBR (Ponce-Reyes et al., 2012; Rojas-Soto, Sosa, & Ornelas, 2012), represents the northernmost limit of cloud forest (Halffter, 1962; Ponce-Reyes et al., 2012).

Figure 1.

Distribution of different forest types in the “El Cielo” Biosphere Reserve in Tamaulipas, Mexico. The semideciduous forest (300–400 m a.s.l., triangle) and the tropical montane cloud forest (800–1600 m a.s.l., square) are distributed along the windward slope; and the pine-oak forest (1,500–2,000 m a.s.l., circle) and the submontane scrub (700–1,500 m a.s.l., diamond) are distributed along the leeward slope of the reserve.

In this reserve, there is a mixture of tropical and temperate species. Temperate species are dominant overall, but tropical species are dominant at lower elevations (Arriaga, 2000). Climatically, this area is considered to be temperate (warm, temperate with a hot summer climate; Beck et al., 2018; Peer et al., 2007). At a local scale, the topography varies, resulting in a wide diversity of climates: The windward slope of the Sierra Madre Oriental is wet and influenced by marine winds, whereas the leeward slope has a dry climate (Sánchez-Santillán, Binnqüist, & Garduño, 2018) and vegetation types.

We sampled four vegetation types along a climatic gradient from 300 to 2,000 m a.s.l. (Table 1 and Figure 1; de la Rosa-Manzano et al., 2017): SDF, tropical montane cloud forest (TMCF), submontane scrub (SMS), and pine-oak forest (POF). The semideciduous forest (SDF), (300–400 m a.s.l.) is characterized by a warm, subhumid climate and is dominated by Bursera simaruba (L.) Sarg., Casimiroa greggii (S. Watson) F. Chiang, Guazuma ulmifolia Lam., Lysiloma divaricatum (Jacq.) J.F. Macbr., Aphanthe monoica (Hemsl.) J. F. Leroy Sharp, Cascabela thevetia (L.) Lippold., and Trema micrantha (L.) Blume, among others, with heights ranging from 8 to 15 m. The TMCF (800–1,600 m a.s.l.) is mainly characterized by a warm climate, and common tree species include Magnolia tamaulipana A. Vázquez, Quercus germana Schltdl. & Cham., Q. sartorii Liebm., Podocarpus matudae Lundell, Acer skutchii Rehder, Cercis canadiensis L., Carpinus caroliniana (Marshall) Fernald, Ostrya virginiana (Mill.) K. Koch, with heights between 18 and 25 m. These forests (SDF and TMCF) are distributed along the windward slope of the mountain range (González-Medrano, 2005).

Table 1.

Summarized Data of the Vegetation Types Along the Elevational Gradient in the “El Cielo” Biosphere Reserve.

The SMS (700–1,500 m a.s.l.) is characterized by a dry climate and species such as Acacia berlandieri Benth, A. rigidula Benth, Acanthocereus tetragonus (L.) Hummelinck, Helietta parviflora (A. Gray ex Hemsl.) Benth., Cordia boissieri A. D.C., Havardia pallens (Benth.) Briton & Rose, Gochnatia hypoleuca (SC.) A. Gray, and Yucca treculeana Carrière, among others, which reach up to 4 m in height. And, the last forest type, POF (1,500–2,000 m a.s.l.), is characterized by a dry and cool climate and dominated by Pinus teocote Schltdl. & Cham., P. pseudoestrobus Brongn., C. caroliniana Walter, Carya myristiciformis (F. Michx.) Nutt., Carya ovata (Mill.) K. Koch, Q. germana Schltdl. and Cham., Q. affinis Scheidw., and Q. sartorii Liebm., which reach up to 25 m in height (González-Medrano, 2005). These forests (SMS and POF) are distributed along the leeward slope of the mountain range (González-Medrano, 2005).

We sampled a total of 32 plots, each 400 m², including 5 plots in SDF, 15 in TMCF, and 6 in both SMS and POF between 2017 and 2018. The number of plots per vegetation type was determined according to the cover of each vegetation type. All plots were selected in areas with a relatively high density of epiphytes (e.g., old growth vegetation). In each plot, host trees with a diameter at breast height (DBH) of more than 1 cm were sampled; hosts other than trees were also sampled (e.g., agaves and cacti, particularly in the SMS). Host traits such as height, measured using a distance meter (D210, Leica); DBH; and species identity were also recorded. Epiphytic seedlings were not surveyed because they are difficult to identify and have high mortality rates (Winkler, Hülber, & Hietz, 2005). When epiphytic individuals were infertile, we identified them by comparing them to nearby fertile individuals of similar appearance. Epiphytic individuals were distinguished as epiphytic stands (sensu Sanford, 1968; i.e., groups of spatially aggregated stems), and vegetative clones from rhizomatous plants or clumped individuals were considered as a single individual. Sampling using binoculars (EO-D102; Eagle optics) was carried out in the SDF and SMS and in the TMCF and POF using the single-rope climbing technique (Perry, 1978). It is possible that some species were missed in the SDF even with the use of binoculars, as the species accumulation curve did not reach the asymptote (Figure 2 and Online Appendix 2; de la Rosa-Manzano et al., 2017). Nonetheless, the use of binoculars to survey epiphytes is a common practice (Leimbeck & Balslev, 2001). Epiphytic species were divided into six functional groups considering the general ecological characteristics of their families: Bromeliaceae were divided into atmospheric and tank forms, and Orchidaceae into those with or without pseudobulbs. Cactaceae were classified as succulent, whereas Piperaceae were left as such (Table S2).

Figure 2.

Inter- and extrapolation analysis (iNEXT) of epiphyte diversity according to forests type: semideciduous forest (300–400 m a.s.l.), tropical montane cloud forest (800–1,600 m a.s.l.), submontane scrub (700–1,500 m a.s.l.), and the pine-oak forest (1500–2,000 m a.s.l.; Table 1). The interpolated (observed) species richness of vascular epiphytes was obtained by merging abundance and species richness of plots per forest type. Extrapolated (expected) species richness for each forest type was based on the highest maximum number of individuals (ca. 7,500). SDF = semideciduous forest; TMCF = tropical montane cloud forest; SMS = submontane scrub; POF = pine-oak forest.

Voucher specimens of epiphytes and hosts were deposited in the Francisco González Medrano herbarium of the Applied Ecology Institute of the Autonomous University of Tamaulipas.

Data Analyses

Variables determining vascular epiphytic species richness, abundance, and composition

To assess whether differences in the number of sampled plots affected sample completeness, we calculated the proportion of the species pool found in each forest type at the plot level by means of interpolation and extrapolation analyses (Chao et al., 2014), which evaluate sample completeness based on sample coverage. We used the iNEXT function in the iNEXT package (Hsieh, Ma, & Chao, 2016).

To explore how the 3 structural variables (height, host DBH, and surface area), 14 climatic variables (Table 2), and elevation determined epiphytic species richness and abundance, we used linear and generalized linear models, respectively. We ran a total of four models: a structural and climatic model at the tree and plot level for each type of response variable (species richness and abundance). Then, we compared the explanatory power of the structural and climatic models according to the response variables.

Table 2.

Output of the Linear and Generalized Linear Models Assessing the Effect of Structural (DBH, Height, and Surface of Host Tree) and Climatic Variables (MTWeQ, PCQ, and APsq) on Species Richness and Abundance.

At the tree level, the structural linear models included height, DBH, and surface area of host trees. Host tree surface area was calculated from the DBH and tree heights using an equation based on that described in Kershaw, Ducey, Beers, and Husch (2016). Elevation rather than forest type was included as a fixed factor to account for variation in forest structure among forest types, to economize the degrees of freedom, and to avoid multicollinearity. To achieve normality, we log-transformed the response variables (Crawley, 2005). To check whether the assumption of homoscedasticity was fulfilled (function bptest), we used the Breusch-Pagan test (Breusch & Pagan, 1979).

At the plot level, we fitted the generalized linear models with Poisson distributions. These climatic models included climatic variables and elevation as response variables. Climatic variables were retrieved from CHELSA version 1.2, which provides climatologies at high resolution for the Earth’s land surface, and a high-resolution climatic dataset (Karger et al., 2017). We retrieved a total of 19 bioclimatic variables but only used 14; redundant variables (e.g., medium diurnal range, which is defined as the difference between the daily maximum and minimum temperature) and those with a correlation equal to or above 90% (Pearson correlation) were removed to avoid multicollinearity. The utilized climatic variables were annual mean temperature (AMT), maximum temperature of the warmest month, minimum temperature of the coldest month, mean temperature of the wettest quarter (MTWeQ), mean temperature of the driest quarter, mean temperature of the coldest quarter, mean temperature of the warmest quarter, annual precipitation (AP), precipitation of the wettest quarter, precipitation of the driest quarter, precipitation of the warmest quarter, and precipitation of the coldest quarter. We also included the quadratic terms of annual mean temperature (AMTsq) and annual precipitation (APsq) because the relationship between these variables and species richness may be nonlinear (e.g., Tang, Li, Li, & Meng, 2014). In each model, the dependent variables were standardized to assess the importance of their effects.

For the structural models, the best model was selected using backward stepwise selection based on the χ² test (function drop1). For the climatic models, the best model was selected using forward stepwise selection (function add1) based on the AIC values. We first ran a full model (where all explanatory variables were included) and a null model (without explanatory variables). Then, we added one variable at a time, starting with those with the smallest AIC value, to the null model. The best model was determined based on the variance inflation factor (function vif) calculated for each updated model (terms with values above 5 were excluded; Online Appendix 6). The significance of each variable in the final models was calculated using the analysis of variance function in the car package.

General Results

We registered a total of 12,103 epiphytic individuals belonging to 4 families, 15 genera, and 30 species. The highest species richness was found along the windward slope in the TMCF (16 species, or 33% of the total species richness on 670 trees), followed by the SDF (15 species, or 31% of the total species richness on 240 trees; Table 1 and Figure 2). Only about half this number of species, or 8 and 9 species, was found along the leeward slope in the POF and the SMS, respectively (17% on 215 trees and 19% on 384 trees, respectively; Table 1 and Figure 2). Overall, TMCF was significantly more diverse than POF and SMS. However, because the diversity profiles crossed at several points (q = 0–2), no single forest type can be considered more or less diverse than the others, that is, they had comparable diversity. The SDF was the most heterogeneous or, in other words, had the highest number of rare species and lowest number of dominant species (Online Appendix 1). Ultimately, the species richness among forest types was comparable despite differences in sampling effort (Figure 2 and Online Appendix 2). Coverage values were very high, and extrapolated species richness indicated that only a single additional species could be found in SDF and SMS by increasing the sample size (number of individuals per plot; Figure 2 and Online Appendix 2). Finally, species richness slightly decreased with elevation, although this trend differed according to family. For example, orchid species richness showed a steep decline with elevation, while bromeliad and other species richness showed no trend (Online Appendix 3). Only the epiphytes Tillandsia pringlei and Tillandsia usneoides were shared by all vegetation types.

Interestingly, there was no correspondence between epiphytic density (ind/cm²) and host tree species richness or between host tree abundance and surface area per forest type (Table 1). Species richness and abundance (per tree) generally increased with host tree surface area, although the relationship varied according to forest type and response variable (Online Appendix 4). Species richness steadily increased with host tree surface area, with the exception of TMCF, which showed large variability in tree size (Table 1). The relationship between abundance and host tree surface area was highly dependent on forest type, especially in the POF, where abundance largely varied (Online Appendix 4). At the plot level, the number of epiphytic species per plot was positively correlated with the number of host tree species per plot (Online Appendix 5).

Climatic and Structural Variables Determining Species Richness and Abundance

Both structural and climatic variables significantly determined species richness and abundance along the climatic gradient, but their effects varied according to response variable and scale (Online Appendix 6). For total species richness, the precipitation of the coldest quarter and MTWeQ were the most important climatic variables, having a positive effect on species richness in comparison to the structural variables (DBH and height), which also had a significant positive effect but, unlike climatic variables, only explained half of the variation in species richness (Table 2 and Online Appendix 6). For total abundance, the converse pattern was found: The structural variables largely determined abundance, and only one climatic variable had a significant negative effect, explaining a very low percentage of variability. The DBH had the strongest positive effect on abundance and was more than double the effect of height (Table 2 and Online Appendix 6). These findings suggest that climatic variables have by far the largest effect on species richness, while structural variables mainly determine abundance.

Climatic and Structural Variables Determining Dissimilarity in VEAs Among Forests

All dissimilarity analyses showed a very clear differentiation of VEAs according to forest type in terms of both taxonomic and functional composition (plot-based comparisons; Figure 3 and Table 3) as well as a clear difference in composition between slopes (Figure 3, Table 3, and Online Appendices 8 and 9). For instance, bromeliad indicator species clearly showed differences between slopes. In the forests (POF and SMS) along the leeward slope, the atmospheric species were found to be the indicator species. In the forests (SDF and TMCF) along the windward slope, mostly tank species were indicator species, although we also found a high diversity of indicator species in the SDF (six species belonging to five functional groups; Figure 3 and Table 3).

Figure 3.

Nonmetric multidimensional scaling ordination showing dissimilarity of VEAs based on relative abundance (Chao index) data between forest types (at the plot level). Species shown bigger and in bold are indicator species which are categorized according to functional group (color coded, upper-right corner legend, Table 3). Polygons reflect the plot coordinates for each forest type and the border line type indicates the forest type (upper-left corner legend). Gray filled polygons indicate forest types located along the windward slope and white polygons indicate forest types located along the leeward slope. Total dissimilarity percentage among plots (within-forest type) is indicated above each forest type polygon. NMDS = nonmetric multidimensional scaling.

Table 3.

Indicator Species Analysis Values.

Within-forest dissimilarity (among plots) was highest in TMCF and lowest in SDF (windward slope), whereas in the forests along the leeward slope (POF and SMS), dissimilarity was comparable (Figure 3 and Online Appendix 8). Moreover, dissimilarity between the forests on the leeward slope was lower (33%) than that of those on the windward slope (54%; Figure 3 and Online Appendix 8). Dissimilarity among the four forest types was ca. 70% and almost entirely due to species turnover. In other words, dissimilarity was due to actual differences in species composition (species turnover) and not merely to species loss from site to site (Baselga & Orme, 2012), highlighting the specificity in the epiphytic composition of each forest type (Figure 4 and Online Appendix 9).

Figure 4.

CCA assessing the variation of VEAs according to climatic variables (at the plot level). Climatic variables comprise the mean temperature of the wettest quarter and annual precipitation. The proportioned explained by constraining variables was of 30.58%. Gray filled polygons indicate forest types located along the windward slope and white polygons indicate forest types located along the leeward slope. CCA = constrained correspondence analysis.

While structural variables explained 5% of the total variability of VEAs (at the tree level; Online Appendix 7), climatic variables comparatively explained the variability by more than six fold (30.58% at the plot level; Figure 4), highlighting the differential effects of the drivers of the variability of VEAs at different scales. In the climatic model, the first axis explained over 76% of the variation and was related to AP; in contrast, the second axis explained 24% and was related to the MTWeQ (Figure 4). In the structural model, the first axis explained 68% of the variation and was related to variation in height; meanwhile, the second axis explained 32% and was related to changes in DBH. High structural variation was observed in TMCF and SDF (both along the windward slope), as observed along the second axis gradient (Online Appendix 7).

Epiphytic Richness in a Transition Zone

Vascular epiphytic richness in the CBR was low (30 species in 1.28 ha) in comparison to that of tropical rainforests (e.g., Gentry & Dodson, 1987; Kreft et al., 2004). Sampling four vegetation types resulted in double the number of species compared with de la Rosa-Manzano et al. (2017), who only sampled two vegetation types. Overall, the low species richness encountered herein is comparable to that of some subtropical forests (e.g., Atlantic and broadleaf temperate forests; Barbosa, Becker, Cunha, Droste, & Schmitt, 2015; Ceballos, Chacoff, & Malizia, 2016; Díaz et al., 2010; Hsu et al., 2012; Xu & Liu, 2005) but lower than that of other subtropical forests in the southern hemisphere (e.g., Brown, 1990; Hofstede, Dickinson, & Mark, 2001). This low richness may possibly be due to the decreasing number of species per unit area toward the poles, as observed in coniferous forests (e.g., Khine, 2018). In epiphytes, the drop in the number of species seems distinctively sharp (Moran, 2008; Karger, Kluge, & Kessler, 2016; Salazar et al., 2013).

High epiphytic diversity in tropical forests is generally thought to result from the positive effect of precipitation (Kreft et al., 2004; Wolf & Alejandro, 2003). Also, a general pattern of high richness is observed at middle elevations (Gentry & Dodson, 1987), which can be attributed to climatic conditions such as cloud formation, temperature, and precipitation (Carvajal-Hernández & Krömer, 2015; Gentry & Dodson, 1987; Hemp, 2001; Kessler, 2000; Kessler, Herzog, Fjeldsa, & Bach, 2001; Kluge, Bach, & Kessler, 2008; Kluge, Kessler, & Dunn, 2006; Stevens, 1992; Wolf, 2005). In contrast to this general trend found throughout the tropics, total species richness in the CBR decreased with elevation. However, this trend was largely due to the behavior of Orchidaceae in the transition area. Although the study site does not comprise an elevational gradient in the strict sense, the lack of a middle elevational peak in the CBR could be related with the narrow elevation range (the highest elevation reaches 2,200 m a.s.l.; González-Medrano, 2005) or the strong seasonality of the area (García-García & Zarraluqui, 2008, Salazar et al., 2013).

The overall diversity of the VEAs only differed between the TMCF, the POF, and the SMS; the TMCF was more diverse than the POF and the SMS, yet no single forest type was statistically more diverse than all others, meaning that they had comparable diversity. The SDF and the TMCF along the windward slope, which corresponded with warm, subhumid, and humid climates, had higher species richness in contrast to the POF and the SMS along the leeward slope, which corresponded with temperate, subhumid, and dry climates (Sánchez-Santillán et al., 2018).

Climatic and Structural Factors as Determinants of Epiphytic Diversity

Species richness was largely determined by climatic variables in the study area, while abundance was mainly determined by forest structure. Climatic variables associated with seasonality were important determinants of species richness, while abundance was determined to a lesser degree by APsq (quadratic term). Although precipitation is known to be the main determinant of epiphytic species richness in tropical forests (Kreft et al., 2004), in this transition area, seasonality seems to be crucial for vascular epiphytic richness. The very complex topography of the CBR with two opposing slopes (windward and leeward) in conjunction with pervasive summer droughts, which are especially intense in dry areas (Sánchez-Santillán et al., 2018), and strong seasonality of precipitation (Karmalkar, Bradley, & Diaz, 2011) may explain why climatic variables related to seasonality, and not simply AP, determine vascular epiphytic richness in this transition area.

Species abundance was largely determined by structural variables, that is, DBH and height, which were also important for species richness, although to a lesser degree (Table 2). The determination of epiphytic distribution by forest structure is one of the best-known spatial patterns in vascular epiphytes (Zotz, 2007) that exists independently of forest type (e.g., Flores-Palacios & García-Franco, 2006; Hirata, Kamijo, & Saito, 2009; Laube & Zotz, 2006). This seems to be an intrinsic characteristic of vascular epiphytes given their structural dependence on a dynamic substrate, that is, host trees (Spruch, Hellwig, Zotz, & Blasius, 2019). These results suggest that, in the transition zone of the CBR, climatic seasonality at a regional scale determines species richness and at the tree level (host tree size) determines the abundance of VEAs.

The peak in orchid species richness on the windward slope (wet side) reflects the more balanced climatic conditions in contrast to the leeward slope, where extremely low temperatures occur in the POF and severe droughts in the SMS. In fact, in the latter forest type, we found a single orchid species. The TMCF on the windward slope registered the highest rainfall in the area (Table 1) and the frequent presence of fog (Hamilton et al., 1995). This fog or “horizontal precipitation” (Vogelmann, 1973) is likely to offset the short period of drought that occurs in the CBR (3–4 months; Rzedowski, 1996). Fog is known to contribute up to 20% or more of the total water input in cloud forests (Juvik & Ekern, 1978; Stadtmüller, 1987), and it may facilitate the presence of species that require wetter environments, such as orchids. On the other hand, very low temperatures at high elevations limit the distribution of orchid richness (Zotz, 2005). For instance, in the POF, only three species were registered at 1,800 m a.s.l., similar to another pine forest in Veracruz, Mexico, where only a single species of orchid was found (Hietz & Hietz-Seifert, 1995). Meanwhile, the high orchid richness in a pine forest in Chiapas in southern Mexico (Wolf & Alejandro, 2003), located seven degrees south of the study site, highlights the strong role of temperature as the driving force of orchid richness. The high richness in this latter site may result from the beneficial effects of temperature and high humidity, which likely bypassed the negative effects of the phenolic or resinous composition of pine bark that hamper orchid establishment (Wolf & Alejandro, 2003).

On the leeward slope, atmospheric bromeliads such as Tillandsia recurvata and T. usneoides prevailed, similar to what is reported by Rzedowski (1996) for the arid regions of Mexico. In this case, mountains may prevent the passage of moist currents from the coast, thus interrupting the flux of moisture beyond the high areas and causing semi-arid regions, like SMS, where fog banks create favorable environments for epiphytes. The present findings confirm our hypothesis that climatic variables have an important effect on total epiphytic species richness and abundance.

VEAs Are Differentiated According to Forest Type

The variation in VEA composition among forests was related to gradients of structural (DBH and tree height) and climatic variables (temperature and precipitation). The strong structural and climatic differences among forests were clearly reflected in the compositional differences of the VEAs. High dissimilarity was mostly due to turnover, especially in the forests along the leeward slope (dry and cold climate). The exception was the TMCF, where within-forest dissimilarity was very high (higher than among-forest dissimilarity), largely due to nestedness. Such nestedness patterns point out the relevance of dispersal limitations in structuring these assemblages (Baselga & Araujo, 2009), which is highly relevant considering that the TMCF harbors endemic species or those with restricted distributions. Within the small area of the CBR (144,000 ha), the four forest types shared only 2 out of 30 species, even though these species have a very wide latitudinal distribution (GBIF, 2019). These strong dissimilarities in the CBR indicate that the regional diversity is highly partitioned among forest types and trees within-forest types.

The degree of differentiation of VEAs among forest types in the CBR is surprisingly high considering that this site experienced a period of disturbances up until approximately 35 years ago (Arriaga, 2005). In a previous study, such differentiation was not evidenced in epiphytic orchid composition between pastures, regenerating forest, and old growth forest across six reserves, probably due to favorable climatic conditions (Sanchez, Armenteras, & Retana, 2016). In contrast, in another study carried out in a single site of TMCF, there were clear differences in epiphytic composition between isolated remnant trees and forests with a similar elevation that had not experienced disturbances for 50 years, although epiphytes were sampled in the same host species (Larrea & Werner, 2010). While these examples involve the comparison of disturbed sites to old growth forest, where differences in composition are expected, the contradictory results hint at the high complexity of the determinants of epiphytic diversity. Finally, in a lowland forest of the Colombian Amazonia, a similar differentiation of VEAs among landscape units was found, although this was largely explained by differences in host species composition despite structural similarity among landscape units (Benavides, Vasco, Duque, & Duivenvoorden, 2011). Although host species richness was associated with epiphytic richness in the CBR (Wagner et al., 2015), discriminating between the effect of host tree composition and structural variability is difficult and requires a more detailed analytical framework.

Notably, Bromeliaceae and Orchidaceae were specialized to forest type, supporting our prediction. In this study, bromeliad species comprised both tank and atmospheric forms. The former has broad and flat leaves that basally overlap, forming water-tight chambers that store water (Adams & Martin, 1986). The latter has very narrow leaves and abundant foliar trichomes specialized in water and nutrient absorption (Martin, McLeod, Eades, & Pitzer, 1985); these bromeliads obtain moisture and nutrient ions directly from the atmosphere (Benzing, 1990). Atmospheric bromeliads were characteristic of both forests along the leeward slope with a cold and dry climate and low precipitation. Here, drought-tolerant epiphytic species of the genus Tillandsia (atmospheric type) were more prevalent than in areas with high or moderate rainfall (Chilpa-Galván, Tamayo-Chim, Andrade, & Reyes-García, 2013; García-Suárez, Rico-Gray, & Serrano, 2003). In particular, T. usneoides was an indicator of POF and T. recurvata of SMS; both species show adaptations to cope with water scarcity and high temperatures (Benzing, 1990). On the other hand, tank bromeliads and both orchid groups (with and without pseudobulbs) were indicators of forests along the windward slope. In this case, the tank bromeliad Tillandsia viridiflora was an indicator of TMCF, which receives the highest precipitation in the CRB. It is known that tank bromeliads are usually restricted to humid forests (Pittendrigh, 1948), although in dry forests, some tank bromeliads are distributed at the top of their hosts to better access rainfall (Cach-Pérez et al., 2013; Graham & Andrade, 2004; Reyes-García, Griffiths, Rincón, & Huante, 2008).

Almost all indicator species of the orchid groups were found on the windward slope. Indicator species in the SDF have developed either large pseudobulbs (26 × 8 cm for Catasetum integerrimum; Hágsater & Salazar, 1990), deciduousness, or tilting succulent leaves (de la Rosa-Manzano, Andrade, García-Mendoza, Zotz, & Reyes-García, 2015). These adaptations can aid these orchids in storing water for long dry periods, tolerating high temperatures, avoiding water loss, and surviving in harsh environments. Functional group differentiation has been well studied in bromeliads (Benzing, 2000; Pittendrigh, 1948) but is still unclear in orchids; therefore, it is important to continue to study the inter- and intraspecific variability of orchid functional traits in heterogeneous environments.