From Stand to Forest: Woody Plant Recruitment in an Andean Restoration Project

1. Introduction

Rates of tropical forests degradation and deforestation continue growing worldwide [1]. In spite of the well-known effects of tropical forests shrinking in accelerating climate change effects and causing loss of global biodiversity and ecosystem services [2,3,4,5], more than half of the global deforestation rates occurs in the tropics [6,7]. This trend results in decreasing remnants of primary forests, surrounded by a complex landscape matrix of multiple land uses and secondary forests that impacts in primary forests connectivity and resilience [8,9,10]. In this scenario, complementing forest conservation efforts with ecological restoration actions can help to reverse and buffer tropical forest loss [11].

Ecological restoration aims to assist the recovery of an ecosystem that has been degraded, damaged, or destroyed [12]. This discipline relies on ecological succession theories to design restoration strategies that help accelerating (or enabling, in extreme cases) the ecosystem recovery. World tropical regions, harboring the highest levels of biodiversity and ecosystem complexity, and with more than their 50% of forested area belonging to secondary forests, are key to understand the successional pathways leading to a primary forest [11]. However, the multiple successional theories developed during the last century still need further research in order to unravel the mechanisms that drive ecological succession [13].

Forest active restoration consists in implementing specially designed actions for recovering forest ecosystems, and is especially suitable to test successional theories, since it allows to experiment with multiple factors that can affect ecological succession [13,14,15,16]. However, this potential of active restoration projects remains under-utilized due to the lack of rigorous, long-term monitoring [17,18,19]. Linking ecological processes with active restoration techniques through effective monitoring can broaden our knowledge of ecological succession in tropical forests and can contribute to understand the potential of restoration techniques for recovering its function and structure.

Colombia is one of the Earth’s most diverse tropical countries. Unfortunately, it is also an accurate example of the deforestation problems that we described above [20]. The country’s most populated area, the Andean region, entered the XXI century with the 69% of its original forest cover lost [21]. To address the deforestation problem, an increasing number of active restoration projects have been implemented in the region. However, most of them lack of a good monitoring system to prove if restoration objectives are met, and which ecological mechanisms are related to the successional pathways followed by the restored forests [22,23].

In this research, we conducted a woody plant recruitment study in an active restoration project of native Andean forest in northwestern South America (Colombia). Plant recruitment is one of the most important processes for plant succession, since it describes the spontaneous arrival and establishment of new plant individuals in a forest ecosystem [11]. This research aims to characterize woody recruitment density and diversity and to understand the mechanisms driving these attributes by relating them to multiple environmental factors: soil, planted trees structure and diversity, landscape metrics and previous land use. We hypothesized that woody plant recruitment density and diversity will differ across the study area in response to different environmental factors.

In particular, we expect that sites with fertile soils, higher planted trees diversity and development, immersed in a more forested landscape matrix and with less intensive previous land use will exhibit greater density and diversity of recruits. We expect recruitment features will respond to edaphic factors, since they affect plant germination, establishment and survival [11,24]. In that sense, better edaphic conditions will result in greater recruitment density and diversity. Since the recruitment process takes place in a forest plantation, we also expect recruitment features to interact with planted trees. Sites with a better development of planted trees will exhibit a greater recruitment density and diversity, since planted trees will offer more facilitation mechanisms to recruits, like an enhanced microclimate or reducing herb competition [16,25,26]. Also, higher planted trees diversity is expected to harbor a greater recruitment diversity and density, since they offer more diverse environmental conditions and can act as seed sources [27,28]. We expect that those planted areas surrounded by more forested landscapes will harbor greater recruitment density and diversity, since landscape matrix is a proxy of key ecological mechanisms like seed dispersal or availability [8,9]. Finally, we expect that recruitment features will differ between previous land uses, since land use have long-lasting legacies in areas undergoing ecological succession, we expect that those areas with more intensive land uses will have less recruitment density and diversity [24].

We conclude that all of these factors affect woody plant recruitment attributes to some extent: soil properties and planted trees diversity and structure influenced recruitment diversity, while landscape matrix and slope influence recruitment density. Planted trees structure influence both recruitment attributes, suggesting facilitation mechanisms between the planted trees and recruits.

2. Results

2.1. Recruitment Characteristics

We found 1903 recruits in the 22 sampled plots, averaging 80 (±59) individuals per plot and 0.4 individuals·m-² (Figure 1). Recruitment density varied greatly among plots, with a maximum of 201 individuals in plot 25 (0.10 individuals·m-²) and a minimum of 11 individuals in plots 3 and 5 (0.05 individuals·m-²). Density also varied among prior land uses, although these differences were not significant (Kruskal-Wallis, χ² = 5.60, G.L. = 2, p-value = 0.06), there is a marginal effect that indicated less recruitment density in plots with prior presence of exotic pastures (Pennisetum clandestinum) compared with native pastures (Dunn test, Z = -2.73, p-value = 0.02).

For recruitment diversity, we found 131 operational taxonomic units: 99 species and 32 morphospecies (6% of total recruitment), belonging to 56 genera and 34 families (Table A1). The most abundant genus were Miconia, Palicourea, Myrsine and Verbesina. The most abundant families were Melastomataceae, Asteraceae, Rubiaceae and Primulaceae.

Most of recruit diversity and abundance corresponded to native species that were not related to planted species. Only 24 recruited species (18,34% of recruit diversity) were part of planted trees species, while 67.77% of recruited individuals did not belong to planted tree species. Ten most abundant species accounted for 65.45% of recruited individuals, with endozoochory (60%) and anemochory (40%) dominating the dispersal modes (Table 1).

For Hill numbers, minimum sample coverage were 79.75%, which means that 79.75% of recruit community belongs to the 131 species and morphospecies found in this study (Table A2). Mean values of effective species were 7.42 (±3) for species richness (D₀), 4.9 (±1.83) for species diversity (D₁) and 3.63 (±1.42) for species dominance (D₂). None of the three Hill numbers showed significant differences between plots or prior land use, indicating no differences in species richness, diversity and dominance in the 22 sampled plots (Figure A1).

2.2. Recruiment Density and Diversity and Environmental Factors

Denser canopies, steeper slopes and a more forested landscape matrix favored recruits density. Recruitment density was best explained by canopy cover (β = 6.17 ± 2.60), plot slope (β = 0.02 ± 0.01), forest cover proportion in 1 km radius (β = 0.01 ± 0.01) and sand content (β = -1.54 ± 0.90) (Figure 2: Density). This model explains 65% of observed variation in data. Canopy cover, plot slope and forest cover in 1 km present a positive and significant effect in recruitment density (p value < 0.05), while sand content has a negative, non-significant effect (p value: 0.11).

Recruits species richness was also higher in plots with denser canopies. Additionally, physical soil properties interact with chemical soil properties and with planted species richness, establishing a spectrum of habitat conditions that influences recruits species richness. Linear mixed model for species richness (D₀) shows and effect of canopy cover (β = 40.98 ± 9.78) and the interaction of bulk density with planted species richness (β = -2.94 ± 1.01) and CEC (β = 12.24 ± 3.72) (Figure 2: D₀). Those plots with a denser canopy contain a larger number of recruited species. The plots with high bulk density and high CEC and, conversely, low bulk density and low CEC show greater recruitment species richness. Finally, those plots with high bulk density and less planted species, or plot with low bulk density and more planted species also show greater recruitment species richness. The model explains 80.40% ($R_c^2$) of observed data variation. Random factor (5 km buffer) has a variance (${\sigma }_a^2$) of 0.692.

Recruits species diversity (D₁) and dominance (D₂) are affected by similar environmental factors, being favored by sites with acid pH, low aluminum concentration and lower mortality of planted trees (Figure 2: D₁ and D₂). Linear mixed models show that edaphic conditions of pH (β = -37.11 ± 12.05 for D₁, β = -23.97 ± 10.58 for D₂) and aluminum (β = -6.91 ± 2.10 for D₁, β = -5.09 ± 1.89 for D₂) had a negative, significant effect in recruitment diversity and dominance, as well as the maximum mortality rate of planted trees (β = -38.08 ± 9.85 for D₁, β = -24.07 ± 8.81 for D₂).

Canopy cover and planted trees composition also influenced species diversity (D₁). Canopy cover has a significant, positive effect in recruitment diversity (β = 15.30 ± 6.47 for D₁). The second principal coordinate (Figure A2) shows a positive, significant effect in recruitment diversity (β = 3.24 ± 1.27 for D₁). The positive values of this axis represent plots with more abundance of planted species Montanoa quadrangularis, Retrophyllum rospigliosii, Psidium sp. and Quercus humboldtii, while negative values are associated with plots where the species Croton magdalenensis is more abundant. Plots with denser canopy cover and high abundance of M. quadrangularis, R. rospigliosii, Psidium sp. and Q. humboldtii harbor a greater recruitment diversity. Planted species richness shows a positive, non-significant effect on recruitment dominance (p value: 0.07) (β = 0.19 ± 0.09 for D₂).

The models for D₁ and D₂ explain 67.70% and 61.50% ($R_c^2$) of observed data variation respectively. Random factor (5 km buffer) has a variance (${\sigma }_a^2$) of 0.662 and 0.612, respectively.

3. Discussion

3.1. Recruitment Density and Environmental Factors

Our study shows that a ten year-old restored Andean forest can achieve successional trajectories that are very similar to those found in secondary Andean forests. Recruitment density in this study is greater than density found in secondary forests of Medellin municipality (for DBH > 10 cm, 0.2 individuals·m^-2 versus 0.1 individuals·m^-2 found by [35]) and in passive regeneration lowlands of Costa Rica (for DBH > 5 cm, 0.27 individuals·m^-2 versus 0.1 individuals·m^-2 found by [36]). However, recruitment in primary Andean forests is much denser, achieving values of 4 individuals·m^-2 [37]. These values can be interpreted in light of the stand forest development hypothesis, in which vegetation density declines when the canopy closes, light becomes a limiting resource and competition causes the death of many individuals, a phase known as stand exclusion stage [13]. The less dense recruitment in secondary forests could indicate that they have already entered the stand exclusion stage, whereas our restored forest, where canopy closure has not occurred yet, could be in the previous phase, where vegetation is denser. However, understory vegetation density in primary Andean forest seems to be much denser and both secondary and restored forests are far from those values.

Sites with denser canopies and steeper slopes surrounded by a more forested landscape substantially increased recruitment density. Since canopy closure has not occurred yet in the studied plots, it can offer appropriate microclimatic conditions while allowing direct light to reach the understory. Facilitation effects of planted trees canopy over recruited community have been found in other studies and include enhanced microclimate or reducing competition with grass species [16,24,25]. At landscape level, more forested landscape matrices provide a greater seed pool and facilitates faunal movements, which increases seed dispersal, favoring recruitment processes [24,38,39]. Since animals disperse the 60% of the most abundant species in this study, faunal movements play an important role for recruited community in the restored forest. Positive effect of forest matrix in species richness recovery has been demonstrated, but not in vegetation density [9]. Our study shows that forest matrix can also affect recruitment density and suggests that the role that landscape matrix plays in tropical forest secondary succession is broader.

Land use did not explain the differences observed in recruitment density among the plots. However, the positive effect of slope in recruitment density could be a proxy for land use intensity, since areas with steeper slopes tend to be deforested later in Colombia [40,41]. The role that land use plays in forest succession remains one of the biggest gaps in successional studies, mainly because the information about its type, intensity and duration is often scarce [24]. We could not find differences among the three previous land uses present in our study, but the marginal negative effect of foreign pastures (Pennisetum clandestinum) in recruitment density suggests that land use legacies deserve a closer focus in Andean forest restoration studies.

3.2. Recruitment Diversity and Environmental Factors

Recruitment diversity values were also comparable to those found in other secondary forests in Colombian Andes. We found 131 species in total, a value of species richness that is equivalent to secondary forests nearby our study area: 109 species in secondary oak forests (Quercus humboldtii) located in Medellin [42], and 190 woody plant species in a secondary forest of Central Andean Cordillera [43]. The most abundant species, genus and families present in this study are reported as common pioneer species in Andean forest succession [42,44,45], and even in Neotropical succession [11]. These results show that Andean forest restoration efforts in Medellin rural areas can be an effective strategy for recovering the diverse characteristics of secondary Andean forests.

Edaphic factors influenced the three diversity measures considered in this study. Cation Exchange Capacity (CEC) and bulk density had a positive influence on species richness through their interaction. Bulk density was directly associated with clay content, which plays a major role in soil CEC. Recruitment species richness peaked at extremes, while recruitment in plots with average conditions of CEC and bulk density had less species richness. Diversity and dominance were negatively influenced by aluminum concentration and pH, which means that these diversity measures declined with high aluminum content and basic pH. Aluminum has toxic effects for plant root development at high concentration and low pH conditions and increased mortality rates of planted trees in our study area [46,47]. However, in our study low pH increased recruitment diversity and dominance, which could indicate an adaptation of recruitment species to acid soils, since they are common in the Andes [48]. Further research is needed for better understanding the mechanisms that relate these edaphic factors with recruitment diversity. Exploring recruitment species composition and demographic rates could be useful for unveiling these mechanisms.

Our study also demonstrated the presence of facilitation effects between planted trees and recruited individuals. Canopy cover had a positive effect in species richness and diversity, probably through the mechanisms that also affect recruitment density. Planted trees maximum mortality rate had a negative effect on recruitment diversity and dominance, indicating that a loss of tree cover in plantations limited the species that can survive in the understory. Planted species composition was also important for recruitment diversity, which was greater in plots where Montanoa quadrangularis, Quercus humboldtii, Retrophyllum rospigliosii and Psidium sp., while the presence of Croton magdalenensis seemed to reduce recruitment diversity. Another study in this restoration project found different functional traits associated with ecohydrological processes that separated C. magdalenensis from Q. humboldtii and R. rospigliosii, which suggest that planted species composition could affect recruitment community through its functional traits [49].

Planted species richness had an effect in recruited species richness through its interaction with bulk density. These results seemed contrary to the Intermediate Disturbance Hypothesis, which states that diversity peaks at intermediate environmental conditions [50]. Research on recruitment demography rates and species composition may shed light on the mechanisms behind this pattern, for example, intermediate conditions could favor the growth of highly competitive species that displace species richness peaks to the more extreme conditions where it cannot grow.

3.3. Implications for Andean Forest Restoration

This study shows that, in a decade of restoration efforts in tropical Andes, the resultant forest can harbor a recruitment community equivalent to those found in Andean secondary forests, and that this community is mostly composed by native species that are not related with planted trees. Since our study area is located in the periphery of Medellin, the second largest city of Colombia, this research suggest that forest ecological succession can progress in densely populated areas, where landscape matrix is characterized by a great fragmentation and the coexistence of many types of land uses. In this project, forest restoration took place in former croplands and cattle ranches with native or foreign grasses. Previous land uses showed no effect for recruited community attributes. Since recruitment is a strong proxy for forest regeneration potential, it could indicate that restoration activities are suitable for these three types of land use in tropical Andes.

In the study area, seed dispersal, and specially fauna mediated dispersal, played a major role compared to seed bank or planted trees reproduction. Landscape forest cover facilitated faunal movements and we found that it also increased recruitment density, so regarding forest cover nearby areas undergoing restoration efforts could help to choose and adapt restoration techniques.

The poor contribution of planted trees to recruited community could be explained by their young age. Long-term research will help to understand their role in forest succession and if their contribution to plant recruitment increases when they reach their reproductive stage. However, long-term studies of recruitment demography rates are recommended to understand interactions between planted trees and recruited individuals, since the frontier between facilitation and competition mechanisms is complex and determinant for forest community assembly [16].

This study detected a planted trees facilitation effect over recruitment density and diversity, mainly through canopy cover, so choosing species with dense and broad canopies could help the recruitment process in restoration projects, especially in areas where cattle ranching took place and grasses can compete with young trees. Choosing species with these kind of functional traits can allow the implementation of cost effective restoration techniques like applied nucleation [51]. Since this study has also shown that the identity of planted species matter for recruitment processes, further research in functional traits of Andean tree species could help identify which traits enhance recruitment processes, thus accelerating restoration success.

Monitoring restoration activities in tropical Andes is a great opportunity to learn about forest succession in a biodiversity hotspot characterized by highly degraded and populated areas. We encourage long-term monitoring and experimental design of restoration techniques in order to create a synergy between ecological restoration and ecology that broadens our understanding of tropical montane forests while improves their extension and resilience.

4. Materials and Methods

4.1. Study Area

This study took place in rural areas of Medellin, Colombia (6°15′00.0”N 75°34′05.0”W), where an active restoration project of native Andean forest (Mas Bosques para Medellin) is being implemented from 2010, with more than 600 ha restored to date. This project consists in the establishment of mixed plantations of 115 species (91% native from Colombia) in areas with previous land agronomic uses (crops and cattle farming). Planted trees are monitored biannually since 2011 in a net of 70 permanent plots of 707 m² each, in which every planted tree is marked and coded.

We conducted the woody plant recruitment study in 2021 in the 22 oldest permanent plots, which were 10 years old when the recruitment measurement was done (Figure 3). The altitudinal range of selected plots is 2010 to 2690 AMSL, with mean annual values of 14.5ºC for temperature and 2200 mm for precipitation, which places them in the wet montane forest life zone [52]. The dominant soil types in the study area are Andisols. We measured as recruitment all woody plant individuals > 80 cm height present in an 8 m radius from the plot center (201 m²). For each recruit, we measured total height in cm, DBH for individuals taller than 130 cm, and we identified it to the highest taxonomic level possible. Each recruit was marked with an unique ID for future measurements.

4.2. Response Variables

We selected four response variables for describing recruitment characteristics. For each plot, we calculated recruitment density (No. individuals·m^-2), and for diversity we estimated the first three Hill numbers (D₀, D₁, and D₂) [53]. We chose Hill numbers because is a parametric family of diversity indexes that obey a replication principle, they have units (effective number of species) that enable easy comparisons between assemblages, and they have equivalences with classic diversity measures like species richness (D₀), Shannon-Wiener entropy (D₁) and Gini-Simpson dominance (D₂) [53]. We standardized Hill numbers with a measure of sample completeness, the minimum sample coverage found in the plots (79.75% of sample coverage), to reduce abundance bias in diversity measurements [54]. We performed all diversity measurements with iNEXT package in R software [55,56].

4.3. Predictor Variables

4.3.1. Edaphic Conditions

We collected five soil samples per plot and homogenized them for a fertility analysis, in which we measured pH, CEC, organic matter, phosphorus, aluminum, copper, iron, manganese, zinc and soil texture (sand, silt and clay). We also collected one sample per plot with a metallic cylinder for measuring bulk density.

4.3.2. Plantation Structure and Diversity

We used the database with the results of planted trees monitoring in years 2015, 2017 and 2019. We calculated the following measurements for each plot: basal area per plot in 2019, average mortality rate and maximum mortality rate (year 2017) [57], species richness (D₀) in 2019 and the first two axes of a Principal Coordinate Analysis that we performed for planted species composition in 2019 following [15] (p. 4) (Figure A2). Additionally, we included from this database plot elevation, slope and previous land use, which had three categories: crops, native pastures and exotic pastures (Pennisetum clandestinum Hochst. ex Chiov.). We also measured canopy cover on each plot with hemispherical photography with a Nikon Coolpix 950 and an 8 mm Sigma lens. We made five photos per plot and analyzed them with the software HemiView 2.1 to obtain mean canopy cover per plot [58].

4.3.3. Landscape Metrics

We carried out a visual classification of current land cover using a Medellin orthophoto from 2021 (ESRI World Imagery, 15 m of resolution). We distinguished three types of land covers: 1) Forested areas, 2) Crops and pastures and 3) Urban areas. Following the analyses proposed by [9] (p. 4), we classified the land cover in a radius of 1, 2, 3 and 4 km around the plots. We then calculated landscape metrics for each radius with the software V-LATE 2.0 in ArcMap [59]. For each land cover class, we calculated total area, proportion, mean patch size, total edge and mean patch edge. We also calculated total number of patches, dominance and edge density for the three land cover classes combined (landscape level analysis).

4.4. Statistical Analyses

To relate the recruitment density and diversity with soil properties, plantation structure and landscape metrics, we adjusted four lineal models, one for each response variable. For density, we adjusted a generalized linear model (GLM) with Poisson-distributed errors. We detected subdispersion in our standard errors, so we corrected them using a quasi-GLM where the variance is given by φ × μ, where μ is the mean and φ the dispersion parameter [60]. For the three Hill numbers (D₀, D₁ and D₂), we adjusted three linear mixed models. We included a random factor because we detected spatial autocorrelation in our response variables within a 5 km radius, so we classified the plots in groups of 5 km radius and included the resultant five groups as the random factor in our mixed models.

Since we had a large dataset comprising numerous environmental variables, we performed a preselection of the three type of explanatory variables in order to avoid multicollinearity problems. For edaphic variables, we first made a Principal Component Analysis (PCA) and included the first and second axes in our linear models, but these axes were systematically excluded from them. We excluded those edaphic variables that were strongly correlated (clay and silt with sand and bulk density), and then adjust a saturated linear model with all edaphic variables (Equation 1) and perform a backward selection using the stepAIC procedure from MASS package in R [61]. The stepAIC method selects the best model based on AIC values, but we also check the AIC for small samples (AICc) since our number of explanatory variables is higher than the number of observations.

y= β₀ + β₁ logpH + β₂ MO + β₃ logP + β₄ logAl + β₅ logCEC + β₆ arcsenA + β₇ bulk_dens + β₈ Fe + β₉ Cu + β₁₀ logMn + β₁₁ logZn + ε

(1)

Where logpH: log-transformed pH, MO: organic matter (%), logP: log-transformed phosphorus, logAl: log-transformed aluminum, logCEC: log-transformed CEC, arcsenA: arcsin-transformed sand (%); bulk_dens: bulk density, Fe: iron, Cu: copper, logMn: log-transformed manganese, logZn: log-transformed zinc.

Landscape metrics were strongly correlated between them, so we selected the landscape metric that showed the highest correlation with response variables. For this, we constructed a set of linear mixed models including one of the landscape metrics as a covariate, following [15] (p. 3). We selected the model with the lowest AICc, which was the one that included forest proportion within a radius of 1 km. We included this landscape metric along with the plantation variables and performed another linear model backward selection with stepAIC and AICc (Equation 2).

y= β₀ + β₁ PCoA1 + β₂ PCoA2 + β₃ canopy + β₄ BA + β₅ mean_mort + β₆ max_mort + β₇ D₀_19 + β₈ elevation + β₉ slope + β₁₀ land_use + β₁₁ forest_1km + ε

(2)

Where PCoA1: first principal coordinate, PCoA 2: second principal coordinate, canopy: canopy cover, BA: basal area in 2019, mean_mort: mean mortality rate (2015-2019), max_mort: maximum mortality rate (2017), D₀_19: planted species richness 2019, elevation: plot elevation in m, slope: plot slope in %, land_use: prior land use per plot, forest_1km: forest proportion in 1 km radius.

Finally, we adjusted a saturated model with the preselected edaphic variables and the preselected plantation variables and built a saturated model with them and all their possible interactions. We perform the last backward selection with stepAIC and AICc to obtain the final model for each response variable. We also checked that all final models had a Variance Inflation Factor (VIF) lower than 10 [60]. Model assumptions were checked graphically following the method proposed by [60].

All the statistical analyses were performed in R software [56]. We constructed the mixed models with nlme package [62], and the GLM with stats package. We standardized all variables in order to compare the magnitude of its partial regression coefficients with base package. We calculated the AICc with qpcR package [63].

5. Conclusions

This study demonstrates that one decade of active restoration techniques of Andean forests can recover woody plant recruitment attributes that are equivalent to secondary Andean forests. Successional trajectories did not differ among plots or previous land use, indicating that restoration techniques can be effective in multiple deforestation contexts. Recruitment was mostly composed by native species that were not related to planted trees, but planted trees affected recruitment attributes through facilitation processes, which indicates that restoration can enhance natural trajectories of secondary succession in the Andes.