Abstract

Replanting rubber (Hevea brasiliensis) plantations involves clear-cutting, uprooting and soil preparation, which severely disturb soil biota. This study assessed the effects of legume cover and logging residues on soil macrofauna recovery during the first 24 months after replanting. Macrofauna were monitored following clear-cutting of 38 - 40-year-old rubber stands at two industrial plantation sites in Côte d’Ivoire. Four treatments were compared in randomized block design: bare soil, legume cover (Pueraria phaseoloides) alone, legume cover with leaves, branches and stumps, and legume cover with all residues including trunks. Macrofauna were hand-sorted from soil monoliths (0 - 10 cm) collected in the inter-rows and beneath retained residues, before clear-cutting and at six-months intervals thereafter. Six months after clear-cutting, macrofauna density had declined by 64% and 93%, while taxonomic richness had decreased by 63% and 90% at the flat and sloping sites respectively. However, with only three blocks per site, these declines were statistically significant only at the sloping site. Legume cover alone increased total macrofauna density, mainly through Hymenoptera, whereas taxonomic richness, diversity and the densities of earthworms, isopods, millipedes, centipedes and spiders increased mainly in treatments with residue retention. Thus, residue retention promoted broader community recovery than legume cover alone, although recovery remained incomplete after 24 months: community composition still differed from that of the mature stands, and earthworm density remained below its initial level at the sloping site. No additional effect of retaining trunks was detected during the monitoring period. Macrofauna diversity was positively associated with soil health, carbon transformation and soil structure maintenance in simple regressions, but only its association with carbon transformation remained statistically significant at both sites after accounting for treatment and sampling date. Retaining logging residues in combination with legume cover is therefore recommended to promote soil macrofauna recovery after rubber replanting.

Keywords

Hevea brasiliensis, Pueraria phaseoloides, Clear-Cutting, Soil Health, Biodiversity, Rubber Replanting

Share and Cite:

1. Introduction

Soil fauna is a major component of biodiversity in tropical ecosystems and plays an essential role in ecosystem functioning [1]-[3]. Soil macrofauna contribute to litter decomposition, organic matter transformation, nutrient cycling, and the formation and maintenance of soil structure [4] [5]. Their abundance, diversity and community composition are closely related to the availability of organic resources, vegetation cover, soil properties and microclimatic conditions [6]. Because of their sensitivity to environmental disturbance, macrofauna communities are widely used as indicators of soil quality and of recovery following disturbance [7] [8].

Rubber plantations undergo a major disturbance at renewal. In Côte d’Ivoire, replanting generally involves clear-cutting and mechanical uprooting of the mature trees, followed by windrowing and subsoiling. These operations can alter soil physical and biological properties by increasing bulk density and penetration resistance, reducing porosity and water infiltration, and removing or redistributing the litter and organic matter accumulated during the previous cycle [9]-[13]. The temporary loss of vegetation cover and litter also modifies soil temperature and moisture [14]-[16]. Such changes reduce the trophic resources and habitats available to macrofauna and can strongly affect their abundance, diversity and composition, as reported after logging and replanting in other tropical plantation systems [17] [18].

The recovery of macrofauna after such a disturbance depends, among other factors, on the re-establishment of organic resources and suitable habitats. In rubber plantations, legume cover crops contribute to soil protection, biological nitrogen fixation and organic matter inputs [19] [20]. Retaining the logging residues of the previous cycle (leaves, branches, stumps and trunks left between the planting rows) may further promote biological recovery by supplying organic matter and nutrients, protecting the soil surface and creating favourable habitats [21]-[23]. Organic residue retention is known to improve resource availability, soil organic carbon and soil biological communities [24]-[28]. In the experiment studied here, residue retention accelerated the restoration of soil health indicators [29] and favoured the recovery of springtails, nematodes and microbial communities during the first years after replanting [13] [30] [31].

However, how the whole macrofauna community responds to clear-cutting and mechanical soil preparation, and whether its recovery depends on the legume cover, on the logging residues or on the amount of residues retained, remains unknown. The first two years after replanting are particularly informative because they combine marked changes in soil conditions, re-establishment of vegetation cover and progressive decomposition of residues. We therefore monitored soil macrofauna over 24 months after clear-cutting and replanting at two sites with contrasting soils and topography in Côte d’Ivoire. Our objectives were to 1) quantify the immediate impact of clear-cutting and land preparation on macrofauna density, taxonomic richness and diversity; 2) disentangle the respective effects of a legume cover and of logging residues, and of the amount of residues, on the recovery of these descriptors; 3) characterize changes in community composition and in the main taxonomic groups during recovery; and 4) examine the relationships between macrofauna diversity and soil functions. We hypothesized that (H1) clear-cutting would strongly reduce macrofauna density and diversity; (H2) recovery would be faster with a legume cover than on bare soil, and faster still, increasingly so with the amount of residues retained (R0L0 < R0L1 < R1L1 < R2L1); (H3) community composition in residue treatments would converge toward the pre-felling composition faster than in treatments without residues; and (H4) macrofauna diversity would be positively related to soil functions.

2. Materials and Methods

2.1. Study Sites

The study was conducted in two industrial rubber plantations in southern Côte d’Ivoire, at the experimental sites described by Perron et al. [29] and Kouakou et al. [13] [31]. The Bongo plantation of the Société Africaine de Plantations d’Hévéas (SAPH) is located in the south-east of the country (5˚30'32"N, 3˚32'52"W), and the plantation of the Société des Caoutchoucs de Grand-Béréby (SOGB) in the south-west (4˚43'10"N, 7˚06'42"W). Both sites have a subequatorial climate with a major rainy season from May to July and a minor one from October to November, separated by a long dry season (December-April) and a short one (August-September). Mean annual rainfall is 1640 mm at SAPH and 1674 mm at SOGB.

The SAPH site is almost flat (slopes < 5%) and its soil is a yellow ferralic Arenosol with a loamy sand texture in the 0 - 10 cm layer (10% clay, 2% silt, 87% sand). The SOGB site is hilly (slopes of 10% - 25%) and its soil is a red Ferralsol with a sandy clay loam texture (23% clay, 10% silt, 67% sand). The topsoil (0 - 10 cm) is acidic at both sites (pH-H₂O 4.3 at SAPH and 4.7 at SOGB); total C (0.85% vs 1.28%) and total N (707 vs 974 mg·kg⁻¹) are higher at SOGB, whereas available P (Bray II) is higher at SAPH (7.2 vs 5.2 mg·kg⁻¹) [29]. At both sites the clear-cut stands represented the first rubber cycle after tropical rainforest conversion, and their understory consisted of shrubs and grasses.

2.2. Experimental Design and Treatments

The 40-year-old (SAPH) and 38-year-old (SOGB) rubber stands (253 and 233 living trees·ha⁻¹, respectively) were felled with bulldozers in November 2017. Branches and trunks were cut into 1 - 2 m sections with a chainsaw, and a bulldozer with a V-shaped blade built windrows, so that logging residues were placed in the inter-rows. Between November 2017 and January 2018, the same design was set up at both sites: four treatments in four randomized complete blocks, i.e. 16 plots per site (38.5 × 40 m at SAPH and 40 × 40 m at SOGB; Table 1). In the control (R0L0), all residues were exported, no legume was sown and the plots were kept bare by weekly mechanical weeding. In R0L1, all residues were exported and the legume Pueraria phaseoloides was sown. In R1L1, trunks were exported but leaves, fine branches and stumps were retained together with the legume; this treatment represents a conventional approach to residue management in commercial rubber plantations in Africa [29]. In R2L1, all residues, including trunks, were retained together with the legume.

P. phaseoloides was broadcast-sown at 10 kg·ha⁻¹ in February 2018. From August 2018, the inter-rows were weeded twice a month until the legume outcompeted weeds; it eventually colonized the windrows and completely covered the residues. In early March 2018, all plots were subsoiled to 60 cm depth in the planting rows (not in the inter-rows) and received 200 kg·ha⁻¹ of dolomite. Young rubber trees were planted in May 2018 at SOGB (clone IRCA230, 8 × 2.5 m) and in July 2018 at SAPH (clone IRCA41, 6 × 3 m), with 40 g·tree⁻¹ of a slow-release N-P-K-Mg fertilizer (15-8-9-3) in the planting hole [29].

Table 1. Logging-residue treatments and stocks of carbon and nitrogen in the retained residues (t·ha⁻¹).

Legume = Pueraria phaseoloides. Residue C and N stocks from Kouakou et al. [13].

2.3. Macrofauna Sampling and Identification

Macrofauna were sampled during five campaigns: in October 2017 in the mature stands immediately before clear-cutting (T0), and 6 (April 2018), 12 (October 2018), 18 (April 2019) and 24 months (October 2019) after clear-cutting. Campaigns thus alternated between October (onset of the minor rainy season; T0, 12 and 24 months) and April (end of the long dry season; 6 and 18 months). At T0, before the plots were delimited, one monolith was taken in the inter-row of each of the four blocks. Following clear-cutting, macrofauna were monitored in blocks 1-3, with one monolith collected from the middle of the inter-row of each plot within a 15 × 10 m sampling area. This area was delineated to minimize edge effects and was located close to the sampling points used for soil functional indicators. In plots where residues were retained, samples were collected at the residue-soil interface after carefully moving aside sections of trunks, following the same procedure used for nematode sampling in these plots [13], and the litter and residue fragments within the monolith frame were hand-sorted with the soil. In R1L1 and R2L1, all monoliths were therefore taken beneath the windrows, at the residue–soil interface; no monolith was taken in the parts of the inter-row not covered by residues or in the planting rows, which had been subsoiled. Because the residue-free parts of these plots were not sampled and the area covered by residues was not measured at each date, an area-weighted estimate at the plot scale could not be calculated. The results for the residue treatments therefore describe the macrofauna of the residue-covered inter-row, and the treatment comparisons refer to this habitat. One monolith per plot was taken at 6 months and three monoliths per plot at 12, 18 and 24 months; the three monoliths were averaged per plot. In total, 248 monoliths were processed. The T0 monoliths of blocks 1-3 were used as block-specific pre-disturbance references (n = 3 per site), so that each monitored block served as its own reference; the T0 monolith of block 4 was used only to describe the overall fauna.

Macrofauna were collected following the Tropical Soil Biology and Fertility (TSBF) method [32], modified in depth: a 25 × 25 cm soil monolith was excavated to 10 cm depth and all visible macro-invertebrates were hand-sorted and preserved in 96% ethanol. Sampling was restricted to the 0 - 10 cm layer, where the effects of residue management were expected and were observed in the same plots [29], and to match the depth of the soil functional indicators. In the laboratory, individuals were counted and identified under a stereomicroscope intotaxonomic groups, mostly at the order level.

For each plot and date, density (individuals·m⁻²) was calculated as the mean of the monoliths of the plot, whereas taxonomic richness, Shannon diversity (base-2 logarithm) and Pielou evenness were calculated from the counts pooled across these monoliths, not as the mean of monolith-level indices. Each value was therefore based on one monolith at T0 (one per block) and at 6 months, and on three pooled monoliths per plot at 12, 18 and 24 months. Unlike density, richness and, to a lesser extent, Shannon diversity calculated from pooled counts increase with sampling effort. Richness was therefore also expressed as rarefied richness, i.e., the expected number of groups in a random subsample of 20 individuals [33] [34], calculated from the same pooled counts, to assess whether the difference in sampling effort influenced comparisons with T0. Only samples containing at least 20 individuals qualified for rarefaction; samples with fewer individuals were excluded from this analysis.

2.4. Soil Functional Properties and Soil Health Index

Soil functioning was assessed on the same plots and dates with the Biofunctool^® set of indicators [35] [36], as described in Perron et al. [29] and Kouakou et al. [31]. Carbon transformation was assessed with permanganate-oxidizable carbon (POXC [37]), in-situ basal soil respiration (SituResp^® [38]) and soil fauna feeding activity (bait-lamina test [39]); nutrient cycling with extractable ${\text{NH}}_4^{+}$ and ${\text{NO}}_3^{-}$ (1 M KCl) and ${\text{NO}}_3^{-}$ adsorption on anion exchange membranes [40]; and soil structure maintenance with aggregate stability at 0 - 2 cm and 2 - 10 cm [41], water infiltration (Beerkan method [42]) and the Visual Evaluation of Soil Structure score (VESS [43]). Soil samples were taken from the 0 - 10 cm layer and VESS was assessed on the 0 - 25 cm profile.

Following Obriot et al. [44], Thoumazeau et al. [35] and Perron et al. [29], indicators were converted into unitless scores using “more is better”, “less is better” (VESS) or “optimum” scoring curves. The nutrient indicators used an “optimum” curve, with the T0 concentrations of the mature stands as the optimum [45] [46], which also accounts for the risk of nutrient losses by leaching when availability is excessive [44]. For each function, a principal component analysis on the correlation matrix (FactoMineR [47]) was used to weight the indicators according to their contributions to the components with eigenvalues > 1, yielding three sub-indices (carbon transformation, nutrient cycling and structure maintenance). The sub-indices were given equal weight and summed into a soil health index (SHI). All indices were computed separately for each site.

2.5. Data Analysis

All analyses were conducted separately for each site because the sites differed in soil properties, topography, and pre-clear-cutting macrofauna communities, and because the T0 reference conditions were site-specific. Statistical significance was assessed at P < 0.05.

Immediate disturbance effects

The immediate effect of clear-cutting and land preparation was expressed as the percentage decline in mean density and diversity at 6 months relative to T0, according to the following equation:

$\text{decline}\left(\text{\%}\right)=\frac{\left(\text{T0}-\text{T1}\right)}{\text{T0}}\times 100$

T0: density or diversity value before clear-cutting;

T1: density or diversity value 6 months after the disturbance caused by clear-cutting and the establishment of the treatments.

Because each block had a single T0 monolith, the four plots of a block shared the same reference, and plot-to-baseline differences were not independent. The decline was therefore tested at the block level: the values of the four plots of each block at 6 months were averaged and compared with the T0 monolith of the same block, and the three block-level differences were tested with a paired t-test (n = 3 independent comparisons per site; density log (x + 1)-transformed at the plot level before averaging within blocks). A Wilcoxon signed-rank test was not used because, with three pairs, its smallest attainable two-sided P-value is 0.25. Treatment effects at 6 months were tested with two-way ANOVAs (treatment + block).

Temporal dynamics and treatment effects

For the post-disturbance period (6 - 24 months), density, richness and Shannon diversity were analysed with a split-plot-in-time model, in which treatment was tested against the block × treatment (plot) error term, and date and the treatment × date interaction against the residual error. With this balanced design (3 blocks × 4 treatments × 4 dates), the model is equivalent to a linear mixed model with block as a fixed effect and plot as a random intercept accounting for repeated measurements on the same plots. Density was log (x + 1)-transformed and residuals were checked graphically. Evenness was analysed for 12 - 24 months only, since it is undefined for samples containing fewer than two groups, which were common at 6 months. Treatment means were compared with Tukey’s HSD using the appropriate error term. Changes over time within each treatment, including T0, were tested with two-way ANOVAs (date + block) followed by Tukey’s HSD.

Recovery

Recovery at 24 months was expressed as a recovery ratio, RR (%) = X₂₄/X_T0 × 100, where X₂₄ is the value of a plot at 24 months and X_T0 that of the T0 monolith of the same block. RR = 100% indicates full recovery to the initial level, whereas RR > 100% indicates a recovery exceeding the initial level. RR was compared among treatments with two-way ANOVAs (treatment + block) and Tukey’s HSD, and compared with 100% with one-sample t-tests.

Community composition

Densities were square-root transformed and Bray–Curtis dissimilarities computed. Composition over 12 - 24 months (a balanced dataset without empty samples) was analysed with permutational multivariate analysis of variance (PERMANOVA [48]) including date, treatment and their interaction, first for each site and then for both sites together with site as an additional factor. Because the same plots were sampled at each date, the permutation scheme followed the split-plot (repeated-measures) structure of the design, with 9999 permutations. The treatment effect was tested by permuting whole plots, i.e. moving the complete time series of each plot together, within blocks, and its pseudo-F ratio used the plot (block × treatment) mean square as denominator. Because the four treatments can be arranged in only 4! = 24 ways within each of the three blocks, all 13,824 possible whole-plot permutations were enumerated for this test. The date and date × treatment effects were tested by permuting dates within each plot, so that observations were exchanged only within a plot’s time series, against the residual mean square. In the analysis combining both sites, site was treated as an additional whole-plot factor and plots were nested in blocks within sites. Treatment and site × treatment were tested against the plot mean square by permuting whole plots within blocks, date and its interactions by permuting dates within plots, and site against the block-within-site mean square by permuting complete blocks between sites; because only 20 such arrangements exist with three blocks per site, the smallest attainable P-value for site is 0.10. With three plots per treatment, a whole-plot permutation test of a pairwise contrast between two treatments has only 2³ = 8 distinct permutations within blocks, so that its smallest attainable P-value is 0.125. Pairwise contrasts were therefore described by their partial R² and pseudo-F and were not used for formal inference; the univariate analysis of dissimilarity to T0 (see below) provided the treatment comparisons. Homogeneity of multivariate dispersions among treatments was tested with PERMDISP [49], and the groups contributing most to the dissimilarity between treatments were identified with SIMPER [50]. Compositional recovery was quantified, for each plot and sampling date (6 - 24 months), as the Bray-Curtis dissimilarity to the T0 community of the same block (an empty sample having a dissimilarity of 1). The resulting values were analysed using the split-plot model described above. Temporal trajectories were visualized by non-metric multidimensional scaling (NMDS) of the Bray-Curtis dissimilarities between treatment × date mean compositions (Kruskal stress-1 reported). The densities of the nine main groups over 12 - 24 months were analysed with the split-plot model after log(x + 1) transformation.

Links between macrofauna diversity and soil health

The relationships between Shannon diversity and the soil health index (SHI) and its three sub-indices over 6 - 24 months were analysed for each site (n = 48 plot × date observations: 3 blocks × 4 treatments × 4 dates) with linear mixed models in which each index was the response variable, Shannon diversity a continuous predictor, treatment, sampling campaign and block fixed effects, and plot a random intercept accounting for the repeated measurements of the same plots. The diversity coefficient therefore estimates the association between diversity and soil function within treatments and campaigns, i.e. beyond their common response to residue management and season. Models were fitted by restricted maximum likelihood (R packages lme4 and lmerTest), and the diversity coefficient was tested with the Kenward–Roger approximation of the denominator degrees of freedom; adjusted coefficients are reported with their 95% confidence intervals. Simple linear regressions, which ignore these dependencies, are also shown for illustration only, with their adjusted R² values.

All statistical analyses were performed using R 4.2.0 [51]. Univariate, diversity and multivariate analyses, including PERMANOVA, were conducted using appropriate R packages. The split-plot permutation tests of the PERMANOVA and the block-level paired t-tests were computed in Python 3.13 (NumPy, SciPy). Soil indices and their regressions with macrofauna diversity were computed using FactoMineR and ggplot2 [52].

3. Results

3.1. Macrofauna Community Characteristics

In total, 27,368 individuals belonging to 20 taxonomic groups were collected with 11,009 individuals at SAPH and 16,359 at SOGB. Hymenoptera (46.7% of individuals), Isoptera (20.7%), earthworms (13.3%) and Coleoptera (6.3%) accounted for 87% of the individuals, whereas 11 groups each represented less than 1%. Before clear-cutting, mean density was 981 ind.·m⁻² at SAPH and 2389 ind.·m⁻² at SOGB, with an average of 8.3 and 6.7 taxonomic groups, respectively (Table 2). Community composition differed between site: at SAPH, Hymenoptera (44% of individuals) and earthworms (29%) were dominant, whereas earthworms (60%), Isoptera (17%) and Hymenoptera (16%) predominated at SOGB.

3.2. Effect of Clear-Cutting and Land Preparation

Six months after clear-cutting, density had declined by 64% at SAPH and 93% at SOGB, taxonomic richness by 63% and 90%, and Shannon diversity by 54% and 91%, respectively (paired t-tests on block means, n = 3, P = 0.009 - 0.036 at SOGB and 0.059 - 0.32 at SAPH; Table 2). The decline in density was highly variable at SAPH, where two plots of block 1 (R1L1 and R2L1) and the R0L0 plot of block 2 had higher densities than the T0 value of their block, so that the mean density of block 1 exceeded its T0 value (648 vs 448 ind.·m⁻²), whereas no macrofauna were found in any plot of block 3; at the block level, none of the declines at SAPH was therefore significant. At this date, treatments did not differ significantly at either site (two-way ANOVAs, P ≥ 0.09 for density, richness and Shannon diversity).

Table 2. Soil macrofauna parameters before (T0) and 6 months after clear-cutting and land preparation (all treatments pooled), and percentage decline of the site mean.

Values are means ± SE. P: paired t-test on the three block-level differences between the mean of the four plots at 6 months and the T0 monolith of the same block (blocks 1-3; density log (x + 1)-transformed at the plot level before averaging within blocks).

3.3. Recovery of Density and Diversity after Legume and Residue Inputs

Over 6 - 24 months, treatment had a significant effect on density, richness and Shannon diversity at both sites, and date on all three descriptors; the treatment × date interaction was significant only for richness at SOGB (Table 3). Averaged over this period, density was much lower in R0L0 than in the three legume treatments, which did not differ from each other at either site (back-transformed means: SAPH, 121 vs 399 - 950 ind.·m⁻²; SOGB, 53 vs 306 - 740 ind.·m⁻²; Tukey, P < 0.05). In contrast, richness and diversity distinguished the residue treatments from the legume-only treatment. At SAPH, richness was 3.3 groups in R0L0, 6.3 in R0L1 and 8.9 - 9.3 in R1L1 and R2L1 (three significantly different levels), and Shannon diversity was higher in R1L1 (1.92) than in R0L0 (0.98) and R0L1 (1.21), with R2L1 intermediate (1.73). At SOGB, richness was higher in R1L1 and R2L1 (7.5 - 7.8 groups) than in R0L0 (3.0), with R0L1 intermediate (5.3), and Shannon diversity was higher in R1L1 and R2L1 (1.62) than in R0L0 and R0L1 (0.93 - 0.99). R1L1 and R2L1 never differed significantly.

Within treatments (Figure 1), richness in R0L0 remained significantly lower than T0 at all sampling dates at SAPH (3.3 - 3.7 vs 8.3 groups), while both density and richness in R0L0 at SOGB showed only partial recovery. In R0L1, R1L1 and R2L1, density and richness no longer differed significantly from T0 from 12 months onwards at both sites; at SOGB, richness in R1L1 and R2L1 at 18 months even exceeded T0 (11.7 and 10.3 vs 6.7 groups). At SAPH, Shannon diversity in R0L1 tended to decline again after 12 months (0.94 at 24 months vs 2.07 at T0), whereas it did not differ significantly from T0 in the residue treatments. Evenness over 12 - 24 months differed among treatments at SOGB (F_3,6 = 10.18, P = 0.009) but not at SAPH (P = 0.14); at both sites it was lowest in R0L1 (0.46 - 0.48 vs 0.64 - 0.67 in R0L0 and R1L1), where Hymenoptera made up 77-83% of individuals at 24 months.

Table 3. Effects of treatment, date and their interaction on macrofauna parameters over 6 - 24 months after clear-cutting (split-plot-in-time model; P values in brackets).

At 24 months (Table 4), density in the three legume treatments had returned to or exceeded the T0 level (recovery ratio of 115% - 356%), whereas it remained at 67% (SAPH) and 9% (SOGB) of T0 in R0L0. Recovery ratios for richness were highest with residues (112% - 146%) and lowest in R0L0 (44% - 71%). For Shannon diversity, recovery ratios in the residue treatments (74% - 154%) did not differ significantly from 100%, whereas R0L1 remained at 45% (SAPH; significantly below 100%) and 76% (SOGB) of T0. Rarefied richness showed that richness ratios above 100% partly reflected the larger sampling effort from 12 months onwards: at 24 months, rarefied richness in R1L1 and R2L1 (3.8 - 5.4 groups per 20 individuals) did not differ significantly from T0 (5.7 at SAPH and 3.9 at SOGB), whereas R0L0 at SAPH remained lower than T0 (2.6 groups; Tukey, P < 0.05).

Table 4. Recovery ratio (%) of macrofauna density, taxonomic richness and Shannon diversity 24 months after clear-cutting.

Values are means ± SE (n = 3). Different letters indicate significant differences among treatments within a row (ANOVA with block, Tukey’s HSD, P < 0.05); no letters: treatment effect not significant. * Mean significantly different from 100% (one-sample t-test, P < 0.05).

Figure 1. Temporal changes in (A) taxonomic richness, (B) total density, (C) Shannon diversity and (D) Pielou evenness of soil macrofauna in each treatment at SAPH and SOGB. Month 0 = before clear-cutting (T0, shared by all treatments; dotted line). Points are means ± SE (n = 3). Different letters indicate significant differences among dates within a treatment (ANOVA with block, Tukey’s HSD, P < 0.05; density tested after log(x + 1) transformation); ns = no significant date effect. Evenness is shown only for samples with at least two groups and was not tested within treatments. R0L0 = no legume, no residues (bare soil); R0L1 = legume (Pueraria phaseoloides) only; R1L1 = legume + leaves, fine branches and stumps; R2L1 = R1L1 + trunks.

3.4. Changes in Community Composition

In the NMDS (Figure 2(A), Figure 2(B)), all treatments moved away from the T0 composition at 6 months. From 12 months onwards, R1L1 and R2L1 formed a compact group, closest to T0 at SOGB, whereas R0L0 followed an erratic trajectory far from T0 and R0L1 was intermediate. Over 12 - 24 months, treatment was the main driver of community composition (PERMANOVA, R² = 0.50 at SAPH and 0.45 at SOGB, P = 0.007 at both sites), whereas date explained less variation (R² = 0.08, P = 0.002 at SAPH; R² = 0.05, P = 0.048 at SOGB); the date × treatment interaction was significant at SAPH (R² = 0.16, P = 0.011) and at SOGB (P = 0.028) (Table 5; whole-plot permutations within blocks for treatment, permutations within plots for date). Pairwise contrasts between treatments explained 3% - 51% of the variation (partial R²), the smallest values being found for R1L1 vs R2L1 (3.5% at SAPH and 2.9% at SOGB). In the analysis combining both sites, treatment explained 44% of the variation (F = 20.34, P = 0.0001) and site only 3% (F = 4.17, P = 0.10, the smallest attainable value with three blocks per site), and the site × treatment interaction was not significant (P = 0.66), indicating a consistent treatment effect at both sites. Multivariate dispersion differed among treatments at both sites (PERMDISP, F = 15.8 at SAPH and 17.9 at SOGB, P < 0.001) and was highest in R0L0 (mean distance to centroid 0.36 - 0.41 vs 0.19 - 0.28 in the other treatments), so the contrast between R0L0 and the other treatments also reflects a larger heterogeneity among bare plots.

The Bray-Curtis dissimilarity to T0 decreased after 6 months at both sites (date effect: SAPH, F_3,24 = 4.44, P = 0.013; SOGB, F_3,24 = 42.2, P < 0.001) and was lower in the three legume treatments than in R0L0 (treatment effect: SAPH, F_3,6 = 12.46, P = 0.006; SOGB, F_3,6 = 19.96, P = 0.002; Figure 2(C), Figure 2(D)). At 24 months, it was 0.69 (SAPH) and 0.72 (SOGB) in R0L0, 0.46 - 0.48 in the other treatments at SAPH, and 0.51, 0.40 and 0.33 in R0L1, R2L1 and R1L1 at SOGB, where R1L1 and R2L1 were significantly closer to T0 than R0L0 (Tukey, P < 0.05).

Table 5. PERMANOVA of macrofauna community composition (Bray-Curtis dissimilarities on square-root-transformed densities) over 12 - 24 months after clear-cutting at each site.

Split-plot PERMANOVA with 9999 permutations: treatment tested against the plot (B × T) term by permuting whole plots within blocks (complete enumeration of the 13,824 possible permutations); date and D × T tested against the residual term by permuting dates within plots. With this balanced design, the sums of squares of T, D and D × T do not depend on the order of terms; n = 36 plot × date samples per site.

Figure 2. Changes in soil macrofauna community composition after clear-cutting. (A, B) NMDS of Bray-Curtis dissimilarities between the mean compositions of each treatment × date at SAPH and SOGB; lines join T0 (yellow circle, centred at the origin) to successive dates (6 - 24 months) for each treatment. At SOGB, no individuals were found in R1L1 at 6 months. (C, D) Bray-Curtis dissimilarity between each plot and the T0 monolith of its block (mean ± SE, n = 3; 1 = no taxon in common). R0L0 = bare soil; R0L1 = legume only; R1L1 = legume + leaves, fine branches and stumps; R2L1 = R1L1 + trunks.

3.5. Responses of the Main Taxonomic Groups

SIMPER showed that Hymenoptera, earthworms and Isoptera contributed most to the dissimilarity between R0L0 and the residue treatments (together 52% at SAPH and 57% at SOGB), followed by Coleoptera, Scolopendromorpha and Diplopoda at SAPH, and by Isopoda, Diplopoda and Araneae at SOGB. Over 12-24 months, Hymenoptera were 9 to 30 times more abundant in the three legume treatments than in R0L0, with no difference between R0L1 and the residue treatments (Table 6). In contrast, earthworms, Isopoda, Diplopoda and Scolopendromorpha were abundant only where residues were retained. At SAPH, earthworm density was 311 - 337 ind.·m⁻² in R1L1 and R2L1, i.e. similar to T0 (283 ind. m⁻²), versus 39 in R0L1 and 7 in R0L0. At SOGB, it reached 412 - 689 ind.·m⁻² with residues versus 118 in R0L1 and 23 in R0L0, but remained well below the T0 level (1424 ind.·m⁻²). Isopoda, absent from the T0 samples, reached 39 - 214 ind.·m⁻² in the residue plots. Araneae were more abundant in all legume treatments than in R0L0, and most abundant with residues at SOGB. Isoptera were most abundant in the residue treatments at SOGB (991 - 1439 ind.·m⁻²) and made up 48% of the individuals in R2L1 at 24 months, which largely explains the high density of this treatment.

Table 6. Density (ind.·m⁻²) of the main macrofauna groups before clear-cutting (T0) and in each treatment over 12 - 24 months after clear-cutting.

Treatment values are means over 12, 18 and 24 months (n = 9 plot × date observations). F and P: treatment effect in the split-plot model on log(x + 1)-transformed densities. Different letters indicate significant differences among treatments (Tukey’s HSD, P < 0.05); letters are omitted when the treatment effect was not significant, and identical letters despite a significant F indicate that no pairwise contrast was significant after adjustment.

3.6. Relationships between Macrofauna Diversity and Soil Functions

Figure 3. Relationships between soil macrofauna diversity (Shannon index) and (A) the soil health index, (B) carbon transformation, (C) soil structure maintenance and (D) nutrient cycling at SAPH and SOGB over 6–24 months after clear-cutting (n = 48 plot × date observations per site). Lines are simple linear regressions shown for illustration; values are their adjusted R² ∗∗ P < 0.01; ∗∗∗ P < 0.001 for these simple regressions, which do not account for treatment, campaign or repeated measurements of plots (adjusted mixed-model estimates are given in the text). Colours indicate treatments.

In simple regressions, macrofauna Shannon diversity was positively related to the soil health index at both SAPH (R² = 0.31) and SOGB (R² = 0.29; Figure 3(A)). It was also positively related to the carbon transformation index (SAPH, R² = 0.38; SOGB, R² = 0.48; Figure 3(B)) and to the structure maintenance index (SAPH, R² = 0.16; SOGB, R² = 0.15; Figure 3(C)). In contrast, no association was detected between diversity and the nutrient cycling index at either site (adjusted R² ≈ 0; Figure 3(D)). After adjusting for treatment, sampling campaign and block while accounting for repeated measurements within plots, the association between diversity and the soil health index remained positive and statistically significant at SOGB (β = 0.052, 95% CI 0.009 to 0.096, P = 0.020), but was no longer significant at SAPH (β = 0.001, 95% CI −0.027 to 0.029, P = 0.93). Carbon transformation was the only sub-index showing a positively and statistically significant association with diversity at both sites (SAPH, β = 0.020, 95% CI 0.0004 to 0.039, P = 0.046; SOGB, β = 0.031, 95% CI 0.009 to 0.052, P = 0.006), although the lower confidence limit at SAPH was close to zero. The association with structure maintenance was not statistically significant at either site (SAPH, β = −0.007, 95% CI −0.021 to 0.007, P = 0.34; SOGB, β = 0.012, 95% CI −0.004 to 0.028, P = 0.13). Similarly, no statistically significant association was detected between diversity and nutrient cycling (SAPH, β = −0.012, 95% CI −0.032 to 0.009, P = 0.25; SOGB, β = 0.009, 95% CI −0.016 to 0.034, P = 0.46). For the soil health index and the carbon transformation, the adjusted coefficients were smaller than the slopes of the simple regressions (soil health index: 0.001 vs 0.059 at SAPH and 0.052 vs 0.063 at SOGB; carbon transformation: 0.020 vs 0.044 at SAPH and 0.031 vs 0.046 at SOGB). This reduction was particularly pronounced at SAPH. Plot-level variance was estimated to be zero or close to zero in most models. Consequently, the adjusted estimates were close to those obtained from the corresponding fixed-effects linear models.

4. Discussion

4.1. Clear-Cutting and Land Preparation Strongly Reduced Macrofauna

Clear-cutting and land preparation reduced macrofauna density by about two-thirds at SAPH and by more than 90% at SOGB, with similar losses of richness and diversity, and many monoliths were empty six months after clear-cutting. This supports H1 and agrees with the decline in soil fauna reported after mechanized logging and replanting in other tropical and forest systems [17] [18]. With three blocks, however, the decline was statistically significant only at SOGB; at SAPH it varied strongly among blocks, from an increase in block 1 to the absence of macrofauna in block 3. It also parallels the responses of other soil organisms in the same plots, where total nematode abundance fell by 87% (SAPH) and 78% (SOGB) [13] and all soil functions except structure maintenance at SOGB were significantly disturbed [29]. The likely causes are the loss of habitat and organic resources, soil compaction and disturbance by heavy machinery [9] [53]-[55], and the change in microclimate caused by the loss of canopy and understory, which exposes the soil to high radiation and temperature and low moisture [56] [57]. Treatments did not differ at 6 months, as observed for soil microbial communities in the same experiment [31]: at that time the legume had been sown only two months earlier and residue decomposition had barely started.

Two features of the experimental design should be considered when interpreting the magnitude of this decline. First, the 6-month sampling took place at the end of the long dry season, one month after subsoiling, whereas T0 was sampled in October, so part of the decline may be seasonal. The observation that both density and richness in R0L0 remained significantly lower than in T0 at 12 and 24 months, during the same season as T0, provides further evidence of a genuine effect of the disturbance. Second, the stronger decline at SOGB may be related to its steeper slopes, which favour the loss of topsoil and organic matter by erosion [58], to the more intense compaction by machinery reported for this more clayey soil [13], and to a pre-felling community dominated by earthworms (60% of individuals), which are particularly sensitive to soil disturbance. With only two sites, however, these factors cannot be separated from each other or from other site differences.

4.2. Legume Cover Restored Density, Logging Residues Restored Diversity

Macrofauna recovery followed two distinct patterns across treatments. Total density was restored by legume cover alone, as R0L1 reached or exceeded the T0 level at both sites, similar to the residue treatments R1L1 and R2L1. This recovery was driven mainly by Hymenoptera, which accounted for 77% to 83% of individuals in R0L1 at 24 months and contributed to its low evenness. Social insects are highly mobile and rapid colonizers that can exploit the shade, litter and roots provided by P. phaseoloides, whose litter releases half of its carbon within a week [59]. In contrast, taxonomic richness, Shannon diversity and most other groups, including earthworms, Isopoda, Diplopoda, Scolopendromorpha and, at SOGB, Araneae, showed greater recovery in the residue treatments R1L1 and R2L1 than in R0L1. Thus, legume cover was sufficient to restore total density, whereas the retention of logging residues was associated with a broader, though incomplete, recovery of macrofauna diversity and composition. These saprophagous groups depend directly on dead wood and litter as food and shelter, and their predators follow them [4] [60]. Logging residues also buffer soil temperature and moisture and, in these plots, drastically reduced runoff and soil loss [23], which probably favoured the survival and recolonization of less mobile, desiccation-sensitive taxa. Hypothesis H2 is therefore only partly supported: recovery was faster with the legume than on bare soil and faster still with residues, but this ranking applied to diversity rather than to density. It is the contrast between R1L1 or R2L1 and R0L1, rather than between these treatments and R0L0, that isolates the effect of residue retention, because R0L0 represents bare soil subjected to weekly weeding and does not reflect a realistic management practice.

No additional effect of the amount of residues was detected: R1L1 and R2L1 did not differ significantly in density, diversity or composition. With three blocks, however, the power to detect moderate differences between these two treatments was limited, and this absence of a significant difference does not show that they are equivalent. The same result was obtained for soil functions in these plots over 18 months [29]. Trunks, which distinguish R2L1 from R1L1, decompose much more slowly than leaves and twigs (about 29 months vs 6 - 12 months in rubber plantations [61]), so their effect on macrofauna may appear later. The large density of Isoptera in R2L1 at SOGB suggests that trunks are already being exploited by wood-feeding termites. This is consistent with the importance of the quality of organic inputs, and not only their quantity, for macrofauna [62] [63].

Recovery nevertheless remained incomplete after 24 months. The density overshoot observed in the residue treatments, particularly at SOGB, reflects an increase in the abundance of a few social insect groups rather than a return to the T0 community. Earthworms, the dominant group before clear-cutting at SOGB, remained at 29% - 48% of their T0 density even with residues. Moreover, richness values above T0 were largely an effect of the larger sampling effort from 12 months onwards, since rarefied richness in the residue treatments did not exceed T0. The absence of residues did not completely prevent recovery, as Shannon diversity in R0L0 at SOGB returned close to its initial value. However, R0L0 communities remained were impoverished and heterogeneous. The removal of logging residues is known to reduce the density and diversity of several groups of soil organisms by removing habitat and altering understory vegetation and nutrient availability [64].

4.3. Community Composition during Recovery

Treatment explained about half of the variation in community composition at 12-24 months at both sites, far more than date, and its effect did not differ significantly between sites despite their contrasting soils. This consistency contrasts with the strong site dependence of nematode and microbial responses in the same experiment [13] [31] and suggests that, for macrofauna, the management of organic inputs overrides the effect of soil texture at this stage. As predicted by H3, the residue treatments showed the lowest dissimilarity from T0, particularly at SOGB, whereas R0L0 remained distant. However, dissimilarity to T0 remained 0.33 - 0.48 in the residue treatments after 24 months, and the post-disturbance communities included taxa absent before clear-cutting (T0), such as Isopoda, while earthworms remained under-represented at SOGB. Recovery therefore involved a partial reassembly towards a new community rather than a simple return to the mature-plantation state. The larger multivariate dispersion in R0L0 indicates that bare plots developed idiosyncratic communities, possibly because macrofauna are spatially aggregated in scarce carbon-rich hotspots such as roots and organic debris [65].

4.4. Macrofauna Diversity and Soil Functioning

Macrofauna diversity was positively related to the soil health index, carbon transformation and soil structure maintenance at both sites in simple regressions; after accounting for treatment, campaign and repeated measurements, only the association with carbon transformation persisted at both sites, and that with the soil health index persisted only at SOGB. These results partly support H4 and are consistent with the relationship between macrofauna richness and soil health reported in the same experiment over the first 18 months [29] and with the recognized link between soil biodiversity and ecosystem functions [66] [67]. Macrofauna activity contributes to aggregation and macroporosity [4] [68], which is consistent with the relationship with structure maintenance observed in the simple regressions, although this relationship did not persist after adjustment. No relationship was found with nutrient cycling, which in this study was assessed through mineral N availability and mainly depends on microbial processes [69].

These relationships are correlative and should be interpreted with caution. The residue treatments cluster at high diversity and high function values, so both variables may respond jointly to logging residues inputs rather than one driving the other. The mixed models addressed part of this issue by estimating the associations within treatments and campaigns while accounting for the repeated observations of the same plots; the adjusted associations were weaker than those of the simple regressions, and the association with the soil health index was no longer detected at SAPH. Much of the covariation between diversity and soil functions thus reflects their common response to residue management and season. The data therefore do not allow us to conclude that a loss of diversity causes a loss of function.

4.5. Limitations and Management Implications

The primary limitations of this study are outlined as follows: 1) the investigation was conducted at only two locations, which hinders the ability to attribute site-specific differences to particular factors; 2) monitoring involved three blocks per site post-clear-cutting, which reduces the statistical power, especially when comparing treatments at specific times or testing the immediate decline against T0; 3) there was a disparity in sampling effort between the initial time point (T0) and 6 months, where one monolith per plot was used, compared to 12 to 24 months, where three monoliths were utilized. This discrepancy may affect estimates of taxonomic richness but not density, and was assessed through richness rarefaction analysis; 4) the taxonomic resolution was relatively coarse, potentially limiting the sensitivity of diversity indices and obscuring some functional responses; 5) sampling was confined to the 0 to 10 cm soil layer, which might result in underestimating groups residing deeper, such as endogeic earthworms; 6) the variation in sampling seasons introduces complexity to the date effect; 7) there was no parallel monitoring of an undisturbed mature stand; 8) the study’s monitoring period of 24 months is short compared to the time required for trunk decomposition; and 8) in the residue treatments, sampling was restricted to the residue-covered inter-row, so that the results do not provide area-weighted estimates at the plot scale. Future research should aim to extend monitoring throughout the immature phase, incorporate more detailed taxonomic identification and functional group analysis, and sample deeper soil layers.

From a management perspective, sowing a legume cover as early as possible after clear-cutting is sufficient to restore macrofauna abundance, whereas retaining logging residues between the planting rows promoted a broader, though still incomplete, recovery of macrofauna diversity and of decomposer and predator groups in the residue-covered inter-row. No additional benefit of retaining trunks (R2L1) was detected over the common practice of retaining leaves, fine branches and stumps while exporting trunks (R1L1) during the 24 months of monitoring. This suggests that trunks may be valorized without compromising the early recovery of macrofauna, although the two practices cannot be considered equivalent on the basis of three blocks and a later effect of trunks cannot be excluded as they decompose; their long-term contribution to soil fertility also remains to be assessed [22] [29]. The adoption of residue retention should also consider operational and sanitary constraints, including access for planting and maintenance, fire risk and the potential role of stumps and roots as inoculum sources for root rot diseases such as white root disease [70].

5. Conclusion

Clear-cutting and mechanical soil preparation during rubber replanting caused a severe loss of soil macrofauna, with declines in density and richness of about two-thirds at SAPH and more than 90% at SOGB; with three blocks, this decline was statistically significant only at SOGB. Over the 24 months of monitoring, a legume cover was sufficient to restore macrofauna density, largely through social insects, whereas the retention of logging residues promoted a broader recovery of taxonomic richness, diversity and saprophagous and predatory groups such as earthworms, isopods, millipedes, centipedes and spiders, and brought community composition closest to that of the mature plantation. This recovery remained incomplete after two years: community composition still differed from the initial community, and earthworms remained well below their initial density at SOGB. No additional effect of retaining trunks with leaves, fine branches and stumps was detected during the monitoring period, which does not exclude a later effect as trunks decompose. In the residue treatments, these conclusions apply to the residue-covered inter-row where macrofauna were sampled. Macrofauna diversity was positively associated with soil health, carbon transformation and soil structure in simple regressions, but only the association with carbon transformation persisted at both sites after accounting for treatment and sampling date, although these relationships remain correlative. These results support combining logging-residue retention with a legume cover to promote the recovery of soil macrofauna during rubber replanting.

Acknowledgements

We are grateful to the managers of the SAPH and SOGB plantations for providing facilities and labour during the sampling campaigns, and to Daouda Kouassi and Antoine Manizan (SOGB) and Soumahoro Mouman, Angelica N’chot, Richmond Konan and Jean-Marc Toto (SAPH) for their help in supervising the fieldwork.

Funding

This study was carried out within the FERTIM project, funded by the Institut Français du Caoutchouc (IFC) and by the Michelin, SIPH and SOCFIN companies. A.K. Kouakou received additional support from the IRD PhD grant programme (ARTS, 2019).

Author Contributions

Kouakou Aymard: Conceptualization, Methodology, Data curation, Formal analysis, Investigation, Visualization, Writing original draft, Writing, review & editing. ZON Serge: Formal analysis, Writing, review & editing. Silué Simon: Methodology, Visualization, Writing, review & editing. Alain Brauman: Methodology, Conceptualization, Validation, Writing, review & editing, Supervision, Project administration.

Conflicts of Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

References