Materials and methods

Study area: This work was carried out from May 2016 to April 2017. The study area is located northwest of the municipality of Victoria, in the state of Tamaulipas, Northeastern Mexico. Succession gradient was delimited on a small homogeneous plain of approximately 400 ha (4 km²), at 320-350 masl with East exposure. The area is set between the skirts of the Sierra Madre Oriental (SMO) and rural areas of Ejido (Ej.) Rancho Nuevo and Ej. Santa Ana (23°51’00” - 23°51’45” N & 99°13’15” - 99°14’00” W). Dominant vegetation type in the studied fragment is low thorn forest (^{INEGI, 2013}), although it shares elements with submontane scrub and Tamaulipan thorn scrub communities (^{García-Morales et al., 2014}) in areas of secondary vegetation. Low thorn forest in the study area can be designated as deciduous (INEGI, 2013) or semideciduous (^{Treviño-Carreón & Valiente-Banuet, 2005}). It is characterized by thorny tree species, 8 to 10 m height, and a high proportion of species (50 to 75 %) that loses their leaves during the dry season (Treviño-Carreón & Valiente-Banuet, 2005). Climate in the area is warm subhumid with summer rains. Annual average temperature ranges from 18 to 24.3 °C, and total annual precipitation fluctuates between 717.3 to 1 058.8 mm (^{Almaguer-Sierra, 2005}); however, long drought periods are a common pattern during several months.

Successional gradient: The succession gradient was delimited based on the processing, classification and analysis of Landsat satellite images of years 1973, 1986, 2000, 2005 and 2013, which were obtained from the Global Land Cover Facility (^{GLCF-GLS, 2016}) and USGS Global Visualization Viewer (^{GLOVIS, 2016}) websites. Each of the four last images had a spatial resolution of 30 m; thus, it was necessary to modify the cell size of the 1973 image, from 60 to 30 m. We followed a method for delimiting approximate time for chronosequences and successional time using satellite imagery (^{Sánchez-Reyes et al., 2017}), which is briefly described here. First, images were automatically segmented. Then, some segments were manually selected as training fields, from which a spectral signature was extracted to classify the images into each of the following four land cover categories: conserved low thorn forest, secondary low thorn forest, modified areas, and bare soil areas; such classification was done by using the maximum likelihood algorithm. Through a reclassification, the four land cover categories were merged as follows: vegetation (conserved and secondary low thorn forest), and disturbance (modified and bare soil areas). Using these two simplified categories, images were processed by cross tabulation to designate the final successional categories for this study. 1) Conserved areas. Closed canopy: trees with heights up to 15 m, diameters at breast height (DBH) up to 50 cm; an heterogeneous number of shrubs and herbaceous plants; species such as Celtis pallida Torr., Casimiroa greggi (S. Watson) F. Chiang, Ebenopsis ebano (Berland.) Barneby & J. W. Grimes, and Mascagnia macroptera (DC.) Nied. 2) Areas with 31 years of succession. Closed canopy: trees up to 9 m height, DBH up to 25 cm; a low number of herbaceous plants, but a contrastingly higher number of shrubs; species such as Havardia pallens (Benth.) Britton & Rose, Randia obcordata S. Watson, Cordia boissieri A. DC., and Croton cortesianus Kunth. 3) Areas with 17 years of succession. Open to closed canopy: trees up to 7 m height, DBH up to 18 cm; a high number of shrubs and herbaceous plants; species such as Acacia farnesiana (L.) Willd., H. pallens, Mimosa malacophylla A. Gray, and Sida acuta Burm. 4) Areas with four years of succession. Open areas: absence of trees, dominance of shrubs from 1 up to 2.5 m height, and an homogenous herbaceous cover on the forest floor; species such as Calyptocarpus vialis Less., C. cortesianus, Malvastrum coromandelianum (L.) Garcke, M. malacophylla, and Solanum erianthum D. Don. Characterization of each category was done on basis of identification and visual observation of vegetation. The study area historical disturbance is characterized by slash and burn for conversion to agriculture areas, which was corroborated by interviewing owners of local farms; therefore, original type of disturbance was the same for all areas. Procedures were conducted with the aid of IDRISI Selva 17.0.

Microclimate measurement: Eight plots of 10 x 10 m in each of the four successional categories were delimited randomly using Geographic Information Systems (GIS) software (Digital appendix 1). Successional categories consisted of several patches of different forms and areas; thus, for a random distribution of plots, all patches within each category were unified. Although we had several replicates of categories, the random procedure originated an unequal number of plots by patch (Digital appendix 1). Therefore, we considered plots as the replicates of successional categories, since through the used procedure these plots were randomly located. Besides, their independence was guaranteed by locating each plot with a minimum distance of 40 m from other plots or successional categories, in order to avoid edge effect and influence of surrounding microclimate (^{Davies-Colley, Payne, & Elswijk, 2000}; ^{Gehlhausen, Schwartz, & Augspurger, 2000}). Subsequently, microclimate measurements were performed monthly from May 2016 through April 2017 on each plot. Variables measured are as follows: 1) wind speed, 2) temperature, 3) relative humidity, 4) heat index, 5) dew point, 6) evapotranspiration. Heat index as considered here measures thermal stress or thermal sensation by combining temperature and humidity (^{Anderson, Bell, & Peng, 2013}; ^{Lee & Brenner, 2015}); as such, it has not been considered as a microclimatic variable for successional analysis, but recent studies show its association with species distribution (^{Sánchez-Reyes, Niño-Maldonado, Barrientos-Lozano, Jones, & Sandoval-Becerra, 2015}; Sandoval-Becerra et al., 2017). Microclimatic variables were recorded at the center of each plot using a Kestrel 4500 portable weather station, which was placed 1.5 m above the ground, attached to a metal pole, bearing a protection at the top to avoid direct solar radiation. The device was programmed to record variable measurements every 30 min for 12 h (06:00-18:00 h). Three weather stations were placed simultaneously, each one in a different plot. After recording measurements, the devices were transferred to the next three plots, until the 32 plots were evaluated. Therefore, variable measurements required 11 consecutive days monthly, that is, each plot had only one measurement by month. All measurements were conducted in clear or partly cloudy days, with no precipitation.

To assess the influence of seasonality on microclimate variation along the successional gradient, variables measured were analyzed taking into account two major seasons that characterize the study area (^{Secretaría de Gobierno, 2015}). Wet season, from May through October 2016, and dry season from November-2016 through April-2017; i.e., six months each season. For each microclimate variable, we calculated the average monthly value at each plot (the average value of measurements from 12 hours). These values were used to estimate statistically significant differences between succession categories, applying a discriminant function analysis (FD). Data were analyzed separately for the wet and the dry seasons. This method aims to find those factors that contribute most to separate the groups, by searching for the best linear combination of the original variables, which originates new variables known as discriminant functions or roots. Interpretation was accomplished taking into account the first two significant roots of the analysis, which represent the best discrimination between groups (^{Rencher, 2002}). The Wilk´s Lambda value obtained, is a measure of the deviations within each group with respect to the total deviations. It assumes values between 0 and 1. Values close to 0 indicate a high discriminant capacity of the variables; while values close to 1, indicate a low capacity to discriminate (^{Pozo-Díaz & Carrasco, 2005}).

The FD analysis allowed to determine the F and p values between pairwise comparisons of categories, to test for significant differences in succession categories based on microclimatic variables. In addition, the relative contribution of each variable measured to the separation of succession categories was quantified, according to the correlation values between these variables and the new discriminant functions; these correlations are known as loadings or factor structure coefficients (^{Rencher, 2002}).

Finally, the FD showed the optimal graphical distribution for each plot, i.e., it depicts the best discrimination of the groups/succession categories, based on its total microclimatic variation. Such observations are called canonical observations or canonical scores (^{James & McCulloch, 1990}; ^{Rencher, 2002}; ^{Karels, Bryant, & Hik, 2004}). Also, the average canonical values of the observations in each group, or centroids, were obtained; these represented the average microclimatic conditions in each succession category. The interpretation of both graphs (scores and centroids) was made according to the variable with the highest correlation value (factor structure coefficients or loadings) for each root, since they had the greatest contribution to the separation between succession categories. Discriminant function analysis was carried out using STATISTICA 8.0 (^{StatSoft Inc., 2007}); the graphs were generated in ArcGIS 10.9 (^{Esri, 2013}).

Results

The discriminant function analysis indicated significant global differences between the succession categories as a function of the microclimate, both in the wet season (F_{18, 518} = 2.7877; P < 0.0001) and the dry season (F_{18, 518} = 7.1442; P < 0.0000). The discriminant value of the model was higher in the dry season (Wilks ‘Lambda = 0.53414) than in the wet season (Wilks’ Lambda = 0.76991).

In the wet season, both roots of FD explained similar percentages of the microclimate variation between categories, with a 51.98 % (eigenvalue = 0.14422) for the root 1, and a 42 % (eigenvalue = 0.116534) in the root 2. According to pairwise comparisons, the microclimate observed in the category of 17 years does not differ from that registered in areas with four and 31 years of succession; the remaining comparisons were significantly different (Digital appendix 2). Microclimatic variation explained by the two FD roots was higher in the dry season, since 69.70 % (eigenvalue = 0.516675) was described by root 1, and only 22.21 % (eigenvalue = 0.164619) was explained by root 2. The pairwise comparisons during the dry season indicated that the microclimate is different amongst all successional categories (Digital appendix 2).

During the wet season, in root 1, the variable with the highest correlation value was the heat index (-0.60); while in root 2, wind speed was the most important variable (-0.76) (Digital appendix 3). In the dry season, wind speed had the highest correlation value in root 1 (0.89), and relative humidity was the most important variable in root 2 (0.95), followed by temperature (-0.79) and heat index (-0.66) (Digital appendix 3). The magnitude of the correlations for the most important variables in the dry season was higher with respect to the main variables of the wet season (Digital appendix 3).

Microclimate characterization in the wet season: Distribution of canonical scores (microclimate observations in each plot) was plotted only for the most important variable on each root. The microclimatic conditions during the wet season were more homogeneous, since the separation of canonical scores as a function of the variables was not evident (Digital appendix 4). Overall, in the horizontal axis (root 1), both the conserved areas and the sites of four and 17 years of succession presented average heat index values, close to the origin of the graph. However, in conserved areas, there was a tendency towards lower values of the variable, while in categories of four and 17 years of succession plots were located in areas with higher heat index values. In areas with 31 years of succession there was also a trend towards higher values of heat index than those presented in the other categories (Digital appendix 4A). Regarding the vertical axis (root 2), the highest wind speed values were recorded in areas with four, 17 and 31 years of succession, whereas these values were significantly lower in conserved areas (Digital appendix 4B).

A more specific microclimate characterization of the succession categories was performed based on the position of the centroids; that is, according to the canonical average of the microclimatic values observed in the plots, and was also interpreted only with the two most important variables. In the wet season, differences between succession categories were more evident when analyzing the centroids (Digital appendix 5). Thus, sites with four and 17 years of succession had a heat index close to 40.7 °C, but differed in wind speed since in areas of lower succession average values ranged between 1.25 and 1.5 km/h, whereas in the category of 17 years the values were lower, between 1 and 1.25 km/h. The highest heat index value was obtained in areas with 31 years of succession and fluctuated from 42 °C to 43.29 °C, with a wind speed of 1.25 to 1.5 km/h. For both variables, the lowest values were obtained in conserved areas, with an average heat index of 38.09 °C and a wind speed between 0.25 and 0.5 km/h (Digital appendix 5).

Microclimate characterization in the dry season: Separation of canonical scores, and therefore differentiation between succession categories, was more evident in the dry season (Digital appendix 6). In the horizontal axis (root 1), early succession sites were concentrated mainly in the right end of the graph, so they had higher wind speed values; contrarily, areas of 17 and 31 years of succession were scattered towards the left quadrant, indicating lower wind speed. Most of the conserved areas were located in the center and left of the canonical graph, displaying a low or average wind speed (Digital appendix 6A). In the vertical axis (root 2) the separation of the canonical observations was higher, and exhibited a progressive decrease of relative humidity from the conserved areas where the highest values were presented, towards the plots with four years of succession where the lowest humidity was obtained (Digital appendix 6B).

During the dry season (Digital appendix 7), the areas with four years of succession were characterized by a wind speed of 1.75 to 2.25 km/h, and a relative humidity of 60.79 to 62.2 %. Categories of 17 and 31 years presented values of wind speed below 0.25 km/h; however, relative humidity in areas of 17 years ranged from 62.2 to 63.59 %, while in the category of 31 years values were close to 65 %. Finally, most conserved areas recorded a wind speed of 0.25 to 0.5 km/h, and the highest relative humidity value, which oscillated between 66.4 and 67.8 % (Digital appendix 7).

Discussion

According to the observed results, early succession areas were related to higher temperature and lower relative humidity, in contrast to areas that have acquired more vegetation cover, i.e. conserved areas. Changes in relative humidity agree with similar studies accomplished on different vegetation cover types (^{Hojdová, Hais, & Pokorný, 2005}). However, temperature was not a significant microclimatic variable for separation of succession categories in this study, despite that the importance of this variable has been demonstrated in other chronosequences (^{Norris et al., 2012}) and forest studies with different degree of perturbation (Hojdová et al., 2005; ^{Hardwick et al., 2015}). Instead, heat index was one of the variables that contributed most to discriminate between categories in both seasons. Aside to the results here obtained, there is no antecedent (to our knowledge) of inclusion of heat index in other microclimate analysis associated with succession. However, its importance has been demonstrated in ecological niche studies for species of leaf beetles (Chrysomelidae) (^{Sánchez-Reyes et al., 2015}) and disturbance gradients (^{Sandoval-Becerra et al., 2017}). This variable measures thermal stress in humans by combining temperature and humidity, and it is frequently used to calculate the “apparent temperature”, also known as thermal sensation (^{Anderson et al., 2013}; ^{Lee & Brenner, 2015}). Therefore, differences found in the succession categories attributable to this variable, suggest a synergistic effect of humidity and temperature on the environment and consequently, on the species present in the study area.

Variation in microclimatic succession patterns is attributed mostly to the incidence of light, modified in turn by vegetation structure. In recently disturbed or early succession areas, low or absent tree cover allows higher direct solar radiation, rising air temperature and reducing humidity; over time, increase in canopy complexity reduces the incidence of light (^{Lienard, Florescu, & Strigul, 2015}) and generates different temperature and humidity conditions (^{Hardwick et al., 2015}). Consequently, modifications on those variables are reflected in changes in heat index. Wind speed was also an important variable in separating succession categories in this work. This variable is considered of relevance for transpiration of plants under the canopy (^{Renaud et al., 2010}). Its variation stems from the presence of trees and vegetation (Renaud et al., 2010; Hardwick et al., 2015). As a normal pattern, wind speed decreases with the increase in complexity of vegetation structure (Hardwick et al., 2015), being consistent with the results observed in this study since plots with earliest succession time showed highest values in wind speed. The significance of this variable in our study area may be related to its influence on dispersal of plant propagules (^{Damschen et al., 2014}; ^{Schurr, Bond, Midgley, & Higgins, 2005}), which is transcendental for the succession process. Accordingly, it is possible that higher wind values in recently disturbed areas may influence successional routes by introducing plant species characteristic of other vegetation types. Such influence needs to be analyzed in future studies.

Seasonality is a very important factor in succession analysis, since variation in abiotic climate factors in a given area contributes to post-disturbance resilience processes (^{Bernhardt-Römermann et al., 2011}). This work provides evidence on the influence of seasonality in microclimate along a successional gradient. During the dry season, relative humidity (RH) was related to successional differences indicating that different categories of succession may retain different RH conditions. Contrarily, as humidity is not a limiting factor in the wet season, it would be expected that other variables gain importance, as was observed with heat index. Such microclimate variations during succession are of major importance, since these abiotic constraints at a local scale are related to species distribution (^{Checa, Rodriguez, Willmott, & Liger, 2014}). Moreover, microclimatic seasonality is determinant for the development and establishment of vegetation (^{Lebrija-Trejos et al., 2011}).

However, although the effect of seasonality on microclimate during succession was evident in this study, results rejected our hypothesis. Significant differences amongst all succession categories were observed in the dry season, whereas in the wet season 17 years succession areas showed similar microclimate values as those of four and 31 years succession. In addition, variables correlation values were lower during the wet season than in the dry season. Other studies carried out in dry forests in Mexico have shown contrary patterns to our findings, since microclimatic heterogeneity is higher during the wet season (^{Lebrija-Trejos et al., 2011}). These discrepancies may be attributed to minor seasonality conditions in Northeastern Mexico compared to the dry forests of the Mexican Pacific Coast (^{Challenger & Soberón, 2008}; ^{Ceballos et al., 2010}; Trejo, 2010). This suggests a higher heterogeneity in characteristics and vegetation structure between categories, bring about by greater variation in the number of plant species and individuals that resist leaf fall during the dry season along the succession gradient. It has been observed that leaf area has a strong effect on microclimate, and that a high leaf area index gives rise to cooler environments (^{Baker et al., 2014}; ^{Hardwick et al., 2015}). Therefore, the increase in plant productivity and development of new leaves during the wet season may reduce microclimate variation at intermediate areas of succession, because of a more homogenous vegetation structure. ^{Bernhardt-Römermann et al. (2011}) report similar results, they found that resilience after disturbance is lower in dryer or arid areas. Accordingly, in the dry season differences between succession categories must be higher, as shown in this research.

Although vegetation type herein studied is categorized as low thorn forest (INEGI, 2013), because of its characteristics it could also be considered within the dry forest category (^{Challenger & Soberón, 2008}; ^{Ceballos et al., 2010}; ^{Trejo, 2010}). Among the most important characteristics of these ecosystems, in contrast to similar Neotropical forests, is its strong seasonality that causes low or null precipitation during the dry season months (Ceballos et al., 2010; Trejo, 2010). Nonetheless, particular conditions in the study area, such as minor seasonality, presence of certain plant species [i.e. Ebenopsis ebano (Berland.) Barneby & J. W. Grimes or Havardia pallens (Benth.) Britton & Rose] and geographic location (^{Miranda & Hernández, 1963}; ^{Treviño-Carreón & Valiente-Banuet, 2005}) seem to provide for the observed microclimatic variations. This may be because intensity of the gradient depends on similarity of conditions between succession stages as a function of the habitat structure (^{Letcher et al., 2015}).

In both seasons, it is evident that microclimate variables measured at sites with four, 17 and 31 years of succession are different as those measured in conserved areas. It is possible that these fragments may require longer periods to acquire abiotic conditions that resemble those found prior to disturbance. However, the intermediate succession categories (17 and 31 years succession), like categories four and 17 years of succession, exhibited similar microclimatic conditions during the wet season; this decrease in climate heterogeneity suggests that vegetation at the study area may experience a convergent successional trajectory (^{Walker, Wardle, Bardgett, & Clarkson, 2010}), and microclimatic factors could be restored over time (^{Baker et al., 2014}). In this way, it is possible that vegetation also changes with succession time till reaching the most conserved areas condition. This is of relevance for conservation, since distribution of native low thorn forest in the region has been reduced dramatically in the last years, and the existent fragments are surrounded by disturbed areas or sites with different succession time (^{Sánchez-Reyes et al., 2017}). Therefore, it is necessary to assess whether this trend is maintained in plants as well as in other biological groups within the study area, or on the contrary, if current disturbed areas with early succession time may not recover over time to its original condition.

In summary, we found microclimate differences between successional categories, these results agree with similar findings in other works (^{Hojdová et al., 2005}; ^{Hardwick et al., 2015}). However, microclimatic influence was higher in the dry than in the wet season, in contrast to other vegetal communities (^{Lebrija-Trejos et al., 2011}). Although we presented evidence for such findings, it is possible that the sampling design may have influenced the results obtained, and therefore a higher number of replicates may be required for a better microclimatic characterization. Thus, further studies would be necessary in other areas of this endangered ecosystem to corroborate our conclusions. Regardless of the observed patterns, we document the importance of assessing microclimate changes during succession, since this is a key factor on species distribution (^{Swanson et al., 2011}; ^{Arroyo-Rodríguez et al., 2017}). On the other hand, previous evidence suggests that successional trajectories of communities have very different ending points according to several factors (^{Bhaskar, Dawson, & Balvanera, 2014}; ^{Longworth et al., 2014}; ^{Letcher et al., 2015}); contrarily, successional microclimatic changes can be relatively well predicted for various ecosystems. For example, early successional areas are characterized by higher temperature and lower humidity (Hojdová et al., 2005; Hardwick et al., 2015; Letcher et al., 2015), while conserved stages are related to opposing microclimates (^{Chu et al., 2007}; ^{Lebrija-Trejos, Pérez-García, Meave, Bongers, & Poorter, 2010}; ^{Campetella et al., 2011}). In the study area, the low thorn forest patterns found are similar to those reported in other works (Hardwick et al., 2015; ^{Lienard et al., 2015}), but arose because of other variables of higher importance. Indeed, such patterns may originate due to vegetation characteristics, which we point out as a main driver of the observed differences. Besides, the fact that microclimate significant differences were greater during the dry season is of major relevance, since other dry forests shown contrary outcomes (Lebrija-Trejos et al., 2011). It demonstrates that influence of seasonality is distinct according to vegetal communities evaluated. As a consequence, our results suggest that assessment and restoration strategies of low thorn forest in Northern Mexico must be designed particularly for this ecosystem. In this regard, other studies have noted that conservation strategies cannot be the same for all areas, and therefore, each particular ecosystem and its full spectrum of characteristics associated with succession must be taken into account (^{Christensen, 2014}). This work may represent one of the first assessments and characterization of microclimate in low thorn forest in Mexico, demonstrating that observed successional patterns are particular for this ecosystem, and as such provide new information for conservation strategies of this endangered vegetal community.

Ethical statement: authors declare that they all agree with this publication and made significant contributions; that there is no conflict of interest of any kind; and that we followed all pertinent ethical and legal procedures and requirements. A signed document has been filed in the journal archives