Introduction
Tropical montane forests, specifically high Andean ones, entail one of the most important ecosystems on Earth due to their high biodiversity [1, 2] and their key role providing vital ecosystem services to human populations: water, food, and climate regulation. However, high montane forests are one of the most threatened ecosystems and, thus, are target of conservation and restoration efforts promoted in recent years [3-7]. The Andean region of Colombia is the most densely populated of the country. This has led to extensive transformations of its forests. In consequence, only 10 % of the original Andean forests are still conserved and less than 5 % of these are specifically high montane forests [2, 5, 8]. For over a century, Colombian high Andean forests have been transformed into exotic tree plantations. The exotic pine Pinus patula is one of the tree species most used for this purpose [4, 9, 10].
Compared to areas with semi-natural high montane forest vegetation, the soil environment in pine plantations experiences unique dynamics. Pine needles deposited on the soil are rich in lignin and take longer to degrade. This results in the accumulation of a layer of scarcely degraded organic matter (OM). This type of OM is characteristic of coniferous forests [11-14]. Moreover, litter produced by species of Andean forests has a high N and C content, thus its degradation and nutrient turnover rate is accelerated [15-18]. In addition, the low temperature and humidity of the Andean region decreases the speed of OM degradation in the surface horizon. This leads to accumulation of OM in the top soil [19-23]. OM is the basis and the main source of energy of the soil's food chain [24-27], therefore, any change in OM forms will have strong effects on the entire soil biota, mainly on its macrofauna, because their digestive system lacks of lignolytic enzymes to process material with highly polymerized polysaccharides [28-32].
Edaphic macrofauna entail all soil invertebrates larger than 2 mm belonging to different taxonomic groups. Macrofauna are essential elements of the edaphic ecosystem because of their size, whereby they have a unique role in the ecosystem, specifically by regulating soil dynamics such as nutrient cycling and controlling pests and diseases [33-36]. Moreover, they have an indirect effect on OM mineralization and humification processes by enhancing functions of the edaphic meso and microfauna [37-39]. This link with ecosystem functioning makes them indicators of soil quality [40, 41], disturbance [24, 42, 43] and successional processes [44-46]. For these reasons, in the last decade, there has been a growing interest in studying their role in ecological restoration processes [38, 47-50].
In Colombia, studies have investigated the macrofauna composition of Andean forests, exotic tree species plantations, and agricultural lands, revealing a higher soil macrofauna diversity in Andean forests compared to both conifer plantations and agricultural areas [51-59]. However, there are no studies investigating macrofauna succession processes nor the relationships between macrofauna and physico-chemical variables in areas after clearcutting of P. patula. In this study we addressed the questions on how do soil macrofauna composition, structure and function vary among areas with different post clearcutting ages of Pinus patula, and which physico-chemical soil variables better explain the possible variation. To answer these questions, we evaluated the composition and structure of the edaphic macrofauna of three high montane successional areas and a high Andean forest patch. Our leading hypothesis was that areas that experienced pine clearcuttings earlier, and have no recent input of needles onto their soils, were more suitable for the arrival of different macrofauna groups. Consequently, we expected that the similarity with the Andean forest patch will increase with the time elapsed after clearcutting.
Materials and methods
Study site
The survey was conducted at the Parque Forestal Embalse del Neusa (PFEN), a forest reserve established around an artificial lake in the Cundinamarca department, located at 3 100 masl on the eastern Cordillera of Colombia. The average annual temperature in this forest reserve is 10.5 ◦ C, rainfall regime is bimodal, and average precipitation is 1 025 mm/year. The study site is humid (40-70), according to the De Martonne climatic classification method. Soils are classified by their textural class as Silty Clay Loam [60], with predominance of silt and clay. This kind of soil is characterized by its high water retention capacity, high organic matter content, and drainage problems [61]. In 1950, the P. patula plantation was established around an artificial reservoir in the PFEN to avoid any sedimentation and erosion processes that could affect the waterbody [62]. However, because the plantations never received silvicultural management, they caused several environmental problems; consequently, in 2009 the environmental authority Corporación Autónoma Regional de Cundinamarca decided to carry out the removal of these plantations, and to launch an ecological restoration program.
Three areas with different post-clearcutting ages (years after clearcutting Yac) were selected to perform the samplings: Laureles (0 Yac), Chapinero (2.5 Yac) and Guanquica (5 Yac). In addition, a high montane Andean forest patch was selected as the reference forest (RF) (Fig. 1). Since the studied areas are of different sizes, the smallest one (0 Yac 12.6 ha) was chosen to standardize the grid dimension. Then, through satellite imaging, a grid of 300 m2 and 316 points was located, each point was separated by a distance of 20 m (Suppl. 1A). The same grid dimension was used to locate the points in the other areas (2.5 Yac, 5 Yac and RF) also by using satellite images (Suppl. 1B-D). After grids were located, 12 points were randomly selected in each area to establish the experimental plots (10 m2).
The 0 Yac area is almost lacking any vegetation, and only a few pioneer species are present: Phytolacca bogotensis and Muehlenbeckia tamnifolia. The soil in this 0 Yac area is mostly covered by P. patula leftovers (needles, stumps, branches and some logs). Vegetation in the 2.5 Yac area is dominated by the shrub Rubus floribundus and the herbaceous plant Holcus lanatus, the latter is also the dominant species in the 5 Yac area; in this area there is little accumulation of litter on the ground, which consists mainly of herbaceous leaves. The RF is a remnant of tropical high montane forest, with a perhumid regime. It is dominated by Weinmannia tomentosa, however, other species such as Chusquea scandens, Drimys granadensis and Myrsine coriacea are common in this ecosystem.
Soil macrofauna sampling
Sampling of edaphic macrofauna was carried out in July 2014, using the tropical soil biology and fertility method (TSBF), which consists of extracting 25 x 25 x 20 cm soil monoliths [41, 63]. In total, 12 monoliths per area and 48 in total were obtained. Each monolith was stored in a plastic bag to avoid drying. Later, organisms were manually extracted and preserved in 96 % alcohol. Oligochaetes and larvae were previously fixed in 5 % formaldehyde for 72 hours. Indicator morphospecies were identified to the narrowest possible taxonomic level according to existing identification guides [64][65][66][67][68][69][70].
Soil sampling and physico-chemical measurements
Inside each 10 m2 plot of every study area (0 Yac, 2.5 Yac, 5 Yac, RF), a soil sample composed of 3 subsamples was randomly extracted by using the cylinder method [71]. Subsequently, 1 kg of each composite sample was stored in a plastic bag and transported to the Instituto Geográfico Agustín Codazzi [72]. The physical variables evaluated were: 1) slope; 2) bulk density, which was determined by extracting unchanged soil samples at soil depth of 20 cm by using the paraffin cylinder method [73]; 3) real density that was measured by the pycnometer method [74]; 4) texture, which was determined through the Bouyoucos method by using a hydrometer [73, 75]; 5) porosity, that was determined by the mercury porosimetry technique [76]; and 6) temperature, which was measured with a soil thermometer in the monolith.
The chemical variables evaluated were: 1) pH, measured with a potentiometer, for which a water solution (1:1) was prepared [76]; 2) interchangeable acidity, and 3) saturation percentage of interchangeable acidity (S.I.A) were determined by using 1 M potassium chloride extraction and quantification by volumetric titration [76]; 4) soil organic carbon (SOC) was determined by the Walkley-Black volumetric titration [77]; 5) total nitrogen was obtained by Kjeldahl potentiometric titration [76]; 6) cation exchange capacity (CEC) [78]; 7) percentage base saturation (BS) and 8) total bases Ca, Mg, K, Na (BT) were determined by 1 M ammonium acetate extraction (pH 7) and quantification by volumetric titration [76]; and 9) available phosphorous (P) was determined by using the extraction method with Bray II solution and by spectrophotometric quantification in the visible range [76] (Table 1).
Statistical analysis
To characterize each edaphic community, the richness of orders and their abundance, as the number of individuals, were quantified. Orders with fewer than three individuals were excluded from the statistical analysis (Lithobiomorpha, Psocoptera, Mantodea and Pseudoscorpiones). The alpha diversity was measured using the Shannon-Wiener index and species evenness was calculated with the Pielou's evenness index, as implemented in the R-package Vegan [67]. Macrofauna data were tested for normality with the Shapiro-Wilk test and showed a significant deviation from normality (Shapiro-Wilk < 0.05). Therefore, the Kruskal-Wallis test (H : Chi square value; P : significance value; n: sample size) was used to compare the Shannon and Pielou indices among groups. Pairwise comparisons were used to identify significant differences between areas. Moreover, in order to quantify the compositional dissimilarity between areas the Bray-Curtis dissimilarity index and Ward method were used; complete-linkage was implemented as clustering technique. The selection of indicator species in each zone was carried out by using the indicator value index (Ind.Val) [79], using the R-package Labdsv [80]. All statistical analysis were conducted in R [81].
Due to pseudoreplication (lack of spatial replicates per site) among the plots from each area, a non-parametric nested manova [82] was performed using the R-package BiodiversityR [83] with 100 000 random permutations. Because there was not a significant effect of plots on macrofauna abundance (F = 0.00124; P > 0.05), by using generalized linear models (GLMs) for multivariate abundance data [84], the effect of soil physico-chemical variables was evaluated on morphospecies abundance in each of the four areas. GLMs were fitted using the R-package Mvabund [85]. In order to prevent multicollinearity, an individual multicollinearity diagnostic was performed by using the variance inflation factor (VIF) with the R-package mctest [86], as well as a correlation analysis on soil physico-chemical data. pH, interchangeable acidity (A.I), soil organic carbon (CO), nitrogen (N), cation exchange capacity (CEC), Calcium (Ca), Magnesium (Mg), Potassium (K), Sodium (Na), and percentual base saturation (BT). Clusters of variables showing collinearity and high correlation were discarded from the regression analysis. Then, a single variable from each cluster was included in the models. However, soil pH showed a weak correlation with Na and K, therefore, all three variables from this cluster were included in the analysis. GLMs, with the Gaussian error distributions, included the explanatory variables: slope, temperature, bulk density, real density, loam, pH, P, Na and K.
Results and discussion
We collected a total of 3 591 individuals and 24 orders, among which Diptera, Polydesmida, Coleoptera and Isopoda were the most abundant (Table 2). Abundance was significantly lower in the 0 Yac site than in the other sites (RF, 2.5 and 5 Yac) (K-W, H = 11; P = 0.01; n = 48) (Table 3, Fig. 2). This was an expected result because saprophagous groups from the macrofauna show difficulty processing material with high lignin content [28-32]. Consequently, richness and diversity were lowest in the 0 Yac (K-W, H = 28.3; P < 0.00001; n = 48; K-W, H = 32.4; P < 0.00001; n = 48, respectively). In addition, 90 % of the individuals recorded in this area were Dipteran (family Sciaridae) larvae, and 9 % to Coleoptera larvae of species Ancognatha ustulata and A. scarabaeioides. These insects occurr in this area because of their dispersal capacity; as despite their larvae having little mobility, adults have flight capacity [87, 88]. Neita-Moreno & Morón [89] reported that Ancognatha larvae feed mainly on organic matter, but also can feed on vegetation; this versatility in their feeding sources allows Ancognatha to thrive in 0 Yac. Moreover, these larvae were the largest bodied individuals collected, and according to Kalinkat et al. [90] the largest soil organisms are more diet-flexible. These findings agree with those by [91], who recorded that during the first years after clearcutting, saprophagous specialists were most apt to colonize these areas. Overall, few organisms were collected in the 0 Yac, which explains why, according to the indicator value index (Ind.Val), no single species can be regarded as significant indicator of this area (P > 0.1).
Species richness and Shannon index were significantly higher in the forest than in the other sites (K-W, H = 28.3; P < 0.00001; n = 48; K-W, H = 32.4; P < 0.00001; n = 48, respectively). These results agree with other studies in Colombia [51-59]. In contrast, Pielou's index (P) did not show any significant differences between areas with several post clearcutting ages (K-W, H = 6.35; P = 0.09; n = 48), but it was higher and less variable in the forest than the other sites (Table 3). Species richness and diversity were higher in areas with 2.5 and 5 Yac than in 0 Yac, and richness was higher in 2.5 Yac than in the 5 Yac. However, neither abundance nor richness changed significantly between areas (K-W, H = 0.01; P = 0.91; n = 48; K-W, H = 0.04; P = 0.83; n = 48, respectively). Moreover, the 2.5 Yac showed the highest diversity values after the reference forest. The low diversity recorded for the 5 Yac is probably related to the dominance of two morphospecies, which were the best indicators for this area according to the Ind.Val (Table 4).
According to the Bray-Curtis dissimilarity index, 2.5 and 5 Yac were the most similar to the forest, at the same time these two areas showed the highest similarity between them (Fig. 3). The similarity between 2.5 Yac with the forest is probably related to the OM forms, litter quality, and the transition from an initial OM Mor to a Moder-mull humus type [14]; the latter caused by the decomposition of the needles and the addition of litter from different plant species. After pine clearcutting, the accumulation of needles on the ground stops, and the litter layer-composed mainly of needles-is consumed or processed at a higher rate than any litter layer produced by pioneer plant species, facilitating an overall and gradual litter reduction. Also, in open areas, the litter degradation process advances faster than in the native forest because of direct exposure to solar radiation and rainfall. Consequently, the maintenance of the litter and fermentation layers (OL, OF), typical of the native forest, is prevented. As observed in the 5 Yac area, a change to OM without fermentation layer and with a thin litter but with a high humic layer (OH) occurred. On the other hand, the 0 Yac was the most dissimilar to the 2.5 and 5 Yac areas. Although we did not measure the lignin content of the studied areas, we deduced, based on previous studies [92][93], that fresh pine needles contained a higher amount of lignin in comparison to those that have been degrading for more than two years. This finding agrees with the observation that during the first successional stages after clearcutting, the environmental conditions are constantly changing, producing a very stressful environment for colonizing organisms; therefore, during these stages, species turnover is very slow [94-96].
In all four studed areas, the environmental variables with significant effect on morphotype abundance were slope, temperature, bulk density, real density, loam, pH, P, Na and K (GLM; P < 0.05 in all cases) (Table 5). Remarkably, areas with slopes greater than 11° showed the richest macrofauna (Table 2, Fig. 2). This is because slope is one of the factors regulating the amount of OM and humidity of the soil [97-100]; in this study the highest slope values were recorded in the RF and the 2.5 Yac areas. Temperature also had a direct effect on detritivorous fauna, in turn it could have affected OM transformation, soil moisture, nutrient retention [93, 98, 101-103] and the microbiological activity [104].
Macrofaunas were affected by several soil chemical features (Table 5). Specifically, the studied macrofauna showed a positive response to the soil's Na and K content. K content has also been correlated with higher decomposition rates [96, 105]. This, in turn, evidences the importance of organic matter layer thickness. In the organic fermentation horizon of the litter layer, H3O+ ions are released as a result of the hydrolysis of organic acids, explaining the high acidity of the forest [19, 82]. These ions regulate the solubility of elements such as P and S, which have been recognized as limiting elements for the Andean vegetation [1, 108, 109]. In addition, P, Na and pH have been previously reported to determine the presence of some macrofauna groups, as Araneae, Coleoptera, Haplotaxida and Isopoda [110, 111]. Therefore, the formation and depletion of the litter layer are key processes to consider in the restoration of areas where clearcutting of exotic species has been carried out.
The importance of the litter layer for the restoration of high Andean forests was recognized through the identification of the morphospecies showing the highest fidelity and specificity to RF. The presence of a symphylan (Ind.Val = 0.83) and a polydesmid of the genus Phaneromerium (Ind.Val = 1) in the RF soil ecosystem reflects the selectivity of these two species for the RF and highlights the importance of detritivorous fauna in this ecosystem. The genus Phaneromerium is native of South America, its species inhabit exclusively the epigean zone and are characterized by their small size (< 6 mm) [65]. This can be associated with their high population densities in this study. However, two predator morphospecies belonging to the Schendylops and Parajapyx genera were identified as the third and fourth best indicators of the forest (Ind.Val = 0.75 and 0.74, respectively) (Table 4). Overall, the orders most associated with RF were Polydesmida, Symphila, Diplura, Geophilomorpha, Scolopendromorpha and Glommeridesmida (Table 2).
For the 2.5 Yac, the best indicator was a diplopod (Batodesmus sp.) (35 % of the individuals recorded in this area were Dipolopods). A native worm (Glossodriulus palenke) and a cockroach of the Blattidae family were the second and third best indicators ofthis area. All these groups are saprophagous, so they have an important effect on OM fragmentation and soil structure [ 112, 113]. The high occurrence of detritivore groups such as Diplopoda, Isopoda and Oligochaeta in this area could be related to the fact that recalcitrant and woody materials enter their own successional process. As the degradation of these materials proceeds, they are more accessible to a larger number of saprophagous groups [30, 114]. Accordingly, León etal. [18] recorded that the life time of woody material in the soil and the litter produced by P. patula in a low montane forest ranges from 2.1 to 2.5 years. Consequently, the arrival or population explosion of macrofauna detritivorous groups could be associated with the time it takes to the remaining material of the P. patula plantation to reach its degradation point. During this time, the degrading material is accessible to local saprophagous species. This is an important factor to consider for future restoration plans in these areas, because this finding could be interpreted as the first functional change recorded during the early succession of these areas. Also, the increase in saprophagous organisms could indicate the beginning of the heterotrophic phase, characterized by a change from moder type to mull type humus, and an excess in the mineralization rate in comparison to photosynthesis rate [115]. Saprophagous organisms can be enhancers of successional processes due to their role on nutrient cycling and soil formation [93, 116, 117] and on impacting higher levels of the edaphic food chain, promoting post-disturbance recovery of vegetation [118, 119]. The importance of diplopods for the Colombian Andean ecosystems had been previously highlighted [120]. Accordingly, the management of pine stumps, logs and needles should be considered for monitoring restoration and conservation projects in this area; the importance of the management of these plant materials, during the early stages of succession, for the edaphic fauna has been previously stated [94, 121]
For the 5 Yac area, the best indicator was a weevil species of the genus Amphideritus (Coleoptera: Curculionidae), which was more abundant in the larval state, Ind.Val = 0.82 (Table 4). This morphospecies and one belonging to the genus Conoderus (Coleoptera: Elateridae), which was the second-best indicator of this area, were almost absent in the other sites. The third indicator was a spider of the genus Lycosa. These three genera have been associated with grasslands and areas with high anthropogenic intervention [68, 122-126]. Moreover, these morphospecies were not recorded in the forest, and because 5 Yac is the most distant to RF, it is considered that surrounding areas have an important effect on the successional macrofauna process.
Soil compaction is one of the consequences of clearcutting [ 19]. Consequently, the abundance of groups that spend most of their life cycle under the soil could be reduced. Groups that depend on their ability to dig shelters to avoid dehydration, include Geophilomorpha [64, 127] and Diplopoda [128]. In this study, Geophilomorpha, was almost restricted to the forest in which we recorded the highest proportion of loam and the lowest bulk and real density. However, Diplopoda (Polydesmida) were more associated to the 2.5 Yac area (456 individuals), where a higher bulk and real density were registered. This is in concordance with the significant increase in macrofauna abundance as these variables increase (Table 5). Possibly, these variables are more related to the mineral fraction of the soil, and diplopods using the layer of mixed litter (deposited pine needles and litter of the local vegetation) as refuge and feeding substrate. Our results support this statement because the number of Polydesmids recorded in the 5 Yac area was lower (34 individuals) than in the 2.5 Yac area. Even thougth, both areas presented similar bulk and real densities. Moreover, a decrease in the abundance ofthis group would indicate a decrease in the thickness of the OM layer (heterotrophic phase). These results partially agree with findings of a decline in the diplopod community in areas with 5 years after clearcutting [115]. However, this study did not evaluate areas with less than 5 years after clearcutting, therefore it is not possible to establish if the diplopod community experienced a population explosion as we recorded two years after clearcutting.
Our results suggest that passive restoration of areas where clearcutting of P. patula has been carried out is an unviable option, and that at 5 years after clearcut, soil degradation processes may persist or ecological succession may be delayed. These two scenarios could result from vegetation cover removal, direct exposure to sun, wind, and rain, and the arrival of exotic grasses. Accordingly, increases in bacterial and fungal metabolism would result in OL and OF layer reductions while increasing the OH layer, as we observed in the 5 Yac [19, 103, 129]. Furthermore, the isolation of 2.5 and 5 Yac, i.e. the absence of nearby forest patches as propagule sources, is another factor adding to the unviability of passive restoration. The kind of management of leftover P. patula plantations and surrounding areas as well as its effects on macrofauna colonization; ought to be evalueated in future studies.
Overall, our findings support previous studies [95, 96, 130, 131], concluding that edaphic macrofauna are more linked to changes in vegetation cover and to quality and quantity of organic matter. Moreover, evidence fom the present study shows that physico-chemical soil features related to OM have an effect on its macrofauna, and that areas with low OL and OF layers can facilitate the establishment of grasslands [107, 108, 132]. The arrival of grasses may be interpreted as an intermediate stage in the ecological succession; however, changes in the spatial location of carbon stocks and evidence about grassland settlement at the beginning of Andean ecosystems degradation , suggest that this hypothesis is unlikely [43, 97, 107, 131, 132]. Finally, our results support the hypothesis that clearcutting changes the composition of species the function of the ecosystem, leading to random successional trajectories [91]. In our study, these trajectories are determined by the biological offer of nearby areas, and the degradation rate of organic matter. This particularly affects soil macrofauna groups with little mobility, such as Diplopoda and Geophilomorpha, as well as for small-bodied species with specific nutritional diets.
Conclusions
In comparison to the reference forest, the macrofauna composition changed significantly among areas where clearcutting of exotic species has been carried out. Slope, temperature, bulk density, real density, loam, pH, P, Na and K were the physico-chemical soil variables that better explained the observed soil macrofauna differences. These results also suggest that the maintenance of the OM layer should be taken into account when designing future restoration and conservation plans in areas where clearcutting of exotic species has been implemented. On the other hand, these experimental data did not confirm the hypothesis that the similarity with the native forest patch increases with the age after clearcutting. Moreover, since soil macrofauna indicator species (Ind.Val) differed between study areas, the hypothesis of random successional trajectories may be favoured [91]. However, it is necessary to carry out additional samplings to investigate whether the initial soil communities carry on founder effects that determine their unique long-term assemblies [135]. Alternatively, the observed differences in species composition and structure could be transient, occuring at the onset of successional processes, and ultimately leading to assembly convergence with the time.
This research opens new perspectives for the restoration of high montane forests, showing that the pine needle layer is a relevant element of the ecosystem during the stages following P. patula clearcutting. Moreover, this study highlights the importance of different soil macrofauna orders to edaphic processes that will later have an effect on the vegetation. Under the conditions of this study, community assembly analyses based on taxonomic information at the order level proved to be useful to find differences between areas during the early stages of succession; however, analyses carried out at genus or species levels could reveal more detailed patterns. In addition, these results show that groups such as Araneae, Diplopoda, Coleoptera, and even Chilopoda might be useful to monitor and evaluate restoration processes in these areas. Furthermore, counting on species-level information within these orders can help improve the monitoring and diagnosis of these montane forest areas.