Chronic Disturbance in a Tropical Dry Forest: Disentangling Direct and Indirect Pathways Behind the Loss of Plant Richness

Chronic disturbance is widely recognized as one of main triggers of diversity loss in seasonally dry tropical forests (SDTFs). However, the pathways through which diffuse disturbance is acting on the forest are little understood. This information is especially demanded in the case of vanishing Neotropical seasonally dry forests such as the Tumbesian ones. We proposed a conceptual model to analyze the factors behind the loss of woody species richness along a forest disturbance gradient, explicitly considering the existence of direct and indirect causal pathways of biodiversity loss. We hypothesized that the chronic disturbance can act on the woody species richness directly, either by selective extraction of resources or by browsing of palatable species for livestock, or indirectly, by modifying characteristics of the forest structure and productivity. To test our model, we sampled forest remnants in a very extensive area submitted to long standing chronic pressure. Our forests cells (200 × 200 m) were characterized both in terms of woody species composition, structure, and human pressure. Our structural equation models (SEMs) showed that chronic disturbance is driving a loss of species richness. This was done mainly by indirect effects through the reduction of large trees density. We assume that changes in tree density modify the environmental conditions, thus increasing the stress and finally filtering some specific species. The analysis of both, direct and indirect, allows us to gain a better understanding of the processes behind plant species loss in this SDTF.


INTRODUCTION
Chronic disturbance generated by anthropogenic practices of low intensity, such as selective logging (Singh, 1998) and free-range herding of livestock (Méndez-Toribio et al., 2016), is extended widely throughout most forests, but it is especially pressing in the case of seasonally dry tropical forests (SDTFs) (Álvarez-Yépiz et al., 2008;Ribeiro et al., 2015;Jara-Guerrero et al., 2019). However, this persistent and low-intensity degradation is receiving far less attention than deforestation (Álvarez-Yépiz et al., 2008), even though it has been reported a precursor for diversity loss and deforestation (Martorell and Peters, 2005;Ribeiro et al., 2015). This lack of attention is probably a consequence, unlike what happens with acute degradation, e.g., rapid conversion of land use, of the fact that the diffuse effect of chronic disturbance is very difficult to estimate via traditional approaches, often based on satellite imagery (Ribeiro et al., 2015).
In the tropics, the negative effects of this low-intensity disturbance on species richness of woody plants are widely recognized (Gillespie et al., 2000;Anitha et al., 2009;Williams-Linera and Lorea, 2009;Espinosa et al., 2011). Evidence suggests the existence of direct effects on the diversity as a function of the prevalent perturbation (Anitha et al., 2009;Espinosa et al., 2011;Clark and Covey, 2012;Ribeiro et al., 2015;Ribeiro-Neto et al., 2016). For instance, selective logging affects directly species with high wood density (Ribeiro et al., 2019), whereas domestic livestock usually eliminates palatable species (Anitha et al., 2009) and/or species with low wood density, and with no resprouting capacity (Sagar et al., 2003;Clark and Covey, 2012). This selective use of species leads to the elimination of species with particular traits, thus reducing species richness and functional diversity.
There is also a general consensus that diffuse forest disturbance generates changes in most structural features of the forests (Clark and Covey, 2012;Caviedes and Ibarra, 2017) which could exacerbate such diversity loss via indirect pathways. In general, two indirect pathways can be addressed as drivers of the species richness loss, one through changes in density of large trees, and other through the modification and simplification of physical structure of the forest. The first indirect pathway implies changes in stem density of large trees due to selective logging (Fredericksen and Mostacedo, 2000;Chazdon, 2003), together with decays in the smaller stem classes due to bark stripping from livestock during the dry season. These primary effects trigger changes at ground level (Fredericksen and Mostacedo, 2000;Espinosa et al., 2016), i.e., light availability Chazdon, 2001, 2002), soil water retention capacity (Ludwig et al., 2005), or litter decomposition rates (Fujii et al., 2017). Consequently, these changes lead to a loss of diversity due to modification on the dynamics of plant recruitment and mortality of established plants. The amelioration of the microclimatic conditions by tree cover is particularly important in stressful ecosystems both to improve the success of recruitment (Espinosa et al., 2016;Gusmán-M et al., 2018), and to reduce the mortality of saplings due to water limitation, which frequently intensifies in cleared disturbed forests, producing dieback or mortality of poor plastic or adapted tree species (Olson et al., 2018;Hofhansl et al., 2020). All these processes affect the demography of the coexisting tree species, indirectly exacerbating the degradation of the forest, and being able to induce a kind of diversity loss vortex.
The second indirect pathway for the loss of diversity is through the variation of the three-dimensional spatial arrangement of plants in the forests or structural complexity (Ishii et al., 2004;McElhinny et al., 2005). The forest complexity is affected by the chronic disturbance through the elimination of some canopy strata, such as large trees removed by selective logging (Clark and Covey, 2012), and/or modifications of the saplings bank because they differentially consumed by livestock (Méndez-Toribio et al., 2016). The changes in forest complexity can affect the forest productivity because the maintenance of a multi-layered canopy structure promotes complementary resource utilization among species (Ishii et al., 2004). Thus, a simplification of the structure is coupled with a decrease in the primary productivity leading to a reduction of species with different photosynthetic characteristics (Ishii et al., 2004;Pau et al., 2012). Although the causal relationships between species richness and productivity remain under discussion, for the SDTFs, it has been proposed that productivity also determines the plant species richness in these forests (Madrigal-González et al., 2020). Additionally, changes in the structural complexity can also induce a loss of diversity without affecting the primary productivity, due to an increase of generalist species that maintain the productivity but exclude specialist species of some canopy strata (Ishii et al., 2004;Clark and Covey, 2012).
Although there is evidence showing the individual negative effects of the above-mentioned factors on tropical forest diversity, joint analysis of all of them is necessary because the simultaneous action could exacerbate, attenuate or reverse their isolated effects on diversity loss. Here, we evaluated the effect of chronic disturbance of a seasonally dry neotropical forest, probably one of the most endangered forest ecosystems at the global scale (Miles et al., 2006), on the loss of woody species richness. With this in mind, we have built a complete causal scheme in which both direct and indirect relationships on species diversity are considered (Figure 1). We propose that the effect of chronic disturbance on woody species richness is exerted through different pathways acting together. Firstly, direct effects are expected as a result of selective elimination of species through logging and or livestock grazing. Secondly, indirect effects are expected as a consequence of changes both in the stem density but also in the loss of complexity which in advance affects richness through variations in the primary productivity (Pau et al., 2012). We assume that, under chronic pressures and changes in tree density, primary productivity can be maintained as a result of an increase in the density of other life forms, such as shrubs, even in those areas with low density of large tree species, but the productivity can shift in response to changes in structural complexity (Pau et al., 2012). Therefore, we did not consider pathways from tree density to productivity and structural complexity, nor a pathway from chronic disturbance to productivity. Based on this model we hypothesized that indirect pathways could explain the loss of species richness better than direct ones because (i) the loss of species through selective elimination can be compensated with the establishment of generalist or invasive species, such as those dispersed by livestock (Espinosa et al., 2021). (ii) There is evidence that in this stressful ecosystem the environmental filtering and positive plant-plant interactions determine the community structure and plant diversity (Espinosa et al., 2016;Gusmán-M et al., 2018), thus the exacerbation of habitat filtering exerted by changes in the micro-environmental conditions could have strong consequences on the plant species richness. FIGURE 1 | General conceptual causal model for direct and indirect pathways of loss of species richness. Arrows represent causal pathways. Disturbance represents the chronic disturbance measured as the first component of the PCA of distance to populated centers, roads and farms. The model considers that chronic disturbances can affect the species richness of woody plants (Richness) directly, or through three indirect pathways: reduction of large tree density (Density), reduction of the forest structural complexity, which is based on the tree height and basal area (Complexity), or reduction of primary productivity, measured as normalized differentiation of vegetation Index (NDVI) (Productivity), as a consequence of changes in the structural complexity.

Study Area
The study was carried out in the SDTF of southwestern Ecuador, specifically in Zapotillo County, Loja Province (Figure 2). This area is part of the Tumbesian biogeographic region, recognized worldwide for its high level of endemicity, as well as for the extreme loss and forest degradation for decades (Best and Kessler, 1995;Mittermeier et al., 2005;Tapia-Armijos et al., 2015). The mean annual temperature is 22.5 • C and the average annual rainfall is 642 mm/year, which is mainly restricted to the 4 months of rainy season (January-April) . The other 8 months of the year correspond to a dry season, were monthly precipitations are below 10 mm (Maass and Burgos, 2011). As in most SDTFs of the Neotropic, this forest has been used by local people for long for fuelwood and timber extraction and livestock rising, mainly goats (Jara-Guerrero et al., 2019).

Study Design and Forest Sampling
Within an area of 12,400 ha, we haphazardly located three grids of 4.2 × 2.2 km to have a reasonable variation of the forest diversity and heterogeneity in the region. Grids were separated at least 2.5 km from each other and with a distance of 5-8 km from Zapotillo town, with a population of approximately 5,000 inhabitants (GAD cantonal Zapotillo, 2019). Each grid covered a matrix of forest remnants with different histories of land use (Supplementary Figure 1), as well as agricultural areas and populated centers and thoroughly divided in cells of 200 × 200 m. Within each of these grids, we selected randomly 20 cells (Figure 2A). This study design allowed us to avoid the effects of confounding factors such as elevation, soil type, or topography, because this well-replicated and extensive approach explicitly includes areas with high heterogeneity in disturbance levels and current ecological factors. Within each selected cell we located three transects of 50 × 2 m (Figure 2B), separated by 20 m from each other and parallel to the main slope to maintain the highest homogeneity along the transect.
In each transect, we registered all intercepted trees and shrubs with diameter at breast height (DBH) ≥1 cm and height >1.30 m, and sub-shrubs with height ≥30 cm. Shrubs and sub-shrubs were only recorded in the central band of each transect, along 50 × 1 m ( Figure 2B). For each individual we recorded coordinates (x, y), total height, DBH, crown length in two directions (x, y) and crown height. Total height and crown height were measured using a handheld laser hypsometer/clinometer (TruPulse R , Laser Technology Inc.). Crown height was calculated as the difference between the total height and the height of the lowest branch of the individual. All these measurements were made during the rainy season between November 2018 and February 2019.

Forest Characterization at the Cell Scale
For each selected cell we defined an estimate of anthropogenic disturbance based on the distance to populated centers, roads and farms. We adjusted a principal component analysis (PCA) with standardized values with all these variables and used the first component as a chronic disturbance intensity variable. In SDTFs, the ease of access to humans into the forest is a good surrogate of chronic disturbance (Cueva Ortiz et al., 2019). In this region, goats are released every day to forage freely into the forest, and with the sunset, they return to the stables located near the owner's houses (Espinosa et al., 2021). Therefore, the distance of the populated centers to the forest is directly related to the logging and livestock loading history that the forest has supported (Cueva Ortiz et al., 2019). Considering that chronic disturbance implies a change in the intensity of pressures on the forest over time, we considered that the distance to populated centers is a better proxy of chronic disturbance than direct measurements of current pressure, such as livestock loading and/or recent logging. To characterize the forest structural complexity we used the square root index (SQRI) proposed by Barbeito et al. (2009). The SQRI is based on the sum of the square roots of the tree height differences and down weighting by basal area. In this index, the amount of values that are different acquires more importance in the value of the index than the magnitude of the differences and there is no need to define classes when it is used.
To calculate the plant productivity of each cell, we used the Normalized Differentiation of Vegetation Index (NDVI), which allows a better discrimination between plant cover with greater or lesser biomass production. This index yields estimated values of the amount of green foliar biomass (Ke et al., 2015) that are directly related to the photosynthetic capacity of the forest, and therefore, to the absorption of energy by plant cover and the forest structure (Pau et al., 2012). NDVI ranges from −1.0, where photosynthetic activity is absent (areas of bare soil) to 1.0 in areas with maximum green vegetation. We derived the NDVI data from a Landsat image with a resolution of 30 × 30 m. For each cell of 200 × 200 m we calculated the NDVI as the mean of all 30 × 30 m pixels. The last variable was the density of large trees, considering all stems >25 cm DBH in the three transects of each cell.

Structural Equation Modeling
We used structural equation models (SEMs) to test if our conceptual causal model for species richness is supported by our field data. We adjusted the model by using the variables at the cell level, pooling the data of the three transects. This conceptual model represents hypothesized causal relationships between chronic disturbance on woody species richness (Figure 1). SEM allows determining the relative importance of the direct effect of chronic disturbance and its indirect effects through each of the three variables, density of large trees, forest structure complexity and productivity on the species richness of woody plants (see Iriondo et al., 2003). We adjusted a SEM for the species richness of all woody plants and other one for trees using the function "sem" of the package lavaan (Rosseel, 2012) from the R environment. All variables were previously scaled using the "scale" function from the package stats (R Core Team, 2021).
Since the standard SEM assumes normality in the data (Lefcheck, 2016;Gana and Broc, 2019), we used the "mardiaKurtosis" function of seemTool package (Jorgensen et al., 2021) to test the normality of the data. The data of our two models violated the normal distribution assumption, thus we used a robust Maximum Likelihood (MLM). The MLM incorporates a scaling correction of Satorra-Bentler χ 2 (SBχ 2 ) for adjustment of the models (Gana and Broc, 2019). As SEM adjust a linear relationship between variables (Lefcheck, 2016;Gana and Broc, 2019), we evaluated the linearity of relationships between variables of the model through graphical representations. We evaluated the discrepancies between observed and expected covariance matrices using a chi-square goodness-of-fit statistic (χ 2 ). Adequate model fit was presumed with P-values above 0.05, which indicate that the observed and the modeled covariance under the causal relationships imposed by our model were not significantly different. Since this test is prone to several statistical errors, especially when the sample size is large like in our case, we also used two additional measures of the goodness of fit (Lefcheck, 2016;Gana and Broc, 2019), the comparative fit index (CFI) (Iriondo et al., 2003) and Tucker-Lewis Index (TLI) (Gana and Broc, 2019) with values above 0.9 suggesting adequate fit.

RESULTS
Fifty-one of our sixty sampling cells presented natural vegetation. The other nine cells corresponded to croplands and were excluded from the analyses. Overall, we recorded 2,873 woody individuals, belonging to 37 species, from which 19 species were trees and 18 shrubs and subshrubs. The most abundant families were Fabaceae and Malvaceae with 12 and 4 species, respectively (Supplementary Table 1).
Structural equation models did not show discrepancies with the observed patterns neither for total woody species richness (χ 2 = 1.36, df = 3, P = 0.72), nor for tree species (χ 2 = 1.35, df = 3, P = 0.72). The goodness of fit of both models (GOF = 0.99) suggested an adequate fit between our hypotheses and observed data. Our SEMs showed that chronic disturbance has a direct but weak negative significant effect on the total woody species richness ( Figure 3A and Table 1A). We did not detect significant effects on the tree species richness when analyzed separately ( Figure 3B and Table 1B). On the other hand, we found negative direct effects of chronic disturbance on the density of large trees and on the structural complexity (Figure 3). We also found positive direct effects, both from the density of large trees on the total richness, and from density of large trees and productivity on tree richness. We found an indirect effect of chronic disturbance through density of large trees on the total and tree richness ( Table 1). Overall, the indirect pathways better explained the loss  of total species and tree richness than the direct pathway from chronic disturbance ( Table 1). The joint effects explained the 0.30 of the variation for the total species richness, and the 0.53 of the variation for tree richness (Table 1).

DISCUSSION
Although current knowledge highlights the fact that anthropic pressures drive a generalized loss of species richness and forest degradation in STDFs (Gillespie et al., 2000;Álvarez-Yépiz et al., 2008;Anitha et al., 2009;Williams-Linera and Lorea, 2009;Espinosa et al., 2011) the pathways, direct or not, connecting chronic disturbance and diversity loss remained as a challenge. This is especially demanded in SDTFs, and in other forested regions in the Tropics, which are densely populated and dependence on forest ecosystem services is high. This demographic pressure requires a realistic biodiversity management which allows to reconcile conservation and human well-being (Martorell and Peters, 2005). Our work clarifies those direct and indirect pathways and the importance relating chronic disturbance to the loss of plant species richness. Our SEMs evidenced the complex nature of the relationships between the different drivers involved in the chronic disturbance and the richness of total woody species, and particularly of tree richness. The chronic and low-intensity disturbance showed a direct effect on the total woody species richness. This result agree with that proposed by Sagar et al. (2003), who suggested that selective logging generates a specific selective pressure in a SDTF of India, determining losses in diversity. Similar results were reported for a SDTF in Mexico (Williams-Linera and Lorea, 2009), where low richness was associated in a context of free range grazing to areas with high density of cattle trails. Likewise, in other SDTF of India, Anitha et al. (2009) found that high livestock disturbance reduced species richness through the removal of palatable species. Although our results support this issue, dissecting direct and indirect effect gives a better picture of the pathways involved in this loss. We found a significant direct effect of disturbance in species richness at the whole community scale, but not in the case of trees. These results suggest that direct effects of chronic disturbance on species richness is weak, and that its indirect effects as consequence of changes in microclimate conditions, resulting of large tree density reduction, are most important.
A common practice in managed forests is the selection of the largest trees, which can imply the elimination of particular species, as well as changes in the competitive dynamics of the remaining tree species (Clark and Covey, 2012). Williams-Linera and Lorea (2009) showed that canopy tree richness increases in STDFs with the basal area and density of canopy trees and, also, that seedlings and juveniles are highly dependent on water availability, which is directly related to density and complexity. Their results concur with ours because the density of large trees modifies the environmental conditions of the understory, thus generating a very selective environmental filter likely related to the reduction in the soil water availability in open areas (Balvanera et al., 2002). Previous studies in SDTF also demonstrated the close relationship between richness and stem density (Sagar and Singh, 2006;Williams-Linera and Lorea, 2009). The novelty of our result lies in the role of the density of the large trees as a link between chronic disturbance and the loss of species richness. Even more, this indirect pathway was by far a more critical causal pathway for species richness than the direct one (Figure 3 and Table 1).
Contrary to what we expected, our SEM did not show an indirect pathway on species richness through the structural complexity of vegetation. This result can be partially explained by an increase in the density of shrubs in more open canopy areas (for instance Ipomoea carnea). Since we used plant height as the variable of structural complexity index (SQRI), an increase in shrub density can compensate for tree losses. In addition, some disturbed areas are involved in recovery process through secondary succession once the fuelwood and more valuable wood resources have been harvested (Caviedes and Ibarra, 2017). Since many tree species of these forests have resprouting capacity, the structural complexity recovery can occur from a small subset of species very rapidly, i.e., without a necessary recover of species richness (Sagar and Singh, 2005;Álvarez-Yépiz et al., 2008;Balvanera et al., 2011;Clark and Covey, 2012).
Our results also concur with previous studies in SDTFs in which a relationship between primary productivity of vegetation and species richness was found (Sagar and Singh, 2006;Pau et al., 2012). However, the indirect pathway considering both the structural complexity and productivity of vegetation did not explain species richness (Table 1). This result contrasts with that reported by Pau et al. (2012), who found a correlation between NDVI and species richness with a strong effect of structural complexity. However, that study did not evaluate disturbance, and the reported changes were in response to water availability. Thus, harsher sites would be represented by a fewer subset of those species capable of dealing with hydric stress (Balvanera and Aguirre, 2006;Hofhansl et al., 2020), which can explain that significant effect of primary productivity on species richness of trees.
In summary, the chronic disturbances produced by anthropogenic activities of low intensity, such as selective logging and free-range grazing drive a loss of woody species richness mainly through indirect pathways. This means that, chronic disturbances are generating environmental modifications through the lessening of large trees, thus favoring the establishment of a small group of species that better tolerate drought and that direct pathway, which implies exclusion of particular species, showed a weak effect on species richness. These findings contribute to our understanding of how woody plant species are responding to chronic anthropogenic pressures, as well as for accurate predictions of how woody species will respond to future changes and exacerbation of these low intensity perturbation practices. However, these results should be taken with caution since we worked only in a region, being necessary to replicate the study in other areas and develop experimental evidence to test the importance of the different pathways behind the loss of species in the SDTF. This topic is particularly urgent in the Equatorial Pacific region, where the isolated condition of forest remnants, and the differences in the evolutive history of the local flora regarding other areas of Neotropical dry forest, can imply different responses to disturbance and, therefore, to require a different approach to its management.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the article/Supplementary Material.

AUTHOR CONTRIBUTIONS
AJ-G and CE conceived and designed the research and analyzed the data. DG-S collected the data. AJ-G, DG-S, and CE wrote the manuscript. All authors discussed the results and commented on the manuscript.

FUNDING
This work was supported by projects PIC-13-ETAPA-004 financed by the SENESCYT and managed by the Universidad de Cuenca, and PROY_CCNN_1054 financed by the Universidad Técnica Particular de Loja. This research was also partially funded by REMEDINAL-TE CM (P2018/EMT-4338) and PHENOTYPES (PGC2018-099115-B-I00) for AE.

ACKNOWLEDGMENTS
DG-S acknowledges Darwin Valle, María Fernanda Burneo, and Israel Medina for their field assistance. Thanks to anonymous referees whose criticism and suggestions improved a previous version of the manuscript.