How 75 years of rubber monocropping affects soil fauna and nematodes as the bioindicators for soil biodiversity quality index

ABSTRACT

Natural rubber consumption has led to the expansion of rubber (Hevea brasiliensis) plantations which affects the deforestation and ecosystem. However, no study of the long-term effect of rubber plantations on soil biodiversity has been carried out yet. This study aimed to assess the long-term impact of continuous rubber monocropping on soil biodiversity, focusing on soil macrofauna and nematode diversity. Three successive rubber rotations at young and old ages were compared with the adjacent forest in Suratthani province, Thailand. Soil biodiversity quality index was calculated from a set of indicators which were combined into a single score to present a functional assessment of the gradient of disturbance. The results showed three negative effects on soil biodiversity (i) the biodiversity quality index immediately declined after deforestation (ii) the old age rubber plantations had a lower soil biodiversity as the nematodes were a main driver of diversity in the young plantation, and (iii) similarly, for the soil chemical properties, the long-term effect of rubber chronosequence evidenced deterioration in the third rotation. Therefore, two rotations of rubber plantation (around 50 years) seemed to be the maximum length of rubber monocropping in terms of soil biodiversity recovery.

KEYWORDS:

• Rubber (Hevea brasiliensis)
• monocropping
• soil biodiversity
• macrofauna
• nematode

Introduction

The rubber tree (Hevea brasiliensis), which is the unique commercial source of natural rubber, is one of the most important industrial crops in Asia (FAO 2012). The industrial-scale and smallholder monoculture of rubber plantations were developed during the strong demand for natural rubber (Warren-Thomas et al. 2015). This widespread agricultural system has expanded faster than all other tree crops, with a 1.8 fold increase in the last 30 years (Fox and Castella 2013). Several ecological studies have documented a negative impact of deforestation for rubber monoculture on soil quality in terms of soil pH, soil carbon stock, erosion and hydrology (Li et al. 2008; de Blécourt et al. 2013; Ahrends et al. 2015; Drescher et al. 2016; Liu et al. 2019). In addition, the conversion of land use to monoculture plantation leads to the degradation of ecosystem functions and the loss of tropical biodiversity (Fitzherbert et al. 2008).

Nowadays, biodiversity loss has become a major issue in terms of sustainable development, and halting it is one of the United Nations’ sustainable development goals (SDGs) (Bach et al. 2020). Soil biodiversity is important for agricultural ecosystems as the soil is home to more than a quarter of the earth’s total biodiversity (Turbé et al. 2010). Soil biota has three main contributions to the soil. First, soil microorganisms (i.e. bacteria and fungi) and microfauna (i.e. protozoa and nematodes) can transform organic and inorganic compounds which are essential for nutrient cycling (FAO 2020). Second, all soil organisms control the energy cycles by their food webs (Griffiths et al. 2000). And finally, soil macrofauna such as earthworms act as ecosystem engineers to modify soil porosity, water and gas transport (Blanchart et al. 2004; Schneider and Schröder 2012). Therefore, soil biodiversity plays a crucial role in various processes in the soil and has a direct influence on soil fertility (Villenave et al. 2004).

Thailand has been a major producer of the world’s natural rubber since the early twentieth century (IRSG 2015). During the initial of ‘rubber boom’, the rubber plantation in Thailand had started in forest areas. However, most current rubber areas have been in two or three replanting cycles (Chambon et al. 2016). As the typical economic lifespan of a rubber plantation is around 25 years, then two or three rotations of rubber replanting on the same land could lead to overexploitation of soil nutrients and significant loss of soil fertility and biodiversity. This study aimed to examine the long term effect of rubber monoculture on soil biodiversity, focusing on soil macrofauna and nematode diversity as the bioindicators. To address the impact of rubber replanting compared to the deforestation effect, three successive rotations of rubber monocropping (in both young and old age rubber plantations) were compared with the adjacent natural forest, focusing on the main patterns of soil macrofauna and nematodes. Soil invertebrates are the main biological component that regulates many soil processes (Brussaard 1998). They contribute to improving soil aggregate structure and are involved in soil fertility maintenance (Rousseau et al. 2013; Melman et al. 2019). Soil nematodes are the indicators of soil food web conditions which have subsequent effects on ecosystem variables and processes, such as plant and soil microbial biomass, nitrogen mineralisation and decomposition (Neher 2010; Jiang et al. 2017). A soil biodiversity quality index was calculated at the end to integrate the measured soil macrofauna and nematodes parameters into a single score that could be used as an indicator of overall soil biodiversity quality in this context.

Materials and methods

Experimental site and sampling

This study was conducted in smallholding rubber plantations (Hevea brasiliensis) in a major rubber growing area in Suratthani province, southern Thailand and an adjacent natural forest nearby the rubber plantations (8°54′–8°56′N; 99°24′–99°28′E). The initial rubber planting was established by clearing the forests in the area. The area is generally dominated by a tropical monsoon climate with two main seasons: the dry monsoon (December to April) and the wet monsoon (May to November). All samples were collected in September 2015 when the period of light rain precipitation. The mean monthly precipitation in this month was 250.9 mm and the average temperature was 26.8°C (Meteorological Department 2015).

The experiment was designed to assess both the impacts of disturbances following (i) the land clearance – from forest to the first rubber plantation and then to the next second and third rotation – and (ii) the evolution of soil fauna and nematode during the lifespan of each rotation. Therefore, two classes of rubber plantations were chosen according to age ranges: young stage (3–6 years old) and old stage (18–22 years old). The plantations were also categorised based on the number of replanting cycles with three successive rubber plantations: the first rotation – planted after clearing the forest area, the second rotation – planted after the first rubber rotation and the third rotation – planted after the second rubber rotation. Seven treatments with three plots (true replications) of each were assessed. The seven treatments were: Forest (F), the first rotation at the young stage (R1y), the first rotation at the old stage (R1o), the second rotation at the young stage (R2y), the second rotation at the old stage (R2o), the third rotation at the young stage (R3y), and the third rotation at the old stage (R3o) (Figure 1). In the plantations, all samples were collected at 0–10 cm of soil depth because the topsoil layer contains the main part of the soil biota that drive soil processes (European Commission 2010). The samples were taken from the inter-row of rubber trees, as agricultural practices were typically carried out in the inter-row spaces (Kazakou et al. 2016). However, in the forest, all samples were collected between the trees at the same depth.

Figure 1. Experimental design of the study. F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

Soil analysis

In each plot, soil samples were composited from three subsites consisting of 10 rows of rubber trees. Then, composite soil samples were sieved using a 2 mm mesh and dried at room temperature before laboratory analysis performed by the Office of Science for Land Development, Land Development Department, Bangkok, using the following methods: total nitrogen (N)(%) was determined using the Kjeldahl method (AOAC 1984); available phosphorus (P) (mg kg⁻¹) was determined using the Bray II method (Bray and Kurtz 1945); available potassium (K) (mg kg⁻¹), calcium (Ca) (mg kg⁻¹) and magnesium (Mg) (mg kg⁻¹) were extracted by neutral 1 N ammonium acetate (Chapman 1965) and analysed by flame photometer for potassium and atomic absorption spectrophotometer for calcium and magnesium; organic matter (OM) (%) was determined by using the Walkley and Black method (Walkley and Black 1934); and pH was determined in distilled water (1:1 soil-water ratio) (Peech 1965).

Soil fauna assessment

Soil macrofauna

In each plot, soil macrofauna was sampled using the soil monolith collection from a topsoil area of 25 × 25 cm with 10 cm depth (Anderson and Ingram 1993). The samples were individually hand-sorted in the fields and all visible organisms were collected and stored in 70% alcohol solution. All invertebrates were counted and identified in a laboratory according to trophic groups: (a) predator, (b) detritivore, (c) phytophagous, (d) omnivore, and (e) geophagous (Capinera 2008).

Nematodes

Nematodes were extracted from 100 g of fresh soil samples using the elutriation method. With this method, nematodes were separated from soil particles that settled faster and floated in the water. Then, all extracted nematodes were counted and identified under a dissecting microscope (Seinhorst 1956). Nematodes were identified by feeding group which can be deduced from the structure of their mouthparts: (a) bacterivore, (b) fungivore, (c) plant parasitic, (d) carnivore, and (e) omnivore (Bongers and Bongers 1998).

Biodiversity indices

The biodiversity indices were calculated from the soil macrofauna and nematode datasets:

1. Abundance showed the total number of species in the sample.

2. Richness index – which corresponded to the number of taxons detected in the sample (Ludwig and Reynold 1988).

3. Diversity index or Shannon index ($H^{\mathrm{\prime }}$) – which synthesised taxon richness and their relative abundance (Van Strien et al. 2012). $H^{\mathrm{\prime }}{H}^{\mathrm{\prime }}=\sum _{i=1}^s{p}_i\ln (p_i)$where $p_i$ is the proportion of each species in the sample and $\ln (p_i)$ is the natural logarithm of this proportion.

4. Evenness index or Pielou index ($J^{\mathrm{\prime }}$) associated with the Diversity index, assessing the structural diversity from 0, where few species dominate the community, to 1, where all the species are equally represented (Pielou 1966). $J^{\mathrm{\prime }}={\displaystyle \frac{H^{\mathrm{\prime }}}{H_{max}^{\mathrm{\prime }}}}$where $H^{\mathrm{\prime }}$ is the number derived from the Diversity index and $H_{max}^{\mathrm{\prime }}$ is the maximum possible value of $H^{\mathrm{\prime }}$.

Quality index calculation

A quality index was calculated in three steps following Obriot et al. (2016) to combine all biodiversity data sets into one score. The first step was to select the minimum data set. Pearson correlation matrix was calculated and only parameters with correlation coefficients (r) less than 0.8 were kept. The second step was to calculate the normalised indicator score ($Si$). The linear response curve was used to give the indicator. In this case, it was assumed that higher values would result in a better impact on soil quality for all parameters. To normalise data, the observed values ($a$) were divided by the highest value ($a_{max}$) $Si={\displaystyle \frac{a}{a_{max}}}$

In the third step, all selected parameters were weighted by principal component analysis (PCA) to consider the weighted factor ($Wi$) by multiplying the sum of square coordinates of an indicator on each eigenvector ($j$) (only the components with eigenvalue more than 1 were kept) by the percentages of the total variability of each principal component ($fj$). $Wi=\sum _{j=1}^p\lambda j\times fj$where $Wi$ is the weighted factor, $j$ is the sum of squared coordinates of an indicator on each eigenvector, $fj$ is the percentages of the total variability of each principal component and $p$ is the number of principal component.

Then the quality index (QI) was calculated by the weighted factor and the normalised indicator score. $QI=\sum _{i=1}^{n}Si\times Wi$where QI is the quality index, Si is the normalised indicator score, Wi is the weighted factors, and n is the number of indicators in the data set.

Statistical analysis

All statistical analyses were performed using R v3.1.0 (R Development Core Team 2014). Normality and homogeneity of variance were tested in Shapiro and Bartlett tests on model residuals. When one or both conditions were not met, the data were log-transformed. A Linear Mixed Model was generated for quantitative data and a Generalised Linear Mixed Model for count data. One-Way ANOVA was used and followed by Tukey multiple comparisons of means tests to compare differences between all the treatments with a Bonferroni correction. Between-Class Analyses (BCA) were performed on the data based on Principal Component Analyses (PCA) results with consideration of variables. Monte-Carlo permutation tests (999 permutations) were generated to assess the significance of the multivariate representation. Pearson’s correlation coefficients were used to assess the relationships of the data sets.

Results

Soil analysis

All soil properties in every treatment were significantly different (Table 1). The highest values of all soil properties were found in the forest, except the available phosphorus. The results showed that the soil fertility decreased immediately after deforestation and land use change. The decreasing of cation contents (K⁺, Mg²⁺, and Ca²⁺) were similar after the deforestation and each rotation. Old plantations seemed to have a higher value of cation contents than those of young plantations. Furthermore, soil degradation was found in the third rotation.

Table 1. Descriptive statistics of the measured soil properties in each plantation; total nitrogen (N), available phosphorus (P), available potassium (K), calcium (Ca), magnesium (Mg) and organic matter (OM).

Plantation	N	P	K	Ca	Mg	pH	OM
	(%)	(mg kg^−1)	(mg kg^−1)	(mg kg^−1)	(mg kg^−1)		(%)
F	1.7 ± 0.2a	7.0 ± 2.5bc	92.5 ± 46.2a	386.0 ± 225.5a	93.3 ± 10.5a	5.7 ± 0.6a	2.2 ± 0.2a
R1y	1.3 ± 0.4ab	21.1 ± 1.6a	78.7 ± 31.0a	171.8 ± 77.3abc	43.7 ± 28.0b	5.2 ± 0.3ab	1.6 ± 0.4ab
R1o	1.1 ± 0.1abc	4.4 ± 1.1bc	61.0 ± 14.3ab	164.3 ± 97.7abc	43.5 ± 16.0b	5.3 ± 0.1ab	1.7 ± 0.4ab
R2y	1.0 ± 0.6abc	12.5 ± 7.5ab	81.2 ± 34.4a	188.0 ± 106.2ab	39.5 ± 13.3b	5.5 ± 0.5ab	1.7 ± 0.1ab
R2o	1.0 ± 0.4abc	5.5 ± 5.0bc	78.5 ± 97.0a	67.3 ± 16.4bc	25.8 ± 6.3bc	5.1 ± 0.4ab	1.8 ± 0.5ab
R3y	0.4 ± 0.1c	4.3 ± 2.3bc	16.5 ± 9.2b	61.7 ± 32.3bc	15.2 ± 3.5cd	5.0 ± 0.1ab	1.2 ± 0.3b
R3o	0.8 ± 0.1bc	3.2 ± 1.0c	13.8 ± 2.1b	28.7 ± 11.4c	11.0 ± 1.5d	4.7 ± 0.3b	1.4 ± 0.3b

Notes: Each value represents the mean (±SD) from composite soil samples (n = 3). The values followed by different letters indicate statistically different values (P < 0.05) among plantations within columns. F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

Soil fauna diversity

Soil macrofauna

About 30 taxa groups had been found and identified into five groups depending on their trophic groups: predator, detritivore, phytophagous, omnivore, and geophagous. Detritivore was the most abundant group in this experiment. Isoptera was the main fauna in this ecosystem. However, the significant differences were found only in the predator and geophagous groups. The number of predators seemed to be more in young plantations (R1y and R2y) as the habitat of weeds provides more aboveground insects which are their food. With dense shades and higher soil moisture content as the results of close canopies, the forest showed the highest number of geophagous or earthworm, followed by the old plantations (Table 2). However, it was found that the total abundance of soil macrofauna was not significantly different. Only the richness index, linked to the number of all taxa groups, was significant as the third rotation indicated the decreasing of richness index along the chronosequence (Table 3).

Table 2. The abundances of soil macrofauna in each trophic group.

Plantation	Predator	Detritivore	Phytophagous	Omnivorous	Geophagous
	(individuals m^−1)	(individuals m^−1)	(individuals m^−1)	(individuals m^−1)	(individuals m^−1)
F	204.4 ± 26.3ab	419.6 ± 313.4	37.3 ± 14.1	455.1 ± 388.6	222.2 ± 48.4a
R1y	238.2 ± 90.9ab	245.3 ± 261.3	37.3 ± 26.7	186.7 ± 203.3	110.2 ± 44.7ab
R1o	156.4 ± 57.0b	504.9 ± 639.1	5.3 ± 5.3	741.3 ± 1150.5	220.4 ± 29.4a
R2y	508.4 ± 388.1a	229.3 ± 202.6	30.2 ± 8.1	462.2 ± 32.6	40.9 ± 11.1b
R2o	186.7 ± 61.5ab	1050.7 ± 1410.6	21.3 ± 16.0	455.1 ± 542.6	144.0 ± 123.0ab
R3y	103.1 ± 34.7b	581.3 ± 692.3	10.7 ± 0.7	129.8 ± 67.1	56.9 ± 27.4b
R3o	167.1 ± 38.6ab	360.9 ± 306.3	7.1 ± 3.1	138.7 ± 113.5	90.7 ± 54.1ab

Notes: Each value represents the mean (±SD) from soil samples (n = 3). The values followed by different letters indicate statistically different values (P < 0.05) among plantations within columns. F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

Table 3. The total abundances of soil macrofauna and their biodiversity indices.

Plantation	Total abundances	Richness index	Diversity index	Evenness index
	(individuals m^−1)
F	1342.2 ± 721.1	11.7 ± 0.9a	1.7 ± 0.3	0.7 ± 0.1
R1y	837.3 ± 290.4	9.8 ± 2.0ab	1.6 ± 0.2	0.8 ± 0.1
R1o	1664.0 ± 1053.6	7.8 ± 0.8bc	1.4 ± 0.1	0.7 ± 0.1
R2y	1278.2 ± 390.3	9.2 ± 1.3abc	1.5 ± 0.1	0.7 ± 0.0
R2o	1950.2 ± 1248.0	10.4 ± 1.1ab	1.5 ± 0.3	0.6 ± 0.2
R3y	901.3 ± 675.8	6.6 ± 1.3c	1.2 ± 0.5	0.7 ± 0.2
R3o	768.0 ± 477.4	7.7 ± 0.7bc	1.6 ± 0.1	0.8 ± 0.0

Notes: Each value represents the mean (±SD) from soil samples (n = 3). The values followed by different letters indicate statistically different values (P < 0.05) among plantations within columns. F: Forest, R1y: 1st rotation at the young stage, R1o: 1st rotation at the old stage, R2y: 2nd rotation at the young stage, R2o: 2nd rotation at the old stage, R3y: 3rd rotation at the young stage, and R3o: 3rd rotation at the old stage.

Nematodes

Bacterivore was the most abundant trophic group of nematodes in the study. The number of bacterivore nematodes was the highest in the forest but decreased in the old plantation. This decreasing pattern was also found in omnivore, plant parasitic, the total abundance, and richness index of nematodes. The decreasing of richness gradually reduced along the rubber chronosequence after deforestation, and the lowest number of the richness index was found in the old stage of the third rotation (Tables 4 and 5).

Table 4. The abundances of nematode in each trophic group.

Plantation	Bacterivore	Carnivore	Fungivore	Omnivore	Plant parasitic
	(individuals 100g^−1 soil)	(individuals 100g^−1 soil)	(individuals 100g^−1 soil)	(individuals 100g^−1 soil)	(individuals 100g^−1 soil)
F	850.1 ± 492.3a	1.2 ± 2.0	60.1 ± 14.8	9.2 ± 2.5abc	197.0 ± 47.5a
R1y	436.9 ± 180.2ab	11.6 ± 13.4	66.3 ± 31.2	45.5 ± 44.6a	237.9 ± 77.3a
R1o	196.0 ± 72.6bc	4.7 ± 4.2	16.9 ± 5.1	2.3 ± 2.5c	125.4 ± 48.3ab
R2y	379.7 ± 65.5ab	7.1 ± 3.0	53.0 ± 37.5	16.9 ± 6.0ab	266.6 ± 189.2a
R2o	286.5 ± 58.7ab	0.9 ± 1.6	18.8 ± 12.3	6.3 ± 2.9abc	54.9 ± 29.1b
R3y	73.6 ± 30.5c	1.2 ± 1.2	45.6 ± 32.1	4.0 ± 3.2bc	86.0 ± 65.5ab
R3o	189.7 ± 142.4bc	2.3 ± 3.9	15.2 ± 8.5	3.3 ± 1.6bc	79.4 ± 32.6ab

Notes: Each value represents the mean (±SD) from soil samples (n = 3). The values followed by different letters indicate statistically different values (P < 0.05) among plantations within columns. F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

Table 5. The total abundances of nematodes and their biodiversity indices.

Plantation	Total abundances	Richness index	Diversity index	Evenness index
	(individuals 100g^−1 soil)
F	1117.5 ± 540.3a	18.7 ± 3.8a	1.8 ± 0.5	0.6 ± 0.1
R1y	798.2 ± 66.3ab	18.7 ± 2.1a	2.1 ± 0.4	0.7 ± 0.1
R1o	345.3 ± 48.3bc	19.3 ± 2.1a	1.9 ± 0.3	0.6 ± 0.1
R2y	723.3 ± 274.2ab	18.0 ± 1.0a	2.0 ± 0.2	0.7 ± 0.1
R2o	367.4 ± 89.1bc	15.0 ± 3.6ab	1.4 ± 0.4	0.5 ± 0.1
R3y	210.4 ± 44.9c	14.3 ± 0.6ab	1.8 ± 0.3	0.7 ± 0.1
R3o	289.8 ± 184.7c	12.0 ± 1.0b	1.6 ± 0.2	0.6 ± 0.0

Notes: Each value represents the mean (±SD) from soil samples (n = 3). The values followed by different letters indicate statistically different values (P < 0.05) among plantations within columns. F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

Biodiversity quality index

For biodiversity quality index calculation, seven variables from the soil macrofauna and nematode measurement were selected as the minimum of the data sets; the abundance of predator and geophagous feeding group and the richness index of soil macrofauna, the abundance of omnivore and plant parasitic trophic group in nematode with the total abundance and the richness index. The results indicated that the biodiversity quality index was significantly decreasing along the rubber plantation succession with the lowest value in the third rotation. And the richness index of nematode seemed to be a main driver of the biodiversity ecological system in this context, followed by the abundance of earthworm which was the geophagous feeding group (Figure 2(a)).

Figure 2. Biodiversity quality index along rubber chronosequence compared to forest. Part (a) was the biodiversity index from the quality index calculation (mean ± SD, n = 3). The different letters indicate statistically different values (p-value < 0.01). Part (b) was the individual representation and correlation circle of each biodiversity quality index variable from the Between Class Analysis (BCA). F: Forest, R1y: 1^st rotation at the young stage, R1o: 1^st rotation at the old stage, R2y: 2^nd rotation at the young stage, R2o: 2^nd rotation at the old stage, R3y: 3^rd rotation at the young stage, and R3o: 3^rd rotation at the old stage.

The percentages of inertia of two first dimensions of BCA showed the significantly different values (P-value < 0.05) with 54.12% of variability in the first axis and 24.84% in the second axis (Figure 2(b)). The BCA explained the correlation of each biodiversity quality index variables and the plantation into 4 clusters; (1) forest, (2) young age of the first and the second rotation, (3) old age of the first and the second rotation, and (4) young and old age of the third rotation.

Discussion

Relationships between soil biodiversity and soil properties

The results of this study showed a similar pattern of soil properties and soil biodiversity decreasing with time. A rapid decline occurred immediately after deforestation, while the deterioration was continued along the rubber chronosequence. This result can be explained by the fact that soil biodiversity is affected by the interaction of abiotic and biotic factors, on which land use change and degradation may have a negative effect (Tibbett et al. 2019). The explanation could be verified by Van Noordwijk et al. (1997) who found that the conversion of natural forest into an agricultural area leads to soil degradation due to the reduction of aboveground biomass and soil organic carbon. A similar result was also found after rubber monoculture perturbation, which is the cause of soil nutrient loss (Zhang et al. 2007). As soil hosts the largest diversity of organisms, it plays a key role for biogeochemical cycles which are related to regulating and supporting ecosystem services (Smith et al. 2015). In addition, soil biodiversity loss will affect soil processes such as the nitrification or litter decomposition, and also reduce some capacity to recover from the disturbances (Nielsen et al. 2015). Such observations were reported in Mulder et al. (2005) and Sinsabaugh et al. (2008) that soil pH seemed to strongly affect soil biodiversity as it relates to the soil fauna community by influencing the relative abundance and biomass of soil microbes and the activity of soil enzymes. However, this study revealed that the soil degradation varied by available phosphorus. The available phosphorus was less in the forest and decreased along with the plantation ages, because it was gradually depleted in soil solutions with time (Elser et al. 2007; Turner and Condron 2013).

Changing of soil macrofauna and nematodes in rubber monocropping plantation

Soil macrofauna

Soil macrofauna was selected to be a bioindicator for monitoring the effect of land use change in this study because of their sensitivity to the habitat change and the effect on the material cycle and energy flow in ecosystems (Tan et al. 2013; Wu and Wang 2019). The results showed the significant differences only in predator and geophagous trophic groups and the richness index. The increasing of predator in the young age of the first and the second rubber rotation was a result of the larger weed area that serves as food sources for more aboveground insects. Thus, the limitation in food sources is the main cause of biodiversity loss (Semper-Pascual et al. 2019). However, all predators decreased after the third rotation which can be explained by the decreasing of prey and the retardation in decomposition rate (Kajak 1995). Regarding the earthworm assemblages, Guéi and Tondoh (2012) reported that organic carbon was the most important variable for earthworm abundance, followed by nitrogen and pH. Thus, the higher value of these soil variables in the forest and the old plantation caused the higher abundance of the geophagous group. The forest and the old rubber plantation were significant in soil fertility richness through enhanced litter decomposition with more litter (Sofo et al. 2020; Xu et al. 2020). This reason could explain the increase of soil fauna richness index in the forest and the old plantation. However, it was observed that the soil degradation in the third rotation caused a decrease in soil fauna diversity.

Nematode

Nematode, one of the most numerous members of the soil fauna, is an indicators of the soil food web (Herren et al. 2020). It is sensitive to organic amendments and contributes to soil nutrient turnover linked to land use changes (Treonis et al. 2010). For example, the shift toward monoculture was a kind of soil disturbance that reduced the nematode growing (Duddigan et al. 2020). Like other previous studies of the land use changes effect on soil organisms, this study found significantly different values in the abundance of bacterivore, omnivore and plant parasitic feeding nematodes with the total abundance and the richness index. This study also indicated the decrease of nematode abundance in old rubber plantation because weeds grow less well under dense shades causing the loss of plant biomass, an important carbon source for the growth and activity of nematodes (Ferris 2010; Li et al. 2021). Also, the cover plants, directly and indirectly, affect the microclimate which influences the above and belowground organisms (Belnap et al. 2005).

Loss of biodiversity quality index along successive rubber rotations

The results showed that the deforestation, the age effect, and the rotation of rubber plantation were all related to the loss of biodiversity quality index. The results highlight the effect of the deforestation which caused the fragmentation and habitat loss of the biota (Ehrlich and Pringle 2008). The high abundance of earthworms, as the main driver of the biodiversity quality index in the forest, was due to more litter biomass with low soil temperature and higher moisture (Tien et al. 2000). However, Ganesh et al. (2009) reported that the litter with higher polyphenol and lignin observed in forests was not preferred by earthworms.

Since most of the other variables from all minimum of the data sets linked to nematode factors, the age effect showed the difference in the young stage of the first and the second rotation. The decreasing of the biodiversity quality index in the old stage of the plantation was caused by a decreased in plant cover, which affected the richness index of soil biota because of the decreasing ecosystem primary productivity (Ye et al. 2020). The study also found the opposite abundance of earthworm and nematodes during the ageing of the rubber plantation. Demetrio et al. (2019) explained this negative effect of earthworms on nematodes by suggesting that it related to food competition.

Concerning the rotation effect, the results showed that the biodiversity quality index was not significantly different in the first and second rotation, but considerably reduced in the third rotation. Fiorini et al. (2020) observed that soil fauna worked as a sensitive device to detect the soil quality change by soil tillage. This decrease may be due to the strong oxygenation from the tillage which helped accelerate the rapid mineralisation of organic matter, reduced the organic carbon content, and caused the nutrient leaching (Micucci and Taboada 2006; Angers and Eriksen-Hamel 2008). These affiliations were also reported in a study by Melero et al. (2009) that the increasing of microorganism biomass in soil under no-tillage cultivation enhanced soil organic carbon content. Thus, these earlier reports had explained the decreasing of biodiversity quality index after the land preparation by tillage in each rotation.

Moazzam et al. (2014) suggested that the longer cycle in each rotation (40 years) might allow the recovery of maximum carbon content in the rubber plantation as the organic carbon encouraged the soil nutrient cycling and linked to the food resources of soil organisms, Thus, the loss of soil biodiversity could be also restored. Meanwhile, agroforestry is an option to improve ecosystem services that play a critical role in sustainable agriculture (Haggar et al. 2019; Santos et al. 2019).

Conclusion

This study indicated the long-term impact of rubber monocropping on topsoil biodiversity particularly soil macrofauna and nematodes. In addition to the obvious degradation after deforestation, the study observed the pattern of soil biodiversity loss after three rotations of rubber plantations. The soil biodiversity quality index still maintained after the first rotation, and evidenced the deterioration after the second rotation, which is after about 50 years of rubber monocropping. Thus, to reduce this impact on the sustainable land use, alternative cultivations such as agroforestry or intercropping systems could be an attractive option to ensure the environmental benefits.

Acknowledgements

The authors would like to thank all our participants, Prince of Songkla University (PSU), the French Institute of Research for Development (IRD) and CIRAD for their support. Additionally, the authors are grateful to Khon Kaen University (KKU) and the Land Development Department (LDD) for all technical and logistical support, as well as the Department of National Parks, Wildlife and Plant Conservation (DNP), and all farmers for their helpfulness in the fields.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was funded by the Graduate School and Natural Rubber Innovation Research Institute, Prince of Songkla University (PSU) [grant number NAT590338S].

Notes on contributors

Phantip Panklang

Phantip Panklang is a Ph.D. student and an agricultural research officer who has 15 years of experience working with farmers in southern Thailand, mainly focusing on soil quality improvement and soil ecology.

Philippe Thaler

Philippe Thaler is a senior researcher who contributes mainly to the ‘plant-soil interactions and biogeochemical cycles’ research topic. His research focuses on the ecophysiology of the rubber tree and also leads multi-disciplinary studies on the sustainability of natural rubber production.

Alexis Thoumazeau

Alexis Thoumazeau is a researcher, working on assessing the performance of agroecological practices following multicriteria approaches to disentangle tradeoffs between agro-environmental and socio-economic aspects.

Rawee Chiarawipa

Rawee Chiarawipa is an assistant professor in the department of plant science, where he teaches plant physiology, agricultural meteorology, and tropical agricultural resource management. His research interests include rubber plantation management, plant ecophysiology, and the management of fruits and medicinal plants.

Sayan Sdoodee

Sayan Sdoodee is a professor with more than 40 years of teaching experience and research. His research has focused on crop physiology, especially the physiology of tropical fruits and rubber trees.

Alain Brauman

Alain Brauman is a senior scientist. His current work concerns the impact of land use changes on soil functioning and he recently conducted this research to better assess the relationship between soil functional assemblages and soil ecosystem services in agro-ecological systems.