Plant Biodiversity Homogenization across the Chronosequence in Highly Fragmented Landscapes in the Colombian Andean–Amazonian Transition

Abstract

Succession in abandoned pastures in the tropics can progress along different pathways, and the changes in biodiversity on local and landscape scales, as well as in species turnover (β diversity), are still unclear. This study assessed the role of secondary forests as a plant biodiversity reservoir or as evidence of a pattern of biodiversity homogenization in a chronosequence of abandoned pastures in two highly fragmented landscapes (hills and mountains) in the Colombian Amazon. In each landscape, the plant community composition, growth habits, species richness accumulation, indicator species identification, composition dissimilarity, and influence of above- and below-ground environmental filters were evaluated in four successional stages: (i) degraded pastures (DP) (<3 years old), (ii) early forests (EF) (10–20 years old), (iii) intermediate forests (IF) (25–40 years old), and (iv) old-growth forests (OF) (>90 years old). A total of 918 species, 78 families, and 11,284 individuals were recorded. The most representative growth habits were trees and shrubs, while palms and lianas were minimal. The species accumulated rapidly in DP and EF, contrary to that observed in IF and OF; thus, DP and EF achieve inventory completeness faster than IF and OF. It was necessary to include more plots to obtain inventory completeness in IF and OF. OF had a high α diversity with similar species composition throughout (low β diversity) and high homogeneity, suggesting biotic homogenization. On the other hand, IF presented the highest species differentiation (high β diversity) and a higher divergence in species composition than OF. The spatial distance and environmental dissimilarity had the most important role in determining species composition. Finally, our results show divergence in the plant composition across the chronosequence, where DP was opposite from OF in hills. However, in mountains, DP followed the EF and IF categories. The deviation in the species composition in EF and IF suggests an exchange of species in intermediate forest ages.

1. Introduction

Secondary forests cover ~77% of the forested global area, are more common than primary forests in tropical regions and will be the main tropical forest cover in the future [1,2,3]. In addition, because of the global distribution and ecosystem services, secondary forests are becoming more important for biodiversity conservation and sustainable management than primary forests [3]. However, the increase in secondary forests could be facilitating the decline of local biodiversity through a local individual abundance decline [4]. Empirical evidence has shown that species composition may take centuries to recover [5], and the potential of recovery varies as the result of landscape and regional influences [3] and species-specific attributes [6].

Recent efforts to synthesize the available evidence suggest fast tropical forest recovery during secondary succession with variation across regions where the composition of species in slow plant communities is recovered (>12 decades) [7]. Despite the rapid recovery of site-level diversity, plant community assemblages may take much longer [8].

The average range size and diversity, evenness, abundance, and species identity between communities could be distinct, especially with more contrasting ages in secondary forest. Thus, an expected pattern was that wide-ranged plant species colonized secondary forests in contrast to narrow-ranged species (resident species), leading to system homogenization [9].

Biotic homogenization was defined as “the increase in the taxonomic similarity of two or more biotas over a specified time interval” [10]. Such definitions could be interpreted as the reduction of species turnover (β diversity) without the influence of differences in species richness [11,12] and are a major signal of the start of the Anthropocene, the current human-dominated geological epoch [13]. The decline of β diversity between sites with different abandoned use ages could express the degree of biotic homogenization [11]. In addition, the reduction in the spatial turnover in species composition (low β diversity, more than local diversity loss) was related to a low spatial turnover in functions that showed negative consequences for large-scale ecosystem functions [14]. Biotic homogenization has been reported in aquatic and terrestrial taxa and in most of the world’s ecosystems [11,12,15], but it has been poorly studied in highly species-rich plant communities such as the Amazonian Forest.

Secondary forests in the Amazon represent ~0.5% (234.795 km²) of global secondary forests and are barriers to forest recovery in highly deforested areas [16]. In addition, ~78.2% of secondary forests are 20 years old (i.e., secondary intermediate forest), and the median age is 8 years (i.e., secondary early forest) [16]. Despite evidence of local diversity across the chronosequence in tropical forests [5,17,18], knowledge on replacement in plant diversity is scarce; it is not possible to trace signals of biota homogenization.

In the Andean–Amazonian region of Colombia, forests are changing continuously, and processes such as deforestation have caused changes in the climate region and land use [19,20]. Thus, extensive livestock farming increased by 31.22%, accompanied by an expanding agricultural frontier, the establishment of illicit crops, mining, and illegal logging [21]. Another determining factor in the regeneration pattern of vegetation cover is soil. The evidence suggests broad variation in the physical, chemical and biological conditions of the soil during secondary succession [22,23]. The wide heterogeneity of soils and their response to changes in cover generates a multifactorial scenario that may govern recovery trends for vegetation cover during natural regeneration. Deforestation acts synergically with climate change in the Amazon, pushing the highest biodiversity reservoir to potentially irreversible biodiversity losses.

Overall, in tropical regions, there is a trend to abandon the land, which gives way to natural regeneration. Although secondary forests have the potential to buffer ecosystem decline [24], the turnover of and the effect on local species diversity remains ambiguous, and the full stage of the deforestation effect is unclear. In addition, the heterogeneous degree of soil degradation generates a continuum of secondary forests and areas with evidence of arrested succession. Thus, the effect of secondary forests as evidence of biotic homogenization is still unclear. Understanding succession processes and the pattern of plant species recovery will provide tools for conservation efforts in core remnant forest areas that maximize beta diversity (i.e., as a proxy for species turnover) through natural regeneration, along with incorporating enrichment efforts with species that show limited ability to reach some remnant forest areas.

We are not focusing on secondary forest trajectories, but on the pattern of biodiversity α and β to trace the evidence of plant biodiversity homogenization. We considered the limitation in natural systems by following a single recovery trajectory, even within the same forest or region [25,26]. Thus, we observed the multiple sites with different ages of abandonment and multiple successional stages [27,28].

This study evaluated the partition between local and regional diversity to determine the role of secondary forests as evidence of plant biodiversity homogenization, showing a pattern of forest succession in the chronosequence framework. Two widely distributed landscapes (i.e., mountains and hills) were evaluated. The descriptions focused on diversity attributes, along with dissimilarity measurements of the plant community across the chronosequence of abandoned cattle pastures in two highly fragmented landscape units (hills and mountains) in the Colombian Amazon. Given the evidence of the secondary forest biodiversity and recovery capacity, a high local and regional plant biodiversity was expected in the old forest, independent of the forest landscape, maintaining regional biological plant biodiversity.

2. Materials and Methods

2.1. Study Area and Sampling Design

The study area was located in the northwestern part of the Colombian Amazon, in the municipalities of San José del Fragua (1°′19‴52″ N; 75°′58‴28″ W, at an elevation of 540 m above sea level), Belén de los Andaquíes (1°′24′ 5″.1″ N; 75′52′ 2″.2″ W, at an elevation of 320 m above sea level), Morelia (1°′29‴09″ N; 75°′43‴28″ W, at an elevation of 258 m above sea level), and Florencia (1° 36′50″ N; 75° 36′ 46″, at an elevation of 242 m above sea level), in Caquetá, with two highly fragmented landscapes (hills and mountains) in the Andean–Amazonian transition (Figure 1). Hills corresponded to an undulated topography in areas with mainly livestock use and some forest relicts, while mountains had higher slopes and a mosaic of natural forests, pastures and agricultural land, as described in previous studies [23]. The climate in the study area was humid and warm (Caldas–Lang classification), characterized by a monomodal regime where the rainy period is between March and June, with an average annual precipitation of 3376 mm and an average temperature of 25.04 °C. This area has oxisols and ultisols—soils with low fertility and aeration, acids (pH < 6), high aluminum saturation (>60%), and a mainly clayey texture [29].

In each landscape unit, a chronosequence with four different successional categories was established as described in previous studies [23]: (i) degraded pasture (DP) (<3 years old), composed of a mixed cover with degraded Brachiaria spp. grasses (e.g., koronivia grass and signal grass), weeds from the Poaceae and Cyperaceae families, and some shrub species, mainly dominated by Miconia spp.; (ii) early forest (EF), forests between 1 and 10 years old; (iii) intermediate forest (IF), forests between 25 and 40 years old; and (iv) old-growth forest (OF) or mature forests over 90 years old. In each age category, 50 × 50 m plots (0.25 ha) were established, for a total of 33 plots in the area of study (14 plots in the hills, and 19 plots in the mountains) (Table 1).

The chronosequence approach used in this study, grouping plots by age categories given the time of abandonment, has been widely used and analyzed by analysis of variance (ANOVA) in various studies [17,30,31]. For example, Rocha-Ortega and García-Martínez [32] defined a chronosequence with three successional categories (SC): SC l (6–12 years), SC 2 (13–27 years), and old-growth forest (OGF). Abbas et al. [30] established a chronosequence of five successional categories: 0–14, 14–26, 26–52, 52–70, and >70 years, respectively. Martínez-Ramos et al. [17] defined five successional categories: Pasture (<1.5), Early (3.5–5.5), Mid (6–8), Advanced (A, 13–15), and Old-Growth Forest (OGF, forest sites without signs of human disturbance). Hu et al. [31] defined a chronosequence with three categories: the early successional stage (present day), middle successional stage (hundreds of years), and late successional stage (ca. 1300 years). In addition, it is also important to highlight that the study of chronosequence can follow two approaches: (i) successional category (grouping plots with ages within the same category or age range of abandonment), from which the effect of chronosequence on the different attributes of biodiversity is analyzed by means of ANOVA [17,30,31,32]; or (ii) plots with different ages of abandonment along the successional trajectory, which are analyzed using regression models [33,34]. Therefore, both approaches are valid for analyzing the effect of chronosequence on changes in floristic attributes through secondary succession.

2.2. Floristic Data Collection

The measurement of species composition and growth habits was undertaken from February to December 2017. All plots and trees were marked in the field. In each plot, all trees with a diameter at breast height greater than 1 cm (DBH) were counted and identified at the species level, and the diameter at DBH was measured. All specimens were processed, identified, and deposited in the “Herbario Amazónico Colombiano—COAH” of the Sinchi Institute in Bogotá DC, Colombia.

2.3. Environmental Parameters

In the present study, 28 environmental parameters were analyzed: two topographic, three climatic, ten soil chemical, eight soil physical, and five soil macroaggregation. The slope was measured on each plot using a laser hypsometer (Forestry Pro II), and the plot elevation was registered using a digital altimeter with GPS (GPSMAP 64CSX). The climatic variables were temperature and precipitation extracted from a WorldClim dataset [35] to each plot. The third climatic variable was environmental stress [36] defined by $\mathrm{E}=1.\mathrm{e}-3\times \left(0.178\times \mathrm{TS}-0.938\times \mathrm{CWD}-6.61\times \mathrm{PS}\right)$, where TS was temperature seasonality as defined in the WorldClim dataset (bioclimatic variable 4). CWD was the climatic water deficit (in mm/year, measured as above). PS was the precipitation seasonality as defined in the WorldClim dataset (bioclimatic variable 15). For the soil variables, five sampling points were selected in each plot to collect different soil samples and to analyze the physicochemical and morphological properties according to previous studies [23]. Thus, for the physicochemical variables, a 30 × 30 × 30 cm trench was made in each sampling point to collect soil samples at 10 cm increments at a depth of 30 cm. Then, the soil data from 0–10, 10–20, and 20–30 cm layers were pooled and averaged. The following physical parameters: penetration resistance (MPa), texture (sand, clay, and silt) (%), bulk density (g cm⁻³), total porosity (%), soil moisture (%), and structural stability index (SI) (%) were evaluated. The chemical parameters were pH, cation exchange capacity (CEC) (meq 100 g⁻¹), electric conductivity (EC) (dS m⁻¹), available phosphorus (P) (mg kg⁻¹), exchangeable acidity (EA) (mg kg⁻¹), soil organic carbon content (SOC) (%), total nitrogen (N) (%), calcium (Ca) (mg kg⁻¹), magnesium (Mg) (mg kg⁻¹), and potassium (K) (mg kg⁻¹). The soil macroaggregates were evaluated in each sampling point with a 10 × 10 × 10 cm soil monolith, and the aggregate morphology was classified as follows: biogenic macroaggregates, root macroaggregates, physical macroaggregates, and non-macroaggregates. Other soil components such as leaves, roots, and vegetable fragments were classified as organic material [23,37]. The soil macroaggregation data for the 0–10 cm layer were based on an average value of the five sampling points [23].

2.4. Data Analysis

The sampling intensity efficacy was evaluated with rarefaction analysis using the function rarefaction from the R package mobr [38]. The species abundance distribution (SAD) was also analyzed. Two rarefaction methods were used: individual-based rarefaction [39] and spatial sample-based rarefaction with the k-nearest centroid neighbor (kNCN) method [40]. Then, the biodiversity indices: observed species richness (S), rarefied species richness (Sn), effective number of species (ENS) based on PIE (probability of interspecific encounter) (SPIE), and number of individuals or total abundance (N) were calculated with get_mob_stats from the R package mobr [38]. General linear models (GLM) using two-way ANOVA (α = 0.05) were implemented for testing the differences in the ecological indices, where the analyzed factors were successional category (DP, EF, IF, and OF), landscape (hill and mountain) and their interaction, taking into account the replications (plots) used for each factorial combination (Table 1). The response variables were fit for normal distribution with square root transformation. The assumptions of normality and homogeneity of variances were validated with the exploratory residual analysis of the model. The Fisher’s LSD post hoc test was used to create confidence intervals for all pairwise differences between the factor-level means, while controlling for the individual error rate at a level of 0.05. In addition, permutation tests to compare the indices at the alpha scale, beta scale, and gamma scale were also performed using the get_mob_stats function. Indicator species values of each successional category in both landscapes were calculated using the indval function from the R labdsv [41].

Nonmetric multidimensional scaling (NMDS) analysis based on Bray–Curtis distance and 50 iterations was used to compare the dissimilarity within and between each successional category and each landscape unit. NMDS was performed with the metaMDS function from the R vegan package [42]. The analysis of multivariate homogeneity of group (successional category or landscape unit) dispersions based on Bray–Curtis distance (with α = 0.05 and 9999 permutations) was evaluated with the betadisper function from the R vegan package. Bray–Curtis distance was used since it considers abundance and can indicate changes in species composition as well as abundance. Furthermore, it can be accommodated by reducing the distances to principal coordinates, which embeds them within a Euclidean space. Therefore, it is a better measure of distance than the Euclidean distance [42]. The betadisper test has been useful in assessing beta diversity [43]. A permutational multivariate analysis of variance using distance matrices (adonis) [44] with Bray–Curtis distance and p-value (α = 0.05) based on 9999 permutations was used to test whether there were significant differences between the successional categories, landscapes, and their interaction. The adonis function from the R vegan package was used for this purpose. Then, multilevel pairwise comparisons between groups were performed with the pairwise.adonis function from the R pairwiseAdonis package [45].

Mantel correlation pairwise tests [46] between the species dissimilarity matrix (with Bray–Curtis distance) and different environmental matrices (standardized data and Euclidean distance) using the Spearman method and 9999 permutations were performed with the mantel function from the R vegan package. Two sets of analyses were employed on the 28 environmental parameters to remove the effect of autocorrelation and multicollinearity to analyze the relationships of environmental variables with the studied biodiversity parameters. First, the effect of environmental variables on plant communities was evaluated with a redundancy analysis (RDA) [47] using the rda function from the R vegan package. The plant community data were Hellinger-transformed. The species vectors and environmental variables that explained more of the variation between the plots were selected with the envfit function from the R vegan package.

Second, the relative importance of the chronosequence, landscape, spatial distance, and above and below ground environment variables (four explanatory matrices) on species composition was estimated with variation partitioning of the community data with adjusted R-squared in distance-based redundancy analysis (db-RDA) (with Bray–Curtis distance) [48,49] using the varpart function from the R vegan package. Thus, before variation partitioning, the following analyses were performed: (i) principal component analysis (PCA) was performed on 28 environmental variables using the dudi.pca function from the R ade4 package [50]. The permutation test for the RDA was significant. Therefore, the forward and backward selection (999 permutations) based on RDA was used to identify significant components using the ordistep function from the R vegan package. Seven significant components were selected. (ii) The spatial data matrix was transformed to the principal coordinates of neighbourhood matrix (PCNM) suitable for constrained ordination [49]. A permutation test (α = 0.05 and 999 permutations) for the global model based on RDA was performed; since it was not significant, the PCNM vector selection was not done. Therefore, a total of 16 PCNM vectors were used.

All analyses in R language were performed in R. 4.0.3 [51] using the interface in RStudio v.1.3.1093 [52]. The ANOVA and Fisher’s LSD tests were carried out in InfoStat software v.2020 (InfoStat Transfer Center, CBA, Argentina) [53].

3. Results

3.1. Plant Community Composition and Growth Habits

In this study, a total of 918 species, 78 families, and 11,284 individuals (Supplementary Table S1) were identified and recorded in 33 plots belonging to four successional categories in two landscapes (Table 1). In hills, 493 species, 63 families, and 4620 individuals were recorded. The five species with a higher number of individuals were: Siparuna guianensis Aubl. (361), Adenocalymma cladotrichum (Sandwith) L.G.Lohmann (188), Miconia elata (Sw.) DC. (167), Pseudosenefeldera inclinata (Müll.Arg.) Esser (136), and Henriettea fascicularis (Sw.) M.Gómez (128). The five families with more individual abundance were: Melastomateaceae (711), Monimiaceae (375), Rubiaceae (238), Euphorbiaceae (230), and Arecaceae (223). In mountains, 727 species, 75 families, and 6664 individuals were reported. The five species with more individual abundance were: H. fascicularis (461), Graffenrieda conostegioides Triana (196), P. inclinata (171), Wettinia praemorsa (Willd.) Wess.Boer (131), and Virola elongata (Benth.) Warb. (130). In addition, the five families with a higher number of individuals were: Melastomataceae (1500), Rubiaceae (543), Euphorbiaceae (461), Myristicaceae (348), and Flacourtiaceae (284).

The most representative growth habits in both landscape types were trees and shrubs, while palms and lianas were minimally represented (Supplementary Figure S1). In hills, the proportion of tree species increased from 25.45% in the degraded pasture (DP) to 65.84% in old-growth forest (OF), while shrub species decreased from 63.64% to 30.96% from DP to OF. Palms and lianas together accounted for less than 5% of the total species recorded. No lianas were recorded in DP, and no palms were recorded in early forest (EF). A very similar pattern was observed for the proportion of individuals reported. Thus, OF had the highest proportion of arboreal individuals (72.47%), while palms and lianas together were the least abundant (<10% of the total individuals).

In mountains, a similar pattern was evident in the proportion of species and individuals according to growth habit. The proportion of tree species varied between 46.05% and 57.40% for DP and EF, respectively, while shrubs decreased from 52.63% in DP to 42.01% in EF. Palms and lianas had a proportion of less than 6% of the total species recorded in all successional categories. In relation to the proportion of individuals, arboreal individuals were the most dominant (58.47% in OF to 78.03% in EF), contrary to lianas and palms, where palms reached a higher proportion in OF (8.72%), and lianas did so in intermediate forest (IF) (2.88%).

3.2. Species Richness Accumulation

Species accumulation was analyzed with a successional category and landscape unit (Figure 2). Along the chronosequence, the rarefaction based on individual and spatial samples (Figure 2a,b) showed that the species accumulated rapidly in DP and EF, contrary to that observed in IF and OF. At the landscape level, species accumulation was faster in hills than in mountains (Figure 2c); however, considering the spatial proximity of the samples, it was similar in both landscape types (Figure 2d). Additionally, the species abundance distribution (SAD) curves showed that some degraded pasture plots located in the mountains were highly dominated by a single species (Supplementary Figure S2). However, the evenness between successional categories as landscape units was very similar.

A two-way ANOVA showed that the interaction between chronosequence and landscape was not significant for the four diversity indices derived from the rarefaction analysis (p > 0.05). No significant effect of the landscape was evident on four indices (p > 0.05), and chronosequence had no effect on abundance (N) (F = 1.11, p > 0.05). In contrast, the species richness differed significantly between successional categories for observed species richness (S) (F = 12.24, p < 0.001), rarefied species richness (S_n) (F = 13.32, p < 0.001), and effective number of species based on the probability of interspecific encounter (S_PIE) (F = 5.71, p = 0.004). According to Fisher’s LSD test, OF showed significantly higher mean values of species richness (S = 116.42, S_n = 67.14, and S_PIE = 37.15; p < 0.01) than DP (S = 28.25, S_n = 20.62, and S_PIE = 6.17) (all p < 0.05).

Permutation tests on the alpha scale also confirmed these findings (Figure 3a,d,g,j; Supplementary Figure S3). On the beta scale, the higher richness values were only evidenced for S in IF (p < 0.05) (Figure 3b,e,h), while on the gamma scale, higher richness values were recorded for S and S_n both in OF and IF (both p < 0.01) (Figure 3c,f,i). The N did not vary significantly among the successional categories (p > 0.05) (Figure 3k). Furthermore, none of the indices varied significantly between landscape units on the beta and gamma scales (p > 0.05) (Supplementary Figure S3).

3.3. Indicator Species Identification

The indicator value (IV) analysis showed variation in the indicator species according to the landscape (Figure 4). In hills (Figure 4a), two indicator species were characterized in the plots with more contrasting ages (i.e., DP and OF). The DP plots were represented by Aegiphila parviflora Moldenke and Vismia baccifera (L.) Triana & Planch. (both, IV > 0.74, p < 0.05). For OF, ten species were identified, where Compsoneura capitellata (A.DC.) Warb., Hevea pauciflora (Benth.) Müll.Arg., Oenocarpus bataua Mart., Micropholis guyanensis (A.DC.) Pierre, Iryanthera laevis Markgr., Virola pavonis (A.DC.) A.C.Sm., and Protium leptostachyum Cuatrec. were more significant (all IV > 0.80, p < 0.01).

In mountains (Figure 4b), there were no indicator species in DP, while in EF, three indicator species were identified (Inga thibaudiana DC., Miconia sp., and Senna bacillaris (L.f.) H.S.Irwin & Barneby; all IV ≥ 0.50, p < 0.05). In IF, there was one species (Chrysochlamys dependens Planch. & Triana; IV = 0.63, p = 0.034), and in OF, there were 18 species, where Neea ovalifolia J.A.Schmidt, P. inclinata, and Endlicheria sprucei (Meisn.) Mez were more significant (all IV > 0.77, p < 0.01).

3.4. Plant Community Composition Dissimilarity

The adonis test and NMDS ordination (Figure 5) showed significantly different compositions in the plant communities for different successional categories (F = 1.87, R² = 0.15, p < 0.001) (Figure 5a) and landscape types (F = 1.70, R² = 0.05, p = 0.009) (Figure 5b). The plant communities from the DP and OF plots were more dissimilar (F = 2.67, R² = 0.18, p < 0.001), while the EF and IF plots were very similar (F = 0.82, R² = 0.05, p > 0.05) (Figure 5a). The OF plots were dissimilar from EF and IF (both, p < 0.001), as well as DP from IF (F = 1.38, R² = 0.08, p < 0.019). There were no differences between DP and EF (F = 1.11, R² = 0.10, p > 0.05).

In addition, the interaction between chronosequence and landscape was not significant (F = 0.96, R² = 0.08, p > 0.05). However, the partitioned adonis and NMDS analyses by landscape revealed some specific patterns in the community composition across the chronosequence (Supplementary Figure S4). In hills, a very similar dissimilarity pattern was observed between successional categories (Figure 5a and Supplementary Figure S4), while in mountains, although OF plots were dissimilar to the remaining successional categories (all p < 0.05), the DP plots did not differ from EF and IF (both p > 0.05).

The homogeneity of multivariate dispersions analysis (beta diversity) for each group based on average distance to centroid (Figure 6) showed that the species composition within each successional category and landscape unit varied in the chronosequence and topography, respectively.

Consequently, significant differences in the mean dissimilarity variance between successional categories (F = 4.32, p = 0.012) (Figure 6a) and landscape units (F = 6.77, p = 0.014) (Figure 6b) were observed. Furthermore, pairwise comparisons revealed a higher dispersion in IF than in OF (p = 0.003). In addition, a beta diversity analysis partitioned by landscape showed that the variances between successional categories were heterogenous in hills (F = 7.24, p = 0.007), but not in mountains (F = 1.52, p > 0.05) (Supplementary Figure S5). Thus, in hills, the variances in the IF plots significantly differed from those in DP (p = 0.011) and OF (p = 0.010).

3.5. Influence of Environmental Filters

Based on Mantel’s tests (Figure 7), the species Bray–Curtis dissimilarity matrix had a significant relationship with the spatial (i.e., physical separation) (Mantel statistic R = 0.14, p = 0.022) (Figure 7a), above-ground environment (Mantel statistic R = 0.13, p = 0.036) (Figure 7b), soil (Mantel statistic R = 0.21, p = 0.003), and age (Mantel statistic R = 0.43, p < 0.001) distance matrices (Figure 7b–d). As the plots became physically more separated or dissimilar in terms of environment or age, they also became more dissimilar in terms of plant community composition.

The differences in community composition at the plots were also significantly correlated with the elevation (Mantel statistic R = 0.11, p = 0.044), precipitation (Mantel statistic R = 0.14, p = 0.039), temperature (Mantel statistic R = 0.10, p = 0.044), pH (Mantel statistic R = 0.21, p = 0.004), K content (Mantel statistic R = 0.16, p = 0.015), cation exchange capacity (CEC) (Mantel statistic R = 0.11, p = 0.036), bulk density (Mantel statistic R = 0.15, p = 0.017), penetration resistance (Mantel statistic R = 0.23, p = 0.001), total porosity (Mantel statistic R = 0.15, p = 0.016), and organic material (Mantel statistic R = 0.25, p = 0.001) (Supplementary Figure S6).

The redundancy analysis (RDA) of the plant communities associated with the chronosequence or landscape, constrained by above- and below-ground environmental factors (constrained inertia = 89.68%) (Figure 8), revealed an evident separation of the OF plots from the other successional categories along the first ordination axis (RDA1 = 11.9%) (Figure 8a), as well as mountains from hills along the second axis (RDA2 = 6.0%) (Figure 8b). Nine environment variables (three above and six below ground) significantly explained more of the variation between the plots (all p < 0.05).

The OF plots were associated with soils with higher values of non-macroaggregates (r = 0.23, p = 0.020) and EC (r = 0.24; p = 0.016), while the plots with lower ages (DP and EF) were more associated with high values of pH (r = 0.26; p = 0.011) and penetration resistance (r = 0.29; p = 0.008). In relation to the landscape, the mountain plots that were associated with higher elevations (r = 0.35; p = 0.001) and slopes (r = 0.29; p = 0.011) also correlated with high K (r = 0.32; p = 0.003) and organic material (r = 0.35; p = 0.001) contents. In contrast, the hill plots with lower values of these parameters were strongly correlated with a high temperature (r = 0.46; p = 0.001). On the other hand, six species had a strong influence on ordination (all r > 0.5, p < 0.05). Thus, Theobroma subincanum Mart., M. guyanensis, and P. inclinata correlated strongly with OF, while S. guianensis, V. elongata, and Eschweilera coriacea (DC.) S.A.Mori were mainly associated with EF and IF.

In addition, an RDA analysis partitioned by landscape confirmed the separation of the OF plots from the other successional categories in both landscape units, and, in the case of hills, the DP plots separated from the EF and IF plots (Supplementary Figure S7). Additionally, in hills (with α = 0.1), the DP plots were directly associated with root macroaggregates (r = 0.27, p = 0.098), IF plots with biogenic macroaggregates (r = 0.40, p = 0.069), and OF plots with non-macroaggregates (r = 0.41, p = 0.059). The EF plots had lower penetration resistance values (r = 0.34, p = 0.092). In addition, seven species had a strong effect on plot ordination (all r > 0.7, p < 0.05). S. guianensis was more correlated with DP and EF; Gustavia hexapetala (Aubl.) Sm. and V. elongata were more correlated with EF and IF; and Miconia chrysophylla (Rich.) Urb., M. guyanensis, O. bataua, and Tetragastris panamensis (Engl.) Kuntze were more correlated with OF.

In mountains (with α = 0.05), the OF plots were highly associated with soils containing higher concentrations of organic material (r = 0.39, p = 0.016) and non-macroaggregates (r = 0.32, p = 0.034), while high Ca contents were evident in the EF plots (r = 0.32, p = 0.033). In general, the DP and IF plots had lower values for these three soil parameters. For species, six were more important in the ordination (all r > 0.7, p < 0.05). T. subincanum, M. guyanensis, P. inclinata, and E. sprucei correlated strongly with the OF plots, while Alchornea triplinervia (Spreng.) Müll.Arg. and Piper arboreum Aubl. were more associated with the other successional categories.

Finally, variation partitioning of the community showed that spatial distance and environmental (above- and below-ground) variables each explained 12% of the total variation in the species composition. In contrast, the chronosequence and landscape only explained 7% and 6% of the variance, respectively. In total, these four explanatory matrices together explained 25% of the variations, with 75% unexplained (i.e., residuals) (Figure 9).

4. Discussion

The highly fragmented landscape of Caquetá is home to almost 1000 species of trees and shrubs, reaching 111.27 per hectare. This diversity is an order of magnitude lower than that reported by Duivenvoorden [54] (313 spp/ha, dbh ≥ 2.5 cm) and three times less than that reported by Duque et al. [55] (183 spp/ha, dbh ≥ 2.5 cm) for forests with good conservation status in Caquetá. The low diversity showed the effect of the transformation of forests on the capacity for biodiversity.

The variation in the tree recovery rate from degraded pasture (25% and 46% in hills and mountains, respectively) to OF (65.84% and 48.47%) shows a bigger recovery capacity in the hill landscape than the tree species community in the mountains. Despite the differences in the recovery rates, they were in the range of values reported by Poorter et al. [7] for tropical forests. Species removal by growth habit showed the expected pattern, with an increase in the number of tree species with increasing forest age, with a decrease in shrubs [25].

The rapid initial accumulation of species in DP showed several woody species were able to occupy pastures, although there were high values of pH and soil penetration resistance and lower values of organic matter and non-macroaggregates [30]. However, the richness in DP showed lower values throughout the chronosequence, as did the β and γ diversity. This local and spatially reduced species pool represents a challenge for community assembly in this highly transformed region, because only a subset of species is available to contribute to the community ensemble [30]. Old forest (OF) had a high α diversity with similar species compositions throughout (low β diversity) (Figure 3a,b) and high homogeneity (Figure 6), suggesting biotic homogenization [56]. Land cover change and tree species logging could drive such a trend. In Colombia, ~47% of the timber trade consists of illegally felled trees [57] of specific species, which could be structuring the community assembly. The plant community homogenization also arises from the removal of native species, the introduction and propagation of non-native species, the intensification of land use, and the increased homogeneity of environmental resources [58,59,60]; however, the variability in the environmental condition in OF was lower within (OF at a geographical scale) than among forest ages. The convergence in species composition between the OF plots, in addition to the low dispersion of above- and below-ground environmental conditions, showed that the plant community was a “climax” forest [61]. On the other hand, a higher dispersion and β diversity were observed in the intermediate forest (IF), indicating floristic differentiation and divergence in species composition, similar to that reported by Abbas et al. [30] in mid-successional stages of evergreen rainforest in south China. Events such as the isolation produced between vegetation patches from deforestation [62], changes in land use and intensity from cattle ranching activities, and the absence of tree felling [60,63], environmental filters, limitations in dispersion, and recruitment of species [30] could drive such a trend in IF.

When the indicative species were analyzed, most species showed strong preferences in the mountain landscape. The Mantel test showed a significant relationship with the physical separation, above-ground environment, and soil conditions, which could represent pressure from the ecological processes of dispersal and seedling establishment, leading to species preference in the mountains. In hills, the indicator value showed that fewer species are recruited in relatively high abundance in degraded pastures (DP) (Figure 4) compared to the greater species number in the old forest (OF). A. parviflora and V. baccifera are “pioneer” species that positively respond to extreme conditions. Both species are shrubs that produce fleshy fruits with endozoochory dispersal. However, V. baccifera persisted in EF and IF, suggesting that the pioneer title may not accurately define the species preferences [30]. In addition, there was a lack of strong environmental filtering in the secondary forests in the hill landscape. No indicatory species was detected for the plots in EF and IF, suggesting randomness in species occurrences or a high turnover in species composition, as indicated by Abbas et al. [30]. In mountains, EF recruited I. thibaudiana, Miconia sp., and S. bacillaris, and IF did so for C. dependens. These species indicated strong specialization in the colonization in the secondary forests in the mountain landscape. C. dependens persisted in the old forests, reducing its dominance and being replaced by late-successional species. The set of indicator species in OF in both hills and mountains showed an opportunity to select species for enriching treatments in EF and IF.

There was divergence in the plant composition among the different successional categories, following the expected pattern, where DP was the opposite of OF. EF and IF were convergent, falling into a similar category. This pattern persisted in hills; however, in mountains, DP followed the EF and IF categories. The deviation in the species composition in EF and IF suggested an exchange of species that is typical in earlier stages of succession [61]. The RDA analysis indicated that the community composition relationships with the environment, DP, EF, and IF were on the negative side of the RDA1, showing the influence of pH and soil penetration resistance on the community composition. The pH and soil penetrations were typically altered by the cattle pasture use history of those forests [64]. Rodríguez-León et al. [22] reported that soil pH decreases and organic carbon and nutrients such as Ca and Mg increase along the successional trajectory, favoring the establishment of arbuscular mycorrhizal fungal communities that are strongly associated with less-disturbed ecosystems. Likewise, the physical and biological quality of the soil is favored by natural regeneration, which will have a direct impact on the increase of species richness and composition of the plant community [23]. The plant community in OF responded to non-macroaggregates and soil EC, showing the persistent effect of soil-use history on the plant communities’ assembly. The homogeneity analysis showed the same pattern throughout the chronosequence; it is a mature forest (i.e., OF) that is mostly homogeneous when compared to forests in earlier stages of succession (i.e., DP, EF, and IF).

Despite the chronosequence framework applied in this research, the forest age only explained a small proportion of the variation (7%) in the plant composition. In comparison, the proportion of variance jointly explained by the spatial distance and environmental variables was greater than that explained by age (24%), indicating a strong role of dispersal and/or colonization and the environmental filtering of plant community assemblies.

5. Conclusions

This study confirmed that mature forests have a trend of biotic homogenization (i.e., lower β diversity), similar to that observed in degraded pastures with less than 3 years of abandonment, although the species richness was significantly lower than in mature forests. In contrast, the intermediate secondary forest presented the highest species differentiation (i.e., higher β diversity) and a divergence in species composition, mainly because of environmental filters above and below ground, spatial proximity, and species recruitment, as evidenced by the change in indicator species throughout the chronosequence. In particular, a lower similarity in the plant community composition of the intermediate forests with respect to the mature forests was more evident in the hill landscape. This suggests that intermediate secondary forests may follow pathways of succession and species composition that do not necessarily converge with mature forests. In addition, it was also evident that spatial distance and environmental dissimilarity have an important and similar role in determining species composition, although more than half of the variation in the species composition was unexplained, mainly because of stochastic processes. Finally, active restoration plans are needed to conserve diversity in intermediate forests and to promote legal mechanisms for the conservation of forests that favor the conservation of areas with greater biotic homogenization in highly fragmented landscapes in the Colombian Amazon.

Author Contributions

Conceptualization, C.H.R.-L. and L.L.R.-F.; methodology, C.H.R.-L. and L.L.R.-F.; software, A.S.; validation, A.S., C.H.R.-L. and L.L.R.-F.; formal analysis, A.S.; investigation, C.H.R.-L., L.L.R.-F. and A.S.; resources, C.H.R.-L.; data curation, A.S. and L.L.R.-F.; writing—original draft preparation, C.H.R.-L., A.S. and L.L.R.-F.; writing—review and editing, C.H.R.-L., L.L.R.-F., A.S. and J.C.S.; visualization, A.S.; supervision, C.H.R.-L. and L.L.R.-F.; project administration, C.H.R.-L.; funding acquisition, C.H.R.-L. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

Data are available from the authors upon request.

Acknowledgments

The authors thank all of the farmers of the study area for their help and support during the fieldwork; Herminton Muñoz-Ramirez for his support in graphical editing; and Christopher King for reviewing the English of this manuscript.

Conflicts of Interest

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of the data; in the writing of the manuscript, or in the decision to publish the results.

Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Figures and Table

Figure 1. Location of the study zone and plots (Caquetá state, Colombian Andean–Amazonian transition). Abbreviations: DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest.

Figure 2. Species rarefaction curves. (a,c) Individual-based rarefaction (the species richness for each plot was rarefied based on 135 individuals); (b,d) spatial sample-based rarefaction (the species were accumulated by including spatially proximate samples first) with the kNCN method; (a,b) successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest); (c,d) landscapes.

Figure 3. Biodiversity indices derived from the individual rarefaction curves, and permutation tests to compare the indices among the successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest). (a–c) Observed species richness (S); (d–f) rarefied species richness (Sn); (g–i) effective number species (ENS) of probability of interspecific encounter (PIE) (SPIE); (j,k), abundance or number of individuals (N); (a,d,g,j) biodiversity at alpha scale; (b,e,h) biodiversity at beta scale; (c,f,i,k) biodiversity at gamma scale. The upper and lower whiskers represent the maximum and minimum values established at the third (+1.5 IQR) and first quartile position (−1.5 IQR), respectively.

Figure 4. Indicator species values of each successional category (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) in two landscape units. (a) Hill; (b) mountain. Species with a significant indicator value (p < 0.05) were indicator species.

Figure 5. Nonmetric multidimensional scaling (NMDS) based on Bray–Curtis dissimilarity displaying the composition of plant communities in different successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) (a), and landscape units (b). Ellipses represent the standard deviation around the centroid of each category or landscape.

Figure 6. Analysis of multivariate homogeneity of group (i.e., successional category: DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) dispersions (variances) based on the group centroid. (a) Successional category; (b) landscape unit. Asterisk (*) represent mean values. The upper and lower whiskers represent the maximum and minimum values established at the third (+1.5 IQR) and first quartile position (−1.5 IQR), respectively.

Figure 7. Mantel correlation pairwise tests between species dissimilarity matrix and different environmental matrices. (a) Species vs. spatial distance; (b) species vs. above-ground environment (climatic and topographic factors); (c) species vs. soil (i.e., soil physicochemical and macroaggregation parameters); (d) species vs. age (i.e., successional category midpoint). Scatterplots shown in (b), (c), and (d) display the spatial separation gradient.

Figure 8. Redundancy analysis (RDA) of plant communities associated with different successional categories or landscape units, constrained by above- and below-ground environmental variables. (a) Successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest); (b) landscape units. Ellipses represent the standard deviation around the centroid of each category or landscape unit. Blue arrows represent the environmental variables that were significant to constrained ordination (p < 0.05). Red arrows indicate the more important species vectors in the ordination (r > 0.5 and p < 0.05).

Figure 9. Partition of the variation of plant community by four explanatory matrices: (1) chronosequence: representing four successional categories; (2) landscape: corresponding to two landscape units; (3) spatial: representing 16 PCNM (principal coordinates of neighbourhood matrix); (4) environment, including both above- and below-ground variables.

Table 1

Number of plots for each successional category distributed in two landscape units.

Chronosequence	Age Range (Years)	Age Midpoint (Years)	Landscape	Number of Plots
Degraded pasture (DP)	<3	1.5	Hill	2
			Mountain	3
Early forest (EF)	10–20	15	Hill	3
			Mountain	4
Intermediate forest (IF)	25–40	32.5	Hill	6
			Mountain	6
Old-growth forest or mature forest (OF)	>90	90	Hill	3
			Mountain	6

Supplementary Materials

The following supporting information can be downloaded at: https://www.mdpi.com/article/10.3390/f13091422/s1, Figure S1: Proportion (%) of different growth habits along the chronosequence (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) in two landscape units. (a,b) Richness-based calculation; (c,d) abundance-based calculation; (a,c), hill; (b,d), mountain; Figure S2: Species abundance distribution (SAD) curves. (a,c) Standard rank curves; (b,d) cumulative curve; (a,b) successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest); (c,d) landscapes; Figure S3: Biodiversity indices derived from the individual rarefaction curves and permutation tests to compare the indices among the landscapes. (a–c) Observed species richness (S); (d–f) rarefied species richness (Sn); (g–i) effective number species (ENS) of probability interspecific encounter (PIE) (SPIE); (j,k), abundance or number of individuals (N); (a,d,g,j) biodiversity at alpha scale; (b,e,h) biodiversity at beta scale; (c,f,i,k) biodiversity at gamma scale; Figure S4: Nonmetric multidimensional scaling (NMDS) based on Bray–Curtis dissimilarity displaying the composition of plant communities in different successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) in two landscape units. (a) Hill; (b) mountain. Ellipses represent the standard deviation around the centroid of each category for each landscape; Figure S5: Analysis of multivariate homogeneity of group (i.e., successional category: DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest) dispersions (variances) based on the group centroid in two landscape units. (a) Hill; (b) mountain; Figure S6: Mantel correlation pairwise tests between species dissimilarity matrix and different environmental vectors. (a) Species vs. elevation; (b) species vs. precipitation; (c) species vs. temperature; (d) species vs. pH; (e) species vs. K content; (f) species vs. CEC (cationic exchange capacity); (g) species vs. bulk density; (h) species vs. penetration resistance; (i) species vs. porosity; (j) species vs. organic material. All of the scatterplots display the spatial separation gradient; Figure S7: Redundancy analysis (RDA) of plant communities associated with different successional categories (DP, degraded pasture; EF, early forest; IF, intermediate forest; OF, old-growth forest or mature forest), constrained by above- and below-ground environmental variables in two landscape units. (a) Hill; (b) mountain. Ellipses represent the standard deviation around the centroid of each category. Blue arrows represent the environmental variables that were significant to constrained ordination (p < 0.05). Red arrows indicate the species vectors most important in the ordination (r > 0.5 and p < 0.05). Table S1: List of the plant species and their growth habits by plot for each successional category in two landscape units.

References

1. Food and Agriculture Organization of the United Nations. FAO Global Forest Resources Assessment (FRA) 2020; FAO: Rome, Italy, 2020.

2. Pain, A.; Marquardt, K.; Lindh, A.; Hasselquist, N.J. What Is Secondary about Secondary Tropical Forest? Rethinking Forest Landscapes. Hum. Ecol.; 2020; 49, pp. 239-247. [DOI: https://dx.doi.org/10.1007/s10745-020-00203-y] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/33343057]

3. Chazdon, R.L.; Peres, C.; Dent, D.; Sheil, D.; Lugo, A.E.; Lamb, D.; Stork, N.; Miller, S. The Potential for Species Conservation in Tropical Secondary Forests. Conserv. Biol.; 2009; 23, pp. 1406-1417. [DOI: https://dx.doi.org/10.1111/j.1523-1739.2009.01338.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20078641]

4. Newbold, T.; Hudson, L.N.; Hill, S.L.L.; Contu, S.; Lysenko, I.; Senior, R.A.; Börger, L.; Bennett, D.J.; Choimes, A.; Collen, B. et al. Global effects of land use on local terrestrial biodiversity. Nature; 2015; 520, pp. 45-50. [DOI: https://dx.doi.org/10.1038/nature14324] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25832402]

5. Rozendaal, D.M.A.; Bongers, F.; Aide, T.M.; Alvarez-Dávila, E.; Ascarrunz, N.; Balvanera, P.; Becknell, J.M.; Bentos, T.V.; Brancalion, P.H.S.; Cabral, G.A.L. et al. Biodiversity recovery of Neotropical secondary forests. Sci. Adv.; 2019; 5, eaau3114. [DOI: https://dx.doi.org/10.1126/sciadv.aau3114]

6. Martinez-Garza, C.; Howe, H.F. Restoring tropical diversity: Beating the time tax on species loss. J. Appl. Ecol.; 2003; 40, pp. 423-429. [DOI: https://dx.doi.org/10.1046/j.1365-2664.2003.00819.x]

7. Poorter, L.; Craven, D.; Jakovac, C.C.; van der Sande, M.T.; Amissah, L.; Bongers, F.; Chazdon, R.L.; Farrior, C.E.; Kambach, S.; Meave, J.A. et al. Multidimensional tropical forest recovery. Science; 2021; 374, pp. 1370-1376. [DOI: https://dx.doi.org/10.1126/science.abh3629]

8. DeWalt, S.J.; Maliakal, S.K.; Denslow, J.S. Changes in vegetation structure and composition along a tropical forest chronosequence: Implications for wildlife. For. Ecol. Manag.; 2003; 182, pp. 139-151. [DOI: https://dx.doi.org/10.1016/S0378-1127(03)00029-X]

9. Phillips, H.R.P.; Newbold, T.; Purvis, A. Land-use effects on local biodiversity in tropical forests vary between continents. Biodivers. Conserv.; 2017; 26, pp. 2251-2270. [DOI: https://dx.doi.org/10.1007/s10531-017-1356-2]

10. Olden, J.D.; Rooney, T.P. On defining and quantifying biotic homogenization. Glob. Ecol. Biogeogr.; 2006; 15, pp. 113-120. [DOI: https://dx.doi.org/10.1111/j.1466-822X.2006.00214.x]

11. Baeten, L.; Vangansbeke, P.; Hermy, M.; Peterken, G.; Vanhuyse, K.; Verheyen, K. Distinguishing between turnover and nestedness in the quantification of biotic homogenization. Biodivers. Conserv.; 2012; 21, pp. 1399-1409. [DOI: https://dx.doi.org/10.1007/s10531-012-0251-0]

12. Olden, J.D.; Poff, N.L. Redundancy and the choice of hydrologic indices for characterizing streamflow regimes. River Res. Appl.; 2003; 19, pp. 101-121. [DOI: https://dx.doi.org/10.1002/rra.700]

13. Lewis, S.L.; Maslin, M.A. Defining the Anthropocene. Nature; 2015; 519, pp. 171-180. [DOI: https://dx.doi.org/10.1038/nature14258] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/25762280]

14. Van Der Plas, F.; Manning, P.; Soliveres, S.; Allan, E.; Scherer-Lorenzen, M.; Verheyen, K.; Wirth, C.; Zavala, M.A.; Ampoorter, E.; Baeten, L. et al. Biotic homogenization can decrease landscape-scale forest multifunctionality. Proc. Natl. Acad. Sci. USA; 2016; 113, pp. 3557-3562. [DOI: https://dx.doi.org/10.1073/pnas.1517903113] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/26979952]

15. Baiser, B.; Olden, J.; Record, S.; Lockwood, J.; McKinney, M.L. Pattern and process of biotic homogenization in the New Pangaea. Proc. R. Soc. B Biol. Sci.; 2012; 279, pp. 4772-4777. [DOI: https://dx.doi.org/10.1098/rspb.2012.1651]

16. Smith, C.C.; Healey, J.R.; Berenguer, E.; Young, P.J.; Taylor, B.; Elias, F.; Espírito-Santo, F.; Barlow, J. Old-growth forest loss and secondary forest recovery across Amazonian countries. Environ. Res. Lett.; 2021; 16, 085009. [DOI: https://dx.doi.org/10.1088/1748-9326/ac1701]

17. Martínez-Ramos, M.; Barragán, F.; Mora, F.; Maza-Villalobos, S.; Arreola-Villa, L.F.; Bhaskar, R.; Bongers, F.; Lemus-Herrera, C.; Paz, H.; Martínez-Yrizar, A. et al. Differential ecological filtering across life cycle stages drive old-field succession in a neotropical dry forest. For. Ecol. Manag.; 2020; 482, 118810. [DOI: https://dx.doi.org/10.1016/j.foreco.2020.118810]

18. Chazdon, R.L. Tropical Forest Regeneration☆. Reference Module in Life Sciences; Elsevier: London, UK, 2017.

19. IDEAM—Instituto de Hidrología, Meteorología y Estudios Ambientales. Presentación balance deforestación 2019 (1); IDEAM: Bogotá, Colombia, 2019.

20. Gloor, M.; Barichivich, J.; Ziv, G.; Brienen, R.J.W.; Schongart, J.; Peylin, P.; Cintra, B.B.L.; Feldpausch, T.R.; Phillips, O.L.; Baker, J.W. Recent Amazon climate as background for possible ongoing and future changes of Amazon humid forests. Glob. Biogeochem. Cycles; 2015; 29, pp. 1384-1399. [DOI: https://dx.doi.org/10.1002/2014GB005080]

21. IDEAM—Instituto de Hidrología, Meteorología y Estudios Ambientales. Resultados Del Monitoreo de Deforestación: 1. Año 2020. 2. Primer Trimestre Año 2021; IDEAM: Bogotá, Colombia, 2020.

22. Rodríguez-León, C.H.; Peña-Venegas, C.P.; Sterling, A.; Muñoz-Ramirez, H.; Virguez-Díaz, Y.R. Changes in Soil-Borne Communities of Arbuscular Mycorrhizal Fungi during Natural Regrowth of Abandoned Cattle Pastures Are Indicative of Ecosystem Restoration. Agronomy; 2021; 11, 2468. [DOI: https://dx.doi.org/10.3390/agronomy11122468]

23. Rodríguez-León, C.H.; Peña-Venegas, C.P.; Sterling, A.; Castro, D.; Mahecha-Virguez, L.K.; Virguez-Díaz, Y.R.; Silva-Olaya, A.M. Soil Quality Restoration during the Natural Succession of Abandoned Cattle Pastures in Deforested Landscapes in the Colombian Amazon. Agronomy; 2021; 11, 2484. [DOI: https://dx.doi.org/10.3390/agronomy11122484]

24. Pires, A.P.F.; Srivastava, D.S.; Farjalla, V.F. Is Biodiversity Able to Buffer Ecosystems from Climate Change? What We Know and What We Don’t. BioScience; 2018; 68, pp. 273-280. [DOI: https://dx.doi.org/10.1093/biosci/biy013]

25. Chazdon, R. Second Growth. The Promise of Tropical Forest Regeneration in an Age of Deforestation. Available online: https://www.degruyter.com/document/doi/10.7208/9780226118109/html?lang=en (accessed on 12 June 2022).

26. Ferguson, B.G.; Vandermeer, J.; Morales, H.; Griffith, D.M. Post-Agricultural Succession in El Petén, Guatemala. Conserv. Biol.; 2003; 17, pp. 818-828. [DOI: https://dx.doi.org/10.1046/j.1523-1739.2003.01265.x]

27. Uriarte, M.; Condit, R.; Canham, C.; Hubbell, S.P. A spatially explicit model of sapling growth in a tropical forest: Does the identity of neighbours matter?. J. Ecol.; 2004; 92, pp. 348-360. [DOI: https://dx.doi.org/10.1111/j.0022-0477.2004.00867.x]

28. Comita, L.S.; Thompson, J.; Uriarte, M.; Jonckheere, I.; Canham, C.; Zimmerman, J.K. Interactive effects of land use history and natural disturbance on seedling dynamics in a subtropical forest. Ecol. Appl.; 2010; 20, pp. 1270-1284. [DOI: https://dx.doi.org/10.1890/09-1350.1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/20666249]

29. Instituto Geográfico Agustin Codazzi (IGAC). Caquetá, Características Geográficas; Imprenta Nacional de Colombia: Bogotá, Colombia, 2010.

30. Abbas, S.; Nichol, J.E.; Zhang, J.; Fischer, G.A. The accumulation of species and recovery of species composition along a 70 year succession in a tropical secondary forest. Ecol. Indic.; 2019; 106, 105524. [DOI: https://dx.doi.org/10.1016/j.ecolind.2019.105524]

31. Hu, Y.-K.; Pan, X.; Liu, X.-Y.; Fu, Z.-X.; Zhang, M.-Y. Above-and Belowground Plant Functional Composition Show Similar Changes during Temperate Forest Swamp Succession. Front. Plant Sci.; 2021; 12, 658883. [DOI: https://dx.doi.org/10.3389/fpls.2021.658883]

32. Rocha-Ortega, M.; García-Martínez, M. Importance of Nesting Resources and Soil Conditions for the Recovery of Ant Diversity During Secondary Succession in a Tropical Rainforest. Trop. Conserv. Sci.; 2018; 11, 1940082918787063. [DOI: https://dx.doi.org/10.1177/1940082918787063]

33. Zambiazi, D.C.; Fantini, A.C.; Piotto, D.; Siminski, A.; Vibrans, A.C.; Oller, D.C.; Piazza, G.E.; Peña-Claros, M. Timber stock recovery in a chronosequence of secondary forests in Southern Brazil: Adding value to restored landscapes. For. Ecol. Manag.; 2021; 495, 119352. [DOI: https://dx.doi.org/10.1016/j.foreco.2021.119352]

34. Teixeira, H.M.; Cardoso, I.M.; Bianchi, F.J.; Silva, A.D.C.; Jamme, D.; Peña-Claros, M. Linking vegetation and soil functions during secondary forest succession in the Atlantic forest. For. Ecol. Manag.; 2019; 457, 117696. [DOI: https://dx.doi.org/10.1016/j.foreco.2019.117696]

35. Fick, S.E.; Hijmans, R.J. WorldClim 2: New 1-km spatial resolution climate surfaces for global land areas. Int. J. Climatol.; 2017; 37, pp. 4302-4315. [DOI: https://dx.doi.org/10.1002/joc.5086]

36. Chave, J.; Réjou-Méchain, M.; Búrquez, A.; Chidumayo, E.; Colgan, M.S.; Delitti, W.B.; Duque, A.; Eid, T.; Fearnside, P.M.; Goodman, R.C. et al. Improved allometric models to estimate the aboveground biomass of tropical trees. Glob. Chang. Biol.; 2014; 20, pp. 3177-3190. [DOI: https://dx.doi.org/10.1111/gcb.12629]

37. Velasquez, E.; Lavelle, P.; Andrade, M. GISQ, a multifunctional indicator of soil quality. Soil Biol. Biochem.; 2007; 39, pp. 3066-3080. [DOI: https://dx.doi.org/10.1016/j.soilbio.2007.06.013]

38. McGlinn, D.; Xiao, X.; McGill, B.; May, F.; Engel, T.; Oliver, C.; Blowes, S.; Knight, T.; Purschke, O.; Gotelli, N. et al. Package ‘Mobr’: Measurement of Biodiversity Package Version 2.0.2; The Comprehensive R Archive Network: Vienna, Austria, 2021.

39. Chao, A. Estimating the Population Size for Capture-Recapture Data with Unequal Catchability. Biometrics; 1987; 43, pp. 783-791. [DOI: https://dx.doi.org/10.2307/2531532] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/3427163]

40. Chiarucci, A.; Bacaro, G.; Rocchini, D.; Ricotta, C.; Palmer, M.; Scheiner, S. Spatially constrained rarefaction: Incorporating the autocorrelated structure of biological communities into sample-based rarefaction. Community Ecol.; 2009; 10, pp. 209-214. [DOI: https://dx.doi.org/10.1556/ComEc.10.2009.2.11]

41. Roberts, D.W. Package: “Labdsv”: Ordination and Multivariate Analysis for Ecology Package Version 2.0-1; The Comprehensive R Archive Network: Vienna, Austria, 2019.

42. Oksanen, J.; Blanchet, G.; Friendly, M.; Kindt, R.; Legendre, P.; McGlinn, D.; Minchin, P.; O’Hara, R.; Simpson, G.; Solymos, P. et al. Package ‘Vegan’: Community Ecology Package Version 2.5-7; The Comprehensive R Archive Network: Vienna, Austria, 2018.

43. Anderson, M.J. Distance-Based Tests for Homogeneity of Multivariate Dispersions. Biometrics; 2006; 62, pp. 245-253. [DOI: https://dx.doi.org/10.1111/j.1541-0420.2005.00440.x] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/16542252]

44. Anderson, M.J. A new method for non-parametric multivariate analysis of variance. Austral Ecol.; 2001; 26, pp. 32-46. [DOI: https://dx.doi.org/10.1111/j.1442-9993.2001.01070.pp.x]

45. Martinez-Arbizu, P. Package: “PairwiseAdonis”: Pairwise Multilevel Comparison Using Adonis Package Version: 0.0.1; The Comprehensive R Archive Network: Vienna, Austria, 2017.

46. Legendre, P.; Legendre, L. Numerical Ecology; 3rd ed. Elsevier B.V.: Oxford, UK, 2012.

47. Borcard, D.; Gillet, F.; Legendre, P. Numerical Ecology with R; 2nd ed. Springer International Publishing AG: Cham, Switzerland, 2018.

48. McArdle, B.H.; Anderson, M.J. Fitting Multivariate Models to Community Data: A Comment on Distance-Based Redundancy Analysis. Ecology; 2001; 82, pp. 290-297. [DOI: https://dx.doi.org/10.1890/0012-9658(2001)082[0290:FMMTCD]2.0.CO;2]

49. Borcard, D.; Legendre, P. All-scale spatial analysis of ecological data by means of principal coordinates of neighbour matrices. Ecol. Model.; 2002; 153, pp. 51-68. [DOI: https://dx.doi.org/10.1016/S0304-3800(01)00501-4]

50. Dray, S.; Dufour, A.-B.; Thioulouse, J. Ade4: Analysis of Ecological Data: Exploratory and Euclidean Methods in Environmental Sciences, R Package Version 1.7-16; The Comprehensive R Archive Network: Vienna, Austria, 2020.

51. R Core Team. R: A Language and Environment for Statistical Computing v. 4.0.3; The Comprehensive R Archive Network: Vienna, Austria, 2020.

52. RStudio Team. RStudio v.1.3.1093; RStudi—Open Source & Professional Software for Data Science: Boston, MA, USA, 2020.

53. Di Rienzo, J.A.; Casanoves, F.; Balzarini, M.G.; Gonzalez, L.; Tablada, M.; Robledo, C.W. InfoStat v. 2020; Universidad Nacional de Córdoba: Córdoba, Argentina, 2020.

54. Duivenvoorden, J.F. Vascular plant species counts in the rain forests of the middle Caquetá area, Colombian Amazonia. Biodivers. Conserv.; 1994; 3, pp. 685-715. [DOI: https://dx.doi.org/10.1007/BF00126860]

55. Duque, A.; Sánchez, M.; Cavelier, J.; Duivenvoorden, J.F. Different floristic patterns of woody understorey and canopy plants in Colombian Amazonia. J. Trop. Ecol.; 2002; 18, pp. 499-525. [DOI: https://dx.doi.org/10.1017/S0266467402002341]

56. Olden, J.D. Biotic Homogenization. eLS.; Wiley: Chichester, UK, 2008.

57. Calle, H.; Flórez, J. Así Funciona El Tráfico de Madera en Colombia. Available online: https://es.mongabay.com/2018/10/trafico-de-madera-en-colombia-amazonia-bosques/ (accessed on 14 June 2022).

58. Keith, S.A.; Newton, A.C.; Morecroft, M.D.; Bealey, C.E.; Bullock, J.M. Taxonomic homogenization of woodland plant communities over 70 years. Proc. R. Soc. B Biol. Sci.; 2009; 276, pp. 3539-3544. [DOI: https://dx.doi.org/10.1098/rspb.2009.0938]

59. Daru, B.H.; Davies, T.J.; Willis, C.G.; Meineke, E.K.; Ronk, A.; Zobel, M.; Pärtel, M.; Antonelli, A.; Davis, C.C. Widespread homogenization of plant communities in the Anthropocene. Nat. Commun.; 2021; 12, 6983. [DOI: https://dx.doi.org/10.1038/s41467-021-27186-8] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/34873159]

60. Uchida, K.; Koyanagi, T.F.; Matsumura, T.; Koyama, A. Patterns of plant diversity loss and species turnover resulting from land abandonment and intensification in semi-natural grasslands. J. Environ. Manag.; 2018; 218, pp. 622-629. [DOI: https://dx.doi.org/10.1016/j.jenvman.2018.04.059] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/29715671]

61. Clements, F.E. Plant Succession, An Analysis of the Devel-Opment of Vegetation; Carnegie Institution of Washington: Washington, DC, USA, 1916.

62. Arroyo-Rodríguez, V.; Rös, M.; Escobar, F.; Melo, F.P.L.; Santos, B.A.; Tabarelli, M.; Chazdon, R. Plant β-diversity in fragmented rain forests: Testing floristic homogenization and differentiation hypotheses. J. Ecol.; 2013; 101, pp. 1449-1458. [DOI: https://dx.doi.org/10.1111/1365-2745.12153]

63. Rodríguez, C.; Sterling, A. Sucesión Ecológica y Restauración En Paisajes Fragmentados de La Amazonia Colombiana. Tomo 1. Composición, Estructura y Función En La Sucesión Secundaria; Instituto Amazónico de Investigaciones Científicas-SINCHI: Bogotá, Colombia, 2020.

64. Martínez, L.; Zinck, J. Temporal variation of soil compaction and deterioration of soil quality in pasture areas of Colombian Amazonia. Soil Tillage Res.; 2004; 75, pp. 3-18. [DOI: https://dx.doi.org/10.1016/j.still.2002.12.001]