Soil Health Changes Over a 25-Year Chronosequence From Forest to Plantations in Rubber Tree (Hevea brasiliensis) Landscapes in Southern Côte d'Ivoire: Do Earthworms Play a Role?

The agro-ecological drawbacks of the spread of rubber tree plantations in Côte d’Ivoire since the 1990’s are obvious even though they have not been properly investigated. They consist of biodiversity loss, land degradation and food insecurity, which have extended into the existing cocoa-led degraded areas whose rehabilitation have unfortunately not started. This situation increases not only the threat on soil health status but also undermines the capability of soils to deliver ecosystem services that are key to sustainable agricultural production. The current study took advantage of a chronosequence in rubber tree landscapes to assess soil health deterioration in general and possibly earthworm-mediated role in soil health changes. The hypothesis underpinning this study was that earthworms contribute to mitigate soil health deterioration in rubber-dominated landscapes due to their key role in soil functioning. This study conﬁrmed that the conversion of forest to rubber tree plantations signiﬁcantly impaired all soil biological, physical, and chemical parameters at the beginning (7 years) of the chronosequence; followed further by a restorative trend taking place beneath the plantations from 12 years. However, this study failed to ﬁnd evidence of a direct role of earthworms in soil health rehabilitation over time. Mesoscale studies along with the use of appropriate models could help unravel this “black box” and shed some light on the contribution of earthworms as key soil ecosystem engineers.

The agro-ecological drawbacks of the spread of rubber tree plantations in Côte d'Ivoire since the 1990's are obvious even though they have not been properly investigated. They consist of biodiversity loss, land degradation and food insecurity, which have extended into the existing cocoa-led degraded areas whose rehabilitation have unfortunately not started. This situation increases not only the threat on soil health status but also undermines the capability of soils to deliver ecosystem services that are key to sustainable agricultural production. The current study took advantage of a chronosequence in rubber tree landscapes to assess soil health deterioration in general and possibly earthworm-mediated role in soil health changes. The hypothesis underpinning this study was that earthworms contribute to mitigate soil health deterioration in rubber-dominated landscapes due to their key role in soil functioning. This study confirmed that the conversion of forest to rubber tree plantations significantly impaired all soil biological, physical, and chemical parameters at the beginning (7 years) of the chronosequence; followed further by a restorative trend taking place beneath the plantations from 12 years. However, this study failed to find evidence of a direct role of earthworms in soil health rehabilitation over time. Mesoscale studies along with the use of appropriate models could help unravel this "black box" and shed some light on the contribution of earthworms as key soil ecosystem engineers.
Keywords: biodiversity, earthworms, functional groups, land use change, soil degradation, soil threats, rubber tree plantations

INTRODUCTION
In the last two decades, replacement of degraded lands by rubber tree plantations has become common in the humid and subhumid areas of Côte d'Ivoire. The overwhelming presence of these new tree plantations in the agro-ecological landscapes is due to its huge economic returns in lieu of the low profitability of cocoa due to falling world market prices (Ruf, 2012). As a result, smallholder farmers seeking livelihood improvement have entered the scene by (i) replacing their old cocoa or coffee plantations with rubber trees, and (ii) converting the remaining portions of secondary forests and fallow lands into rubber stands. This has significantly contributed to an over 99% rise in the area and production of rubber tree plantations at the national level since 1960 (FAOSTAT, 2016). The immediate outcome is the rise of Côte d'Ivoire to the top as Africa's leading rubber producer with a total production share of 45.9% (FAOSTAT, 2016) The total area and production are estimated at 189,937 ha and 310,655 tons, respectively. It is well-documented that the conversion of natural ecosystems to rubber tree stands resulted in agro-ecological and environmental drawbacks in areas of the world where the cultivation is possible. They include loss of biodiversity (Pia and Konrad, 2015), deterioration of soil quality and reduction of soil organic carbon stocks (Oku et al., 2012;de Blécourt et al., 2013), acute soil erosion, disruption of streams, and risk of landslides (Fox et al., 2014). In other words, as a land-use system characterized by sequential changes in space and time, monoculture rubber tree farming has the potential to increase the magnitude of threats to soil over time (Günal et al., 2015). Soil threats due to land use change and exploitation are summarized in 10 main groups out of which soil erosion, soil organic carbon (SOC) deterioration, nutrient imbalance, loss of soil biodiversity, and soil compaction are the most important (FAO and ITPS, 2015).
It is now well-acknowledged that soils are self-organized ecological systems within which organisms (microorganisms, predators organized in micro food webs, and ecosystems engineers) interact in a nested suite of discrete scales (Lavelle, 1997;Lavelle et al., 2016). Soil ecosystem engineers composed mostly of earthworms, termites, and ants play key roles in creating habitats for other organisms and controlling their activities through physical and biochemical processes (Lavelle, 1997;Jouquet et al., 2006). Furthermore, they contribute to deliver ecosystem services through three different processes (Puga-Freitas et al., 2012;Puga-Freitas and Blouin, 2015;Lavelle et al., 2016): (i) the organization of soil physical structure and associated ecosystem services, (ii) the selection and activation of plant, microbial and smaller invertebrate communities that determine decomposition and nutrient cycling processes, and (iii) the release of hormones that regulate primary production.
All terrestrial ecosystems consist of aboveground and belowground components that interact to influence community and ecosystem-level processes and properties (Wardle et al., 2004). This is particularly true in rubber tree landscapes presumably characterized by the successional replacement of land use and land cover that drive above and belowground interactions. This change also undermines the deliverance of soil-based ecosystem services which are driven by soil organisms among which soil macroinvertebrates play a critical role like any other land uses around the world (Spurgeon et al., 2013;Franco et al., 2016;de Valençia et al., 2017).
One way of assessing the reaction of soil systems to the perturbations brought about by rubber tree plantations is to assess the health of soil in these derived landscapes. Such an assessment should be integrated and holistic including physical, chemical, and biological components. The latter comprise key organisms and community structure with a special reference to their interactive feedback with abiotic parameters. Earthworms are key soil macroinvertebrates due to their contribution to the functioning of ecosystems (Blanchart et al., 1999(Blanchart et al., , 2004Blouin et al., 2013;Fischer et al., 2014) and plant productivity (van Groenigen et al., 2014;Xiao et al., 2018) which are very well-documented. These findings mostly from studies carried out in controlled conditions (laboratory, mesocosms, etc.) are indicative of the potential interactive feedback between earthworms and soil processes. Furthermore, semi-natural studies in mesocosms have revealed that tropical earthworms are organized in ecological groups composed of detritivores and geophageous polyhumics, mesohumics and oligohumics according to their feeding behavior (Lavelle, 1981). There are also two main groups, compacting, and decompacting, based on their physical impact on soils (Blanchart et al., 1999(Blanchart et al., , 2004. To date, studies have rarely attempted to investigate the extent to which interactive feedbacks between earthworms and soil processes influence soil health in field conditions. The successional stages of different land use along a chronosequence in rubber tree landscapes offer the framework to such a study. This study aims to investigate earthworm-mediated role in soil health changes beneath rubber tree plantations as compared to baseline forests. The hypothesis underpinning the current study is that earthworms help mitigate soil health deterioration in rubber-dominated landscapes through their functional impact on soil chemical and physical characteristics.

Study Site
The study was conducted in the village Tiéviessou (Latitude: 5 • 08 ′ 13 ′′ N; longitude: 5 • 01 ′ 26 ′′ W; elevation: 6 m) located in Grand-Lahou District, southwestern Côte d'Ivoire. This region is characterized by a bimodal humid tropical climate marked by two rainy seasons and two dry seasons with steady significant seasonal variation in the past two decades. The total annual rainfall and average temperature of the study year (2013) were 1085.35 mm and 26.9 • C, respectively. The study site is in the Guinean domain and belongs to the ombrophilous area characterized by dense evergreen forests (Guillaumet and Adjanohoun, 1971). The area has experienced tremendous human pressure leading to human-derived landscapes made of degraded secondary forests, plantations of oil palm, rubber tree and cocoa along with food crops. However, a portion of the landscape has been reserved as a protected area known as the Classified Forest of Gobodienou (Figure 1). This landscape is highly irrigated due to the presence

Sampling Design
To meet the objective of the study, sampling plots were selected to capture the most representative features of the rubber tree (Hevea brasiliensis) landscapes together with portions of the baseline ecosystem along a chronosequence: food crops (cassava) of 1 to 3-year-old; 7, 12, and 25-year-old smallholder plantations around Tiéviéssou village and the two settlements, Agnouanssou and Betesso (Table 1). Smallholder famers who are owners of these plantations use very few inorganic inputs in general.
A survey conducted at the onset of this study revealed that they only used inorganic fertilizers including urea and NPK as inputs in the early stages of the plantations to help trees grow smoothly. The chronosequence that is a set of rubber tree plantations in the landscape that share similar attributes but with different ages, was found to be most relevant approach in line with the objective of this study as the initial date of the disturbance and subsequent history of the sites were known (Walker et al., 2010;Zhou et al., 2017). The general assumption supporting chronosequence-based study uses natural forests as baseline ecosystems from which all man-made systems were at first derived. In most cases, the cycle of smallholder plantations starts with food crops that are grown in association with planted rubber trees which will form the basis of monospecific plantations 3 years later. The 7-year-old plantations represent the initial stage of production of the plantation, the 12-yearold ones are halfway of their productive life and the 25-year-old plantations are considered as fully mature because the complete production cycle of a plantation can reach 40 years. In addition to plantations, secondary forests selected from the protected area were considered as baseline ecosystems. By a stratified sampling approach, three sampling plots of size 10 m × 50 m each, were used as replicates and selected afterwards in each land use type (LUT) such that the total number of plots amounted to 15. In each plot, five sampling points of which one was placed at the center point and the remaining at the four corners, were established. Hence, in total, 15 plots were selected along the chronosequence (forest, food crops, rubber 7, rubber 12, and 25 years) and geo-referenced using GPS. Soil health indicators that are sensitive to land use changes and related to some key soil functions such as nutrient cycling, soil structure regulation, soil biodiversity conservation (Schulte et al., 2015) were measured. Chemical parameters namely pH-H 2 O (pH), soil organic carbon (SOC) and macronutrients including total nitrogen (N), total phosphorus (P) and exchangeable potassium (K) were measured. Physical parameters, namely bulk density, aggregate size distribution, and stability were measured. Ecological metrics (density, biomass, species richness, and diversity) and community structure of earthworms were considered as proxies of soil organisms.
Field campaigns to collect data pertaining to soil health indicators were conducted from mid-September to mid-November 2013, corresponding to the short rainy season, which is the most suitable period for earthworm collection.

Soil Sampling and Chemical Analyses
In each plot, soil samples were collected from five distinct points: one at the center, four points at 2 m and arranged in the 4 cardinal points. The samples were collected from topsoil (0-20 cm) using an auger. They were thereafter pooled and thoroughly mixed as a single composite sample. Each composite sample was used as replicates and thus totaling 15 samples per LUT, giving 75 (5 × 15, as we have 5 LUT) in the study area. Samples were air-dried for a week and homogenized using a 2-mm mesh sieve. An aliquot of 100 g of the fine fraction was used and analyzed for pH-H 2 O (pH), soil organic carbon (SOC) and nutrients including total nitrogen (N), total phosphorus (P), exchangeable potassium (K) determination. SOC concentration was measured using the modified method of Anne (Nelson and Sommers, 1982) while N content was extracted by the method of Nelson and Sommers (1980) and determined using the Technicon autoanalyzer (Technicon Industrial Systems, 1977). Phosphorus (P) was measured by colorimetry following nitriperchloric acid digestion and subsequent molybdenum-blue color development (Olsen and Sommers, 1982). Potassium (K) was extracted using ammonium acetate buffer (pH 7) and determined by means of atomic absorption spectrophotometry techniques (Anderson and Ingram, 1993).

Bulk Density
Soil samples were collected using the cylinder in layers 0-10 cm and 10-20 cm. The bulk density (BD) was calculated at laboratory depending on the inner diameter of the core sampler, sampling depth and the oven dried weight at 105 • C. Soil water content was measured gravimetrically and expressed as a percentage of soil water to dry soil weight.

Aggregate Size Distribution and Mean Weight Diameter (MWD)
Soil aggregate distribution was determined using the dry-sieving method (Gilot, 1994) consisting of soil cores collection from each LUT using a cylinder in 0-10 cm and 10-20 cm layers and air-dried to a moisture content of about 5% of the dry weight. The total mass was weighed and identified aggregates further broken by dropping the dry soil blocks from a constant height (1 m) onto a hard surface. Subsequently, the samples were successively placed on a set of six stacking sieves of different meshes (50, 100, 200, 500, 1,000 and 2,000 mµ in ascending order resulting in six aggregate size fractions: <50 mµ, 50-100 mµ 100-200 mµ, 200-500 mµ, 500-1,000 mµ, 1,000-2,000 mµ and FIGURE 3 | Proportional changes in earthworm functional groups' density along the LUTs (forest, food crops, 7, 12, and 25-year old rubber plantations). G, Geophageous. >2,000 mµ). All fractions were weighed, and data analyzed to compute the proportion of aggregate distribution and the mean weight diameter (MWD). MWD is used to measure soil structural stability (Ge et al., 2018). This index is calculated as follows: -MWD is the mean weight diameter (mm), n is the number of aggregate fractions (five), x i is the mean diameter of the ith fraction w i is the weight of soil in the fraction i expressed as a percentage of the dry soil mass.

Sampling and Identification of Earthworms
Earthworms were sampled using a modified method recommended for tropical soils (Anderson and Ingram, 1993), which involves digging of 5 monoliths (25 × 25 × 30 cm) along a transect stretching across the sampling plot.
Since we were not concerned with the vertical distribution of earthworms, the modified size and depth of the monolith of this sampling scheme was used in this study. A soil monolith (50 × 50 × 20 cm) was dug out at each sampling point and used as replicates in each plot with 5 replicates in total per plot and thus 15 per LUT. The monolith was surrounded by a trench of 20 cm depth preventing earthworms from escaping. The sampled soil blocks were deposited in a bucket and specimen were collected by hand-sorting in trays according to the layers (0-10 cm and 10-20 cm) and were preserved in 4% formaldehyde solution. Earthworm specimen were identified to species level (Tondoh and Lavelle, 2005;Csuzdi and Tondoh, 2007), counted, weighed and further allocated into four functional groups. The most common accepted ecological classification is a division in three groups (Bouché, 1977); anecics (vertical burrowers), epigeics (litter layer/surface inhabitants), and endogeics (mineral soil inhabitants). Later on, based on their feeding behavior and ecological characteristics (Lavelle, 1981), provided a nomenclature fitting the tropical context as follows: detritivores and geophageous polyhumics, mesohumics and oligohumics. Detritivores are litter feeders, which feed at or near the soil surface on plant litter. Geophageous earthworms feed deeper in the soil and derive their nutrition from soil organic matter and dead roots ingested with mineral soil (Lee, 1985). Owing to their dependence on soil organic matter (Lavelle, 1981), geophageous earthworms were further divided into three groups: the polyhumics, which feed on decaying residues mixed with little mineral soil, the mesohumics, which feed on soil fairly rich in organic matter, and the oligohumics, which feed on organic matter-poor soil.

DATA PROCESSING AND ANALYSIS
To characterize earthworm community structure in the different LUTs, a Correspondence Analysis (CA) was performed on the matrix of species abundance per LUT using the FactoMineR and factoextra packages in R statistical software. The CA uses a simple Chi square statistic to test for significant dependence between species and LUTs. It also provides the factor scores for species and LUTs, which we used to represent their association graphically. The distance between species is indicative of their similarity (or dissimilarity) in the LUTs space and was thus used to associate each species to each LUT. For each LUT, diversity indices (species richness, Shannon-Wiener and Evenness indices) were computed. The significance of the effects of land use change on earthworm ecological metrics (density and biomass, species richness, diversity indexes) and soil organic carbon (SOC) was tested using a separate Generalized Linear Mixed effects models (GLMM). LUT effects were considered as fixed, and soil chemical parameters (total N, P, K, and pH) as random factors to account for unknown heterogeneity effects. SOC was analyzed as continuous variable, and thereafter applied GLMM to Gaussian distribution after log-transformation. Population density and biomass, and species richness were modeled as count data. The over-dispersion in earthworm density and biomass was tested using the qcc package in the R software. In cases of over-dispersion, the fits of Poisson regression and quasi-Poisson regression were compared with negative binomial regression. Parameters of the mixed-effects models were estimated using lme4 package with the restricted maximum likelihood (REML) method (Bates et al., 2015), and p-values computed based on the Satterthwaite approximations to the degrees of freedom in the lmerTest package (Kuznetsova et al., 2016). The conditional (variance explained by fixed and random factors) and marginal (variance explained by fixed effects only) R 2 -values were calculated following Nakagawa and Schielzeth (2013).
The possible sources of variation in soil chemical and physical characteristics were further investigated focusing on species richness, earthworm density, earthworm biomass and land use type. Using boxplots, we first examined the effect of LUTs on earthworm communities and on soil chemical and physical characteristics. The R package ggpubr was used to compare means. Correlations between earthworm density and soil variables were determined by Pearson's correlation coefficient.
A Biplot Principal Component Analyses were performed using the FactoMineR package to determine the relationships that could exist between earthworm species along with functional groups and soil parameters.

Soil Health Deterioration Index
Changes in soil health along the chronosequence from forest to aged plantations were assessed using degradation/deterioration indices (DIs). Soil degradation index for each soil property was calculated as the difference between mean values of individual soil properties under each LUT and the reference values of corresponding soil parameters in the baseline secondary forest (Lemenih et al., 2005;Dawoe et al., 2014;Tondoh et al., 2015). This index is expressed as a percentage of the mean values under the natural ecosystem. Furthermore, a cumulative DI was obtained by summing up the resultant positive and negative DI's of the individual soil properties for each farm field to be used as an index of soil quality responses (either degradation or improvement) to forest clearing and subsequent cultivation. Soil pH values were not considered in this calculation because the criteria of "more is better" is not true or uncertain over the range of values found in this study (Islam and Weil, 2000).

Interactive Feedback Between Earthworms and Soil Properties
Structural equation modeling was used to investigate the relationships between the abundance of earthworm functional groups (See Dataset S1) and soil properties including pH, soil organic carbon, total nitrogen, total phosphorus, mean weight diameter and bulk density (See Dataset S2, S3). Two models were considered: the first evaluating the direct and indirect effect of soil properties (here considered as reflective indicators) on earthworm density through correlation path analysis; and the second looking at the feedback effect of earthworms on soil health (expressed as latent variable manifested by soil properties). Data were standardized to homogenize measurement scales and to keep the same magnitude order between the observed variances. SEM was carried out with the package "Lavaan" (Rosseel, 2012), and path diagrams were generated with the package "semPlot" (Epskamp, 2015). Evaluation of the model fit was based on the analysis of Chi square, Tucker-Lewis index (TLI), comparative fit index (CFI), and mean square error of approximation (RMSEA).
All statistical analyses were carried out using the R software (R Core Team, 2018).

Community Structure of Earthworms in the Landscape
Up to 24 earthworm species belonging to four families, namely Glossoscolocidae, Acanthodrilidae, Eudrilidae, and Ocnerodrilidae, were collected in the entire landscape ( Table 2). The Acanthodrilidae family harbored the most important community composed of 15 native species while the Eudrilidae's family is composed of 7 species in the entire landscape. It is noteworthy that the whole community is composed of the pantropical species Pontoscolex corethrurus of the Glossoscolocidae family and Eudrilus eugeniae and Hyperiodrilus africanus of the Eudrilidae family whose distribution is spread in degraded lands in Africa. From a functional viewpoint, earthworms are classified into detritivores, geophageous oligo, meso, and polyhumics ( Table 2). The biplot CA performed on the 24 earthworm species (Figure 2) evidenced significant association between species and LUT (χ 2 = 946.54, p < 0.001). As a matter of fact, along the first axis (47.7%, eigenvalue = 0.51), species associated with food crops and the plantations of 7 years, including H. africanus, E. eugeniae, M. omodeoi, Dichogaster terraenigrae, Millsonia sp, Dichogaster mamillata, Gordiodrilus paski, are opposed to those attached to forest lands and the two aged rubber tree stands. On the contrary, axis 2 (33.15%) revealed the opposition of forested areas characterized by Stuhlmannia palustri, D. erhrhardti, D. papillosa, D. baeri, S. compositus, and rubber tree plantations of 12 years where S. zielae, D. eburnea, P. corethrurus, and D. leroyi are common.
With regards to functional attributes, in line with density of the polyhumics, the geophageous oligohumics did not show significant changes in density (Figure 3) and biomass (Figure 4) after the conversion of forest into rubber tree plantations. Conversely, significant changes were found in the density and biomass of detritivores (p < 0.0001) and mesohumics (p = 0.024); and the biomass of polyhumics (p = 0.0001). As for density, detritivore earthworms represented the most important group in several land uses (62.7-63.2%) except the youngest plantations (30.1%) and food crops (2.0%), where they were less present (Figure 3). On the contrary, mesohumics were fairly wellrepresented in cultivated areas including food crops (40.2%), 7year (46. %) 25-year-old plantations (28.3%) while polyhumics were strongly associated with food crops (57.8%) and the youngest plantations (23.0%). Similar trends were found with the biomass (Figure 4) as detritivores accounted for 50.7-94.2% in forest and the two last plantations in the chronosequence, mesohumics being in the range 47.3-66.6% (food crops, rubber 7 and 25 years) and the polyhumics being mostly harbored by food crops at 36.2%.

Changes in Earthworm Communities
Generally speaking, unlike the Evenness index, ecological metrics values revealed significant changes in earthworm communities over 25-years along the chronosequence in rubber-dominated landscapes ( Figure 5). The highest average values of density were

Changes in Soil Chemical Characteristics
The conversion of forest into rubber-dominated landscapes resulted in significant shifts (p < 0.001) in the average values of soil chemical characteristics as portrayed in Figure 6. The greatest value of SOC (22.5 g kg −1 ) was recorded in forest soils, followed by intermediate values in 7-year (12.8 g kg −1 ), 12-yearold plantations (14.5 g kg −1 ), and food crops (12.3 g kg −1 ) and the lowest values were recorded in the 25-year old plantations (9.9). Total N showed a similar trend with the highest value beneath forest (1.93 g kg −1 ) and the lowest in the 25-year old plantations. The value of P significantly increased along the chronosequence with greater ones found in the 25-year old plantations (515.9 mg g −1 ), food crops (391.6 mg g −1 ), and the 12-year old plantations (299.6 mg g −1 ). On the contrary, K did not show a decreasing trend, recording the highest values in the forest and food crops and lowest values over ages of plantations. After an initial low average value in forest (4.43), pH values increased up to the top level in food crops (6.36) before undergoing a steady drop in the plantations (5.71, 5.08, 4.84).

Changes in Soil Physical Characteristics
The distribution of aggregates size significantly (p < 0.0001) varied across LUTs in the 0-10 cm and 10-20 cm soil layers. The impact of forest conversion into rubber tree landscapes was mostly reflected in the macroaggregates (>2,000 µm) that showed a sharp drop in food crops followed by a recovering along with a steady increase in the plantations. The aggregate size classes 1,000-2,000 µm and 500-1,000 µm were the most important in the landscape, although they did not show a clear trend revealing the impact of forest conversion. Overall, the mean weight diameter (MWD) in top and subsoil varied significantly across LUTs (p < 0.001) as shown in Figure 9.
The lowest values of MWD in top (2.47) and subsoil (2.62) were recorded in food crops, while the highest values were recorded in forest (6.59 and 7.98) and last plantations in the chronosequence. Intermediate values were found in the last two plantations of the chronosequence (Figure 7). As for bulk density, values in the 0-10 cm layer was low (1.14) in the forest before a significant rise occurred in the food crops (1.37) with inconsistent trends in plantations (1.36, 1.22, 1.26) thereafter. On the contrary, no significant changes occurred in the 10-20 cm layer.

Soil Health Degradation Indices
The analysis of soil deterioration indices depicted in Table 3 revealed that soil health was severely impaired in food crops (DI = −173.7) and the 7-year old plantations (DI = −345.8). However, this situation reversed after 12 years (DI = −69.9) before full restoration in plantations of 25 years (+84.1). Earthworm density and species richness along with total P were instrumental in reversing the trend of soil health degradation caused by significant drops in SOC, total N and extractable K along the chronosequence ( Table 4).

Relationship Between Earthworm Abundance and Biomass and Soil Organic Carbon
The results of GLMM showed that SOC, earthworm density and biomass varied significantly among LUTs, with 42, 13, and 19% of variance explained, respectively ( Table 5).
For SOC, forest had a regression coefficient which was 14.06 significantly higher than that of the other LUT ( Table 3), indicating that SOC was higher in forest followed by rubber plantations in a decreasing order and food crops. A similar trend was observed for earthworm biomass while density displayed the highest coefficient models in 25 and 12-year-old plantations, indicating higher values in these two plantations as opposed to lower values in 7-year-old stands and in food crops.

Relationship Between Earthworm Species and Soil Parameters
The first two axes of the biplot PCA used to test the interaction between earthworm species and soil parameters accounted for a total variance of 82.1% (Figure 8). Axis 1 (54.8%) represented LUTs rich in SOC, total N, P with moderate values of mean weight diameter which are significantly associated with Dichogaster  (Figure 8).

Interactive Feedback Between Earthworms and Soil Properties
The direct influence of soil physico-chemical parameters on earthworm functional groups, evaluated from the first path analysis, revealed that SOC (T = 0.5, p = 0.0000), MWD_top soil (T = 0.40, p = 0.008), K (T = 0.32, p = 0.024), and P (T = 0.21, p = 0.04) have a strong causal effect on geophageous mesohumic earthworms, while K and MWD_sub show negative correlations with detritivores (T = −0.33, p = 0.039) and polyhumics (T = −0.34, p = 0.04), see Figure 9 and Table 6. However, there is no significant relationship between the abundance of oligohumic earthworms and the physico-chemical soil parameters considered in this study. The second structural model, proposed for assessing feedback effects of earthworms on soil health, was characterized by a poor overall fit. The low values obtained for all evaluation indices including the Comparative Fit Index (CFI: 0.558), the Tucker-Lewis Index (TLI: −1.16), and the Root Mean Square Error (RMSEA: 0.355) indicate that the proposed model was not able

Food crops
Est. and SE represent estimates of regression coefficient and stand error, respectively.
Frontiers in Environmental Science | www.frontiersin.org to accurately reflects the variability observed in the data. Direct effects of earthworms on soil health were not clearly evidenced by the results (Figure 10).

Changes in Earthworm Communities Along the Chronosequence
Apart from the presence of Pontoscolex corethrurus and Eudrilus eugeniae, species composition of earthworm communities was similar to those collected in southern west Côte d'Ivoire (Tondoh et al., 2011(Tondoh et al., , 2015. The conversion of forest into rubber tree landscapes has resulted in a significant shift in earthworm communities both at taxonomical and functional levels. After a significant drop in the food crops and the 7year-old plantations, population density of earthworms showed a significant increase in the 12-year and the 25-year-old plantations. These increases did not differ significantly from that of the forest. Similar trends were reported in previous studies in western Côte d'Ivoire in cocoa growing landscapes (Tondoh et al., 2011(Tondoh et al., , 2015. This trend is likely due to the proliferation and the mixture of the pantropical earthworm P. corethrurus, the African-Wide earthworms Hyperiodrilus africanus and native species (Dichogaster saliens and Dichopaster eburnean) in the two last plantations of the sequence. Growth and expansion of these populations are most likely due to their capability of withstanding degraded agro-ecosystems with low soil organic carbon content (Marichal et al., 2010;. This explains the higher species richness and diversity in plantations due to gradual enrichment or increase in species composition consistently with findings in cocoa landscape in southwestern Côte d'Ivoire (Tondoh et al., 2015). Another explanation is supported by the high soil water content in the 25-year-old plantation characterized by great production of litter (Chaudhuri et al., 2013;N'Dri et al., 2018) of good quality due to its low content in polyphenol, flavonoid and lignin (Chaudhuri et al., 2013). Moreover, the same trend of increased density of earthworm populations was found along rubber chronosequence in southern Côte d'Ivoire (Gilot et al., 1995;N'Dri et al., 2018) and in Tripura, India (Chaudhuri et al., 2013). The trend in biomass was different with higher values in the forest compared to human-derived ecosystems, mostly characterized by the prevalence of small and medium-sized species. It's noteworthy that in the current study, P. corethrurus did not represent the dominant species of the earthworm community as it is in FIGURE 9 | Correlation Path diagram relating soil parameters as exogenous variables and earthworm functional groups as endogenous variables. The value indicates the path coefficients (T). Node opacity indicates their significance (p-value); the width its importance and the color its direction (green for a positive correlation, red for a negative). detri_d, meso_d, oligo_d, poly_d represents respectively density of detritivores, geophageous mesohumic, oligohumic, and polyhumic earthworms. SOC, soil organic carbon; bulk density, soil bulk density, MWDtop and MWDSub, mean weight diameter of top and sub soil; TotalN , total soil nitrogen; K, exchangeable potassium; and P, total phosphorus. disturbed agro-ecosystems. It has been reported in rubber tree plantations in India and across the Amazonia that this species populations represented up to 70% total earthworm populations in density and biomass, respectively (Chaudhuri et al., 2008(Chaudhuri et al., , 2013Marichal et al., 2010). The detritivores and the geophageous mesohumic were the most important functional groups in forestderived landscapes with a marked presence of native species that were composed of earthworms from both functional groups.

Impact of Land Use Change on Soil Properties
Soil organic carbon (SOC), pH, and total phosphorus (P) showed significant variations in their values along the chronosequence, indicating consistent changes after forest conversion into rubber tree plantations. These findings agree with previous research in southwestern Côte d'Ivoire in cocoa landscapes (Tondoh et al., 2011(Tondoh et al., , 2015, the humid forest zone of Nigeria (Oku et al., 2012), in Asia (de Blécourt et al., 2013) and in the Western Kenya highlands (Nyberg et al., 2012). SOC concentration was high and above 20 g kg −1 , which is considered as the threshold value of a soil of good quality (Musinguzi et al., 2013;Lal, 2015), indicating deterioration of SOC by rubber farming. Conversely, the steady increase of P in rubber tree plantations over time was noteworthy and is likely to be factored into fertilizer management. Indeed, with money made out of their plantations, farmers were keen on maintaining their soil fertility status by integrating mineral and organic fertilizers in alleys as shown in legume-based cropping systems in middle Côte d'Ivoire (Koné et al., 2008). Contrasting soil pH values across the landscape revealed that soil acidity should be handled with care as it tended to be low in forest and rubber tree-derived ecosystems. Findings pertaining to SOC, pH and bulk density confirmed results by N'Dri et al. (2018), who found similar trends along another chronosequence in the same study site. The conversion of forest into rubber tree landscapes has resulted in a significant drop in aggregate stability and large macroaggregates (>2 mm) in the food crops prior to a gradual replenishment over time in both layers (0-10 cm and 10-20 cm) of the plantations. In contrast, bulk density was low in the forest but was high under human-derived systems most likely indicating impaired soil porosity in the plantations (Kakaire et al., 2015). Furthermore, the proportion of macroaggregates was severely reduced in the food crop systems and afterwards restored in the plantations. These findings point out the beneficial role of trees through their roots in shaping up soil structure. The following explanations likely account for the following processes: (i) the huge abundance of soil organisms (macro and meso invertebrates, particularly soil mites) that increased at a rate of +121 % in the 25-year old plantations and are known to be active in the breakdown of fresh inputs of organic material abundantly produced in 12 and 25-year-old plantations (N'Dri et al. (2018); (ii) soil aggregation in aged plantations is facilitated by vegetation restoration that caused huge fresh organic matter returns amounting to 3.9 ± 0.1 -5.1 ± 0.6 t ha −1 year −1 (N' Dri et al., 2018) to enhance the aggregation of soil particles with plantation stand age (Bronick and Lal, 2005). In contrast, soil exposure along with lack of residue inputs in the food crops caused declines in aggregation and organic carbon, both of which make soil susceptible to erosion (Pinheiro et al., 2004). Furthermore, the increased fresh organic C in 12 and 25-year old plantation might have enhanced the stability of the aggregates through the binding of mineral particles and the formation of stable aggregates as reported by Demenois et al. (2018) in New Caledonia.

Drivers of Soil Health Change in Rubber Tree Landscapes
The conversion of forest into rubber tree plantations significantly impaired all soil biological, physical and chemical parameters in the first place. As a result, soil health was heavily degraded in food crops and 7-year old plantations. Similar results were reported in cocoa landscapes after 7 years of cocoa cropping in center-west Côte d'Ivoire (Tondoh et al., 2015) and after 3 years in the Ashanti Region of Ghana (Dawoe et al., 2014). The improvement of soil health from 12 years in the rubber tree landscapes is consecutive to the increase in earthworm abundance and P concentration. To some extent, the gentle increase in SOC, the significant rise in large soil aggregate proportion and aggregate stability under the aged plantations contributed to improve soil health status. These findings revealed the beneficial impact of rubber stands on soil health improvement in the long run as the lifespan of the cropping cycle is estimated at 40 years. Indeed, Dawoe et al. (2014) reported an improved soil quality after 30 years of shaded-cocoa farming in the Ashanti region of Ghana owing to the re-accumulation of 85% of the initial forest soil organic carbon stock. Similarly, it was revealed that up to five decades was necessary to significantly improve SOC beneath tree plantations derived from croplands (Sauer et al., 2012). It is therefore not surprising that alternative cropping options, such as rubber-based agroforestry systems (Hevea brasiliensis-Theobroma cacao and H. brasiliensis-Dalbergia cochinchinensis), have the ability to improve soil quality and ecosystem resilience (Chen et al., 2017).

Can Earthworms Play a Troubleshooter's Role in Mitigating Soil Health Deterioration Issues in Rubber Tree Landscapes?
The Biplot-PCA revealed that in the rubber tree landscapes, earthworms' assemblage was featured by land use change, which in turn strongly impacted the soil chemical and physical characteristics. As a result, two functional groups of earthworms were found to be involved. The first group was made up of Dichogaster baeri, Millsonia omodeoi and Stuhlmannia zielae, which are associated either with LUTs rich in SOC, total N and P (i.e., forest and 25-year old rubber tree plantations). The second group consisted of Pontoscolex corethrurus, Hyperiodrilus africanus and Dichogaster papillosa, commonly found in degraded LUTs-food crops and 7-year old plantation-with high pH and bulk density along with moderate total P. Unlike M. omodeoi that was found in forested areas, a similar association was reported in the center west cocoa landscapes of Côte d'Ivoire . Indeed, as highlighted by Koné et al. (2012), M. omodeoi can be viewed as a persistent species due to its ability to adapt to contrasted environments. Furthermore, the strong relationship between earthworms and SOC agrees with previous findings (Tondoh et al., 2011; and confirms the role of SOC as a key driver of earthworms' abundance and community structure in agro-ecological landscapes. Although these findings were expected, they don't provide much clarity on the direction of the interaction between both abiotic and biological components. The structural equation modeling was used to shed light on the relationship between earthworms and soil physico-chemical parameters by testing two hypotheses: the first highlighting the strong influence of soil physico-chemical parameters on earthworm communities, and the second addressing the potential FIGURE 10 | Path diagram relating earthworm functional groups as exogenous variables and soil health as endogenous variables. The value indicates the path coefficient (T). Node opacity indicates their significance (p-value); the width its importance and the color its direction (green for a positive correlation, red for a negative). See Figure 9 for the explanation of abbreviations. feedback effects on these same soil properties, caused by earthworm activities. The results of the correlation path analysis derived from the first hypothesis evidenced the welldocumented direct dependency between soil parameters, notably soil organic carbon on abundance of earthworm functional groups. Unfortunately, the poor overall fit of the measurement model, proposed to test feedback effect of earthworms on soil health, as assessed by the low values obtained for all evaluation indices (0.558, −1.16, and 0.355, respectively for CFI, TLI, and RMSEA), indicated that the suggested model did not properly depict the network of relationships between earthworms and soil health. However, it is now widely recognized that as ecosystem engineers, earthworms modulate the availability of soil nutrients, including organic carbon (Lavelle et al., , 2016. This pattern appears to be reflected in the model results, which indicate a greater influence of mesohumic earthworms on soil organic carbon and total phosphorus than other functional groups. But, the poor fit of the suggested measurement model indicates that the right paths of these interactions are still not wellidentified. The findings suggest a more complex network of relationships than the one proposed in this work. Indeed, one of the weaknesses of this study lies in its observational nature that could not provide insights into causal inference that is mostly driven by land use change to the detriment of soil organisms. Future research should focus on adjusting the model structure by identifying relevant interaction paths that explain the effect of these macro-invertebrates on the soil.

CONCLUSION
The study confirmed the detrimental impact of forest conversion into rubber tree landscapes in southern Côte d'Ivoire on soil health within the first 7 years and the restorative trend taking place beneath the plantations in year 12 and above. This was due to an enabling microenvironment (improved SOC, aggregate stability, exchangeable K, total phosphorus) taking place in rubber tree plantations, which was key for the development of geophageous mesohumic earthworms. In turn, they might have improved soil health through their interaction with soil organic carbon and total phosphorus. However, there were no evidence of direct effect of earthworms on soil heath, suggesting that more investigations at mesoscale is worth being undertaken.

AUTHOR CONTRIBUTIONS
JT conceived the study, designed, and coordinated the writeup of this manuscript. KD and YB contributed to statistical analyses and the review of the manuscript. AG, LA, and JN were involved in data collection and preprocessing and the review of the manuscript. GF was tasked with the proofreading and review of the manuscript.

FUNDING
The field survey and laboratory analyses of this study were supported by the French Ministry of Foreign and European Affairs, through the research program AIRES-Sud, implemented by the Institut de Recherche pour le Développement (IRD-DSF).