Changes in soil chemical properties and their spatial distribution after logging and conversion to oil palm plantation in Sabah (Borneo)

Abstract

Conversion of primary forest into oil palm plantations is common in tropical countries, affecting soil properties, ecosystem services and land-use management. However, little is known about the short-range spatial soil distribution that is important for soil scientists, ecologists, entomologists, mycologists or microbiologists. In this study, seven soil properties (pH, EC (µS/m), P (mg/kg), NO₃^- (mg/kg), N%, C% and C:N) were measured to quantify the spatial autocorrelation across primary forest, selectively logged forest and oil palm plantation in Sabah, Malaysian Borneo. Local variograms were calculated (range ∼5 m) to determine the short-range variation, and a decision tree as well as principal component analysis were implemented to determine if the overall (global) mean differed between land uses. As hypothesised, oil palm soils deviated the most from primary forest soils, which had more fluctuating variograms and in general, a shorter range. Oil palm plantations also showed a difference in the global mean except for electrical conductivity. Selectively logged forests also differed in their short-range spatial structure; however, the global mean and variance remained similar to primary forest soil with the exception of labile phosphorus and nitrate. These results were attributed to initial plantation development, removal of topsoil, fertiliser application and topography.

Introduction

Anthropogenic land-use change often affects soil properties over a shorter time span than natural factors such as climate, topography, vegetation, age and parent material (Amundson and Jenny 1991; Winkler et al. 2021; Yaalon and Yaron 1966). In the tropics, logging of primary forests and conversion into agricultural land is widespread (Benhin 2006; Kleinschmit et al. 2021; World Resource Institute 1991), which reduces the quality and quantity of natural resources (Guerrero et al. 2020). This topic has been widely studied; however, little is known about the short-range variation of soil chemical properties and hence, biodiversity, ecology, entomology, microbiology, mycology and many other ecosystem services.

It is known that the conversion to oil palm plantation reduces numerous ecosystem services provided by soil. For example, after logging, areas that are leveled (along with roads) increase run-off causing erosion (Van Wambeke 1992; Hartanto et al. 2003), and with the removal of topsoil and drainage, the soil organic matter (SOM) content is reduced (Krejčová et al. 2021; Andriesse and Schelhaas 1987). This causes loss of fertility (Dressen et al. 1976) and an increase in soil compaction (Lal 2021; Lyczak et al. 2021). These losses further cause a reduction in biodiversity (World Resource Institute 1991), can affect soil processes such as bioturbation (Tuma et al. 2019), increase evapotranspiration (Uhl et al. 1981) and cause pollution (Henderson and Osbourne 2000; Koh and Wilcove 2007). This reduces the provisioning of ecosystem services such as water purification, food production, livelihood, climate regulation and biodiversity (FAO and UNEP 2020).

The soil properties that provide these functions and services are correlated in geographic space. This spatial autocorrelation reflects processes that help form soils and informs proper land-use management. The most common geostatistical technique to model the spatial structure of soil properties is known as a variogram. A variogram generalises a random process over the distance between point pairs (Matheron 1963). A variogram shows the short-range variance or measurement errors (nugget), the total semivariance (sill) and the distance over which the property is no longer correlated in geographic space and/or time (range).

Oil palm plantations need the least area per litre of any vegetable oil (Basiron 2007; Poore and Nemecek 2018), provide a low financial cost to the consumer and high profitability for the producer, and are a driver of economic development (Qaim et al. 2020). Consequently, oil palm expansion is driving deforestation in Southeast Asia (Henderson and Osborne 2000; Koh and Wilcove 2007; Palm 2005). For example, in Malaysian Borneo, from 1990 to 2018, oil palm was responsible for decreasing the area of primary forest by 13% in Sarawk and 16% in Sabah and decreasing the area of peat forests from by 21% in Sarawak and 19% in Sabah (Jaafar et al. 2020).

In the state of Sabah, Malaysian Borneo, oil palm plantations are generally established from areas that were originally selectively logged forests (Osman et al. 2012), where trees with the highest economic value were harvested first and then fast-growing trees in the following round of harvest. Once logged, the areas were left to naturally regenerate, often for multiple decades. However, even with regulations, logged forest mainly consists of pioneering trees and the forests were left degraded (Reynolds et al. 2011). These areas were converted to oil palm by mechanically levelling, draining and/or removing the topsoil as well as adding infrastructure, fertiliser, herbicides and pesticides.

The aim of this study was to determine if land use alters soil chemical properties. It is hypothesised that the largest change will be seen in the short-range spatial autocorrelation in oil palm and logged forest relative to primary forest soils. Additionally, it is believed that the global mean will be significantly different for oil palm but similar for logged forest and primary forest soils.

Materials and methods

Study sites

The research site (Figure 1a, b) forms part of the Stability of Altered Forest Ecosystems (SAFE) project in the state of Sabah, Borneo Malaysia (Ewers et al. 2011; Tuma et al. 2019) approximately 4 $^\circ $ 40’ 27” N and 117 $^\circ $ 31’ 40” E with an elevation ranging from 370 to 500 m. The climate is characterised by an isohyperthermic soil temperature and an udic soil moisture regime with a mean annual temperature of 26.7 and mean annual precipitation of ∼2,650 mm (Kumagai and Porporato 2012; Kuntashula et al. 2014; Walsh and Newbery 1999). The dominant soil types in the region are Typic/Humic Hapludults and Kandiudults on hilly positions with various types of Histosols, Inceptisols and Entisols in lower regions (Sakurai 1999). The geology is characterised by a metamorphic complex of amphibolites, and gneiss intruded by granite, gabbro, and ultramafic rock (Haile and Lee 1997).

Figure 1. Site within the surrounding area (a), sampling clusters (b) and sampling orientation in each cluster extending from the GPS placed at the black dot (c). Numbers represent how the coordinates were manually entered for each observation. PF: Primary Forest; LF: Logged Forest; OP: Oil Palm.

Three land uses were sampled before any SAFE project-related experimental fragmentation. The primary forest has been under conservation since 1976 and has experienced little to no human disturbances. This area lies on both eastern and western hill slopes and has a mean elevation of 431 m, and Shorea and Diospyros tree genera predominate. The logged forest areas were on western sloping positions, has a mean elevation of 450 m, had been selectively logged at least twice (between 1970–1990 and 1999–2010) and was previously logged for Dryobalanops, Dipterocarpus, Shorea and Parashorea tree genera. Only past signs of logging have been observed and secondary forest now predominates. Oil palm plantations were found almost continuously in slightly lower elevation areas (mean elevation of 385 m) such as foot slopes and valleys, planted in 2006 and are fertilised twice a year with diammonium phosphate ((NH₄)₂HPO₄)), potassium chloride (KCl), ammonium sulphate ((NH₄)₂,SO₄), magnesium sulphate (MgSO₄) and disodium tetraborate octahydrate (Na₂[B₄O₅(OH)₄]·8H₂O) (Elias et al. 2020). However, the amount of fertiliser applied was not known. The oil palm site is managed by the company Benta Wawasan Sdn. Bhd. (Ewers et al. 2011) and consist of Elaeis guineensis monocultures with a low, open canopy and sparse under vegetation.

Twenty-two sampling clusters (Figure 1b) were dug on the three land-use areas (8 in primary forest, 6 in logged forest and 8 in oil palm). Each sampling cluster comprised 10 samples (n = 220), 1 m apart, taken in the shape of a right triangle as shown in Figure 1c. Since the sample locations within each cluster were inside the precision of the GPS (5 m), the coordinates were manually entered.

Composite soil samples (average of total depth) down to 20 cm in depth were taken, were homogenised by hand (in the field by mixing the soil by hand), then sieved (2 mm), and oven-dried at 70–80℃ for 3–5 days. Seven soil chemical properties were measured: pH, EC (µS/m), P (mg/kg), NO₃ ^- (mg/kg), N%, C% and C:N. The pH and EC were measured in a 1:2.5 soil to water solution after mechanically shaking for an hour and filtered. A 1:100 soil to Mechlich-3 solution was shaken for 30 min and filtered to measure labile P with a UV Genesys 10 Thermo spectrometer at an 889 nm wavelength as described by Watanabe and Olsen (1965). To determine NO₃ ^-, samples were extracted in deionised water (1:5 soil:water ratio) and filtered, and the NO₃ ^- was determined using a colorimetric method (Zbíral et al. 1997). The same filtrate was used to determine pH and conductivity using a glass and potentiometric electrode. Total organic C and N percent were measured through a Flash 2000 Thermo Scientific elemental analyser.

Data analysis

A flow chart of the process used to determine the variability of the soil properties on the site is shown in Figure 2. From the observations, variograms were created for each soil property and the model fit was estimated through ordinary kriging. The difference in mean soil properties between land uses was then calculated through a conditional inference tree (ctree). Then, principal component analysis (PCA) was conducted to determine the variance. To determine the correlations between soil properties, a correlation matrix was constructed for all soil properties within a land use. All statistical analysis was conducted in R software (R Core Team 2017).

Figure 2. Flow chart of the process used to capture the variation of measured soil properties.

Spatial autocorrelation

Before developing the variograms, data for each soil property were checked for normality and heterogeneity through a Shapiro–Wilk test (Shapiro and Wilk 1965) and Levene’s test (Levenes 1960), respectively. Variograms were then constructed to determine the spatial structure of individual soil properties within each land use. Variograms are derived from the semivariance between point pairs with distance. Mathematically, this is described as follows (Matheron 1963):

$$\gamma \left( h \right) = {1 \over {2\left| {N\left( h \right)} \right|}}\sum\limits_{N\left( h \right)} {{{({z_i} - {z_j})}^2}}, $$

where $\gamma \left( h \right)$ is the semivariance at distance $h$ , $N\left( h \right)$ is the set of all point pairs with distance, and ${z_i}$ and ${z_j}$ are values of the soil property at distance i and j, Therefore, $h$ is the product of i and j. This is known as an empirical variogram from which theoretical variograms are computed from.

Theoretical variograms were then fit using residual maximum likelihood (REML; Bartlett 1937) to automatically fit the nugget ( $h = 0$ ), the sill ( $h = \;\mathop {\lim }\limits_{h\; \to \infty } \gamma \left( h \right))$ and the range. This method was used because it is an unbiased estimate of variance, it does not need homogenous variance, its residuals can be correlated and it is computationally efficient (Johnston 1972). Four theoretical variogram models were tried for each property including spherical (Sph), exponential (Exp), gaussian (Gau) and wave (Wav) models. A more detailed explanation of these variogram models can be found in the supplementary material A1. The final models were selected through ordinary kriging (Matheron 1963) with leave-one-out cross-validation. The model with the lowest root mean squared error (RMSE) for each land use was selected as the final model.

To compare soil properties on different land uses, the spatial autocorrelation of each variogram was evaluated on their spatial dependence, which involved the ratio of the nugget, sill and range (NSR). In this paper, the NSR is as follows:

$${\rm{NSR}} = {{{\rm{Nugget}}} \over {\left( {{\rm{Sill}} + {\rm{Range}}} \right)}},$$

The nugget, sill and range were scaled from 0.1 to 0.9, so each variable had equal weight (e.g., no large range or sill). The smaller the NSR is, the greater the spatial autocorrelation. This spatial dependency measure was used because it incorporates the three common aspects of a variogram into one simple equation (x and y axis), where the commonly used nugget:sill ratio to compare variograms leaves out the range.

Mean soil properties

The means of the soil properties were evaluated between each land use through a type of decision tree known as a ctree (Hothorn et al. 2006; Zeileis and Hothorn 2015). A ctree recursively splits the data on the independent variable to make the dependent variable more homogenous. It does so until no more splits are possible in the ctree, or a user-defined significance value (p-value) has been reached (Hastie et al. 2009). The ctree was developed with a Bonferroni test at each split with a p-value of 0.01 as the stopping criteria. A ctree was used because they are easy to interpret, do not need linear data, perform well with small sample sizes and are non-parametric (Hothorn et al. 2006; Zeileis and Hothorn 2015).

After the ctree was grown, the residuals for each model were checked for normality, heterogeneity and spatial autocorrelation. If there was spatial autocorrelation, the residuals were kriged and added back to the model, and the distributions were checked again. It should be noted that the models were only used to see the statistically significant differences and not to predict the soil properties.

Soil property correlations

Before PCA analysis (Pearson 1901), the soil properties were normalised $\left( {\mathop \to \limits_{0,1} } \right)$ so that all properties had an equal weight in the covariance matrix. The soil properties were evaluated on how much each contributed to the variability and how well the soils represent the principal components (PCs), also known as Cos². The contribution is the percentage a soil property contributes to the total variance in a PC and ranges from 0 to 100%. The further away the soil property is from the origin (coordinates = 0,0), the more that property contributes to the variance explained. A larger Cos² implies a better representation within a PC and, hence, the more important the property is. Cos² is simply the sum of the squared coordinates of a soil property on a biplot.

A correlation matrix was created to determine how the soil properties correlate with each other within a land use. Correlations were determined using the REML model as described in Sec 2.2.1. p-Values and Pearson’s correlation coefficient (R²) were used to interpret the correlations.

Results

Spatial autocorrelation

The variograms for each land use (Figure 3) and (Table 1) clearly show a large difference in the semivariance for each land use. In primary forest soils, pH (nugget = 0.00, sill = 0.02, range = 0.03 m, NSR = 0.50), EC (nugget = 148, sill = 480, range = 5.85 m, NSR = 0.24), N (nugget = 0.00, sill = 0.008, range = 1.11 m, NSR = 0.29) and C (nugget = 0.50, sill = 0.80, range = 1.50 m, NSR = 0.12) all had Exp theoretical variograms with RMSE of 0.15, 15.2 µS/mg, 0.08% and 1.01%, respectively. On the other hand, labile P (nugget = 2.00, sill = 5.50, range = 3.00 m, NSR = 0.17), NO₃ ^- (nugget = 2.00, sill = 5.50, range = 3.00, NSR = 0.17) and C:N (nugget = 0.00, sill = 1.56, range= 1.11, NSR = 0.29) had a Sph theoretical variogram as their best-fitting model with a RMSE of 2.10 mg/kg, 8.37 mg/kg and 1.31, respectively.

Figure 3. Variograms with the trends for primary forest (PF), logged forest (LF) and oil palm (OP) being green, blue and yellow, respectively. The points represent the semivariance between observation pairs, while the curve is the model of the observation pairs.

Table 1. Theoretical variogram models used, goodness of fit (RMSE) and parameters (nugget, sill, range and normalised spatial dependency (NSR) for primary forest, logged forest and oil palm for each soil property.

For models, “Exp” = exponential, “Sph” = spherical, “Gau” = Gaussian and “Wav” = wave models.

Selectively logged forest changed to more complex theoretical variogram models with only labile P (nugget = 2.00, sill = 5.00, range = 1.50 m and NSR = 0.25) having an Exp variogram with a RMSE of 21.1 mg/kg, while only NO₃ ^- (nugget = 37.9, sill = 79.0, range = 2.80 m and NSR = 1.26) and N (nugget = 4.00×10^-3, sill = 9.00×10^-3, range = 3.80 m, NSR = 14) had Sph theoretical variograms with a RMSE of 8.50 mg/kg for NO₃ ^- and 0.09% for N. Properties with Wav theoretical variograms included pH (nugget = 0.04, sill = 0.10, range = 1.70 m and NSR =0.24) and EC (nugget = 198, sill = 245, range = 1.60 and NSR = 0.63) with RMSE of 0.33 and 16.0 µS/mg, respectively. Carbon (nugget = 1.00, sill = 3.00, range = 2.80 m and NSR = 0.17) and C:N (nugget = 1.00, sill = 3.90, range = 3.20 m and NSR = 0.16) had a Gau theoretical variogram with a RMSE of 1.34% for C and 1.25 for C:N.

In oil palm, the Wav theoretical variogram was dominant with the best-fitting models for labile P (nugget = 200, sill = 400, range = 1.50 m and NSR = 0.49), NO₃ ^- (nugget = 19.4, sill = 27.8, range = 0.88 m and NSR = 0.39), N (nugget = 3.00×10^-4, sill = 9.00×10^-4, range = 1.12 m and NSRE = 0.29) and C (nugget = 0.05, sill = 0.13, range = 1.55 m and NSR = 0.25) with RMSE of 21.9 mg/kg, 5.38 mg/kg, 0.03% and 0.34%, respectively. The Gau model was utilised for pH (nugget = 0.05, sill = 0.06, range = 1.70 m and NSR = 0.23) and C:N (nugget = 1.00, sill = 3.00, range = 4.00 m and NSR = 0.23), which had a Gau theoretical variogram and a RMSE of 0.32 for pH and 1.20 for C:N. Electrical conductivity (nugget = 200, sill = 400, range = 1.00 m and NSR = 0.52) was the only property to have an Exp model and had a RMSE of 20.3 µS/mg.

Mean soil properties

Except for NO₃ ^- and labile P, there was no statistically significant difference between the mean in primary forest and logged forest soil properties (Figure 4) according to the ctrees’ terminal nodes. On the other hand, oil palm soils were statistically different from both primary forest and logged forest for all soil properties except EC. Oil palm plantations had a lower pH (5.25) relative to primary forest (6.00) and logged forest (5.98). Electrical conductivity was the lowest in oil palm (57.5 µS/cm), then logged forest (67.4 µS/cm) and highest in primary forest soils (70.5 µS/cm). Labile P was the lowest in primary forest (5.70 mg/kg), then logged forest (6.61 mg/kg) and highest in oil palm soils (29.1 mg/kg), which was also seen in the PCA. Oil palm plantation soils had low NO₃ ^- (10.4 mg/kg), primary forest soils had relatively high (21.4 mg/kg) and logged forest soils had the highest NO₃ ^- (26.6 mg/kg). Nitrogen was highest in logged forest (0.28%) and primary forest (0.26%) with the lowest in oil palm (0.13%). Carbon follows a similar trend, with logged forest soils having the greatest (3.41%), followed by primary forest (3.12%) and then oil palm with the lowest (1.25%). The C:N ratio indicates that primary forest has the highest (12.3), followed by logged forest (12.0) and then oil palm (9.27) soils.

Figure 4. Boxplot of the significant differences (p<0.01*) of soil properties and land-use types according to the conditional inference trees’ terminal nodes. The horizontal line represents the median, the hinges represent the 25^th and 75^th quantiles, the vertical lines represent the distance between the 1^st and 3^rd quartiles, and the points are outliers. PF: Primary Forest; LF: Logged Forest; OP: Oil Palm.

Soil property correlations

The first two PCs accounted for 64% of the variability in the dataset (Figure 5). Nitrogen and C percent contributed the most to PC1 (darker lines) with 23% and 26% contribution, respectively. These two elements are also best represented in PC1 with a Cos² of 0.75 and 0.84 for N and C, respectively. Nitrate and EC represent the next largest contributors to PC1 with contributions of 16% and 14% with a Cos² of 0.50 and 0.44, respectively. All four of these elements were negatively correlated with PC1 and PC2. In other words, as the two PCs decrease, EC, NO₃ ^-, N and C values increase.

Figure 5. The first two principal components of all soil properties with their contributions to the components (colour of arrows in percent). Ellipsis show the clusters for primary forest (PF), logged forest (LF) and oil palm (OP). The larger the ellipsoid, the more variation in soil properties determined by the observations (points) for the principal components.

Labile P represented more than half the contribution to PC2 (58%) with a Cos² of 0.74. The pH also contributed to PC2 (26%) but had a much lower Cos² (0.32). Unlike the other soil properties, a high P content was correlated with oil palm and had little correlation to primary forest or logged forest. On the other hand, higher values for all other soil properties were more strongly correlated with primary forest and logged forest. This can be seen in the strong directionality for both C and N towards these two land uses.

The R² values (Figure 6) for primary forest soils show that EC was highly correlated with P (0.66), N (0.82) and C (0.82). Labile P correlated with N (0.74) and C (0.74), while N correlated highly with C (0.93). In logged forest soils, EC correlated with N (0.86), C (0.84) and C:N (0.74), while N correlated with C (0.97) and C:N (0.75). Surprisingly, pH negatively correlated with N (-0.69) and C (-0.49) in OP soils and again, N was highly correlated with C (0.91).

Figure 6. Correlation matrix of soil properties within primary forest (PF), logged forest (LF) and oil palm (OP). Blank spaces indicate p-values > 0.01, and therefore, the correlations are not displayed.

Discussion

Logging altered the short-range spatial structure of the soil properties, yet the mean values of most of the soil properties remained similar between primary forest and logged forest soils. This is most notable for pH, where it was observed that the theoretical variogram changed from an Exp model with a high NSR in primary forest soils to a Wav model with a small NSR in logged forest soils. Therefore, the spatial structure changed to a higher degree of spatial dependency that has positive and negative correlations with an increase in range (Mahdi et al. 2020). This agrees with observations made by Lima et al. (2020), that logging substantially alters spatial patterns of many ecosystem parameters, for example, logging could have made the forest less homogenous and/or opened the canopy (Berry et al. 2008). There should be noted that spatial aspects of logged forest is more variable than primary forest; however, the larger differences between logged forest and primary forest in smaller scale variation appear even in situations where aspects is similar. Despite some effects by slope variations and these cannot be excluded, it is believed that selective logging was a major reason for the observed difference between primary forest and logged forest.

Interestingly, EC and NO₃ ^- in primary forest had effective ranges greater than the distance than the soil samples and, therefore, increased without bounds at this scale. This may indicate there is trend in the data such as a slope gradient and sample design, or the scale was too fine (Webster and Oliver 2001). However, this was not found for properties in any other land use, and because of the short distances, it is thought this was due to the fine scale. However, there are many factors that can affect EC such as mineralogy, salts, structure, water content, bulk density and more (Adviento-Borbe et al. 2006). Therefore, further studies into soil physical properties could improve the interpretation of EC and NO₃ ^- on the different land uses.

Labile P had a much greater sill in oil palm than primary forest and logged forest soils. This is not surprising for two reasons: i) P generally binds to sorption sites on minerals and, therefore, gets preferentially eroded causing variability, and ii) the plantation applies P fertiliser, which may have not been precisely applied causing the Wav spatial structure and increasing the variance at this scale.

It was speculated that the low spatial autocorrelation could be due to an unpronounced catena effect so that the logged forest soil properties do not influence each other as much as in primary forest soils (smaller nugget to sill ratio with shorter range). This can make the soil distribution appear random, which has implications for soil mapping in the area. For example, randomness would make field surveys more expensive, or the maps could be less accurate (Beckett and Burrough 1971). This has downstream effects particularly for land-use management, restoration efforts and monitoring environmental resources. However, it is difficult to know if this is a result of other environmental factors such as topography, climate, and parent material. For example, Flynn et al. (2020) found that both physical and chemical factors are strongly affected by landform elements.

Oil palm also shows a difference in the spatial structure of the variograms compared to primary forest. Like logged forest soil, oil palm seem to have changed into a more periodic structure (semivariance fluctuates around the sill over distance), especially for components commonly associated with SOM such as P, N, C and C:N. It is thought that the stripping of SOM during plantation development was the cause of this trend. When stripped, some places are easier to strip than others and may not always be near to each other and, hence, the periodic structure. Therefore, in this area, logged forest and oil palm change the spatial distribution of these seven soil chemical properties, which is important for economic efficiency and sustainability.

Unlike the variogram analysis, primary forest and logged forest soils are similar according to the global mean and variance. Therefore, it is believed that logged forest soils experienced little degradation and/or have partially regenerated. Both land uses show a similar distribution in orthogonally transformed space, indicating that the same soil properties are controlling the global variance. Not surprisingly, N and C account for the most variation in primary forest and logged forest as these land uses were not stripped of SOM. Originally, expected by general patterns of C and N, primary forest would have the highest amount of both C and N; however, trees in the Dipterocarpaceae family associate with ectomycorrhiza (Smits 1994), which help to form microaggregates and root exudes that protect N and C physically and chemically (Center for International Forestry Research 1998). Therefore, N and C were the same in primary forest and logged forest soils even with the more open canopy due to soil stability.

It was unexpected that NO₃ ^- and P were statistically higher in logged forest than primary forest as it was thought that the more open canopy of logged forest would cause more surface erosion or leaching of salts. This could also indicate that conditions for mineralisation are greater in logged forest soils possibly due to better drainage (Stottlemyer et al. 2001) and higher temperatures (Miller and Geisseler 2018), and the litter fall from the pioneering trees could be easier for microbes to break down (Quan et al. 2014). Since logged forest soils have a southeastern slope and a more open canopy, these attributes may have better drainage and higher temperatures and, hence, caused greater mineralisation. Unfortunately, this has an environmental impact as it can cause eutrophication both on- and off-site (Carpenter et al. 1998).

On the other hand, oil palm soils were not similar in the transformed space and besides P, soil properties had a lower global mean and variance in the oil palm cluster. Nevertheless, in the tropics with high temperatures and precipitation, low pH values were expected. The lower pH in oil palm soils is most likely due to nitrification of NH₄ ⁺-based fertiliser application and the subsequent hydrolysis of Al³⁺ both of which release H⁺ (Ragland and Coleman 1960; Sigurdarson et al. 2018). These changes in soil properties have impacts on ecosystem services dramatically such as long-term degradation causing a reduction in growth and yield. For example, it affects nutrient cycles and availability and nutrient toxicity (flora and fauna), decreases water quality, decreases soil fertility, facilitates leaching (ion exchange) and influences SOM accumulation (Lükewille and Alewell 2008).

The results of EC in oil palm were unexpected because it was thought that the application of fertiliser and the position in the landscape (low lying) would cause salt and clay accumulation. However, a pH below 5.50 facilitates the exchange of Al⁺³ and H⁺ for basic cations on the solid phase and, hence, leaching of salts (Haynes & Swift 1986). This reaction would be pushed further with the addition of fertiliser (acidification) and leaching (cations removed). Additionally, this could represent a low water-holding capacity due to the removal of the C-rich topsoil (Brevik et al. 2006). Nevertheless, the global mean EC was similar in all land types, further indicating that there are more processes controlling EC than simply land use. These could include many soil physio-chemical, environmental and anthropogenic factors such as geology, water quality, atmospheric deposition, land management, landscape position amongst others (Corwin and Lesch 2005) and bioturbation/biogenic processes (Tuma et al. 2019).

Not surprisingly, P in oil palm was much greater than in primary forest and logged forest. For oil palm soils, it is speculated that the application of diammonium phosphate increased the P content along planting rows. The high P content in oil palm soils can have a large ecological impact as P often ends up in aquatic systems (Smith et al. 1999) creating hypoxic conditions in lakes and coastal zones, decreases water quality (Khan and Mohammad 2014) and decreases biodiversity (Rabalais 2002). This is of concern as some measurements of P were excessively high, increasing the risk of land degradation.

It was thought that fertiliser in oil palm soils would increase the NO₃ ^- concentration due to the nitrification of NH₄ ⁺ fertiliser; however, this was not the case. The most likely explanations for the low NO₃ ^- are increased leaching from the profile due to the udic soil moisture regime (Spalding et al. 2001), an increase in erosion from heavy rain with less protection from the canopy or the lower pH has been slowing the nitrification of NH₄ ⁺ to NO₃ ^- (Amatya et al. 2011) caused by an alteration of microbe populations (Lee-Cruz et al. 2013; Tripathi et al. 2016). Conversely, in primary forest and logged forest soils, NO₃ ^- was more abundant most likely due to vegetation cover and possibly, and nitrification was not inhibited on these land uses.

Like the spatial distribution of oil palm soils, the loss of N and C most likely originated from the initial plantation development due to oxidation of SOM, anthropogenic pedoturbation and the physical removal of topsoil (Corley and Tinker 2015). Additionally, oil palm do not naturally form symbiotic relationships with mycorrhiza and need to be inoculated with arbuscular mycorrhiza (Sundram 2010). Therefore, N and C do not have the protection that is provided by the fungi in primary forest and logged forest soils. However, C losses were not observed by Khasanah et al. (2015) over long-term monitoring of C in oil palm soils (> 25 years). This is because initially, C declines and then increases from oil palm residues, forming a U-shaped curve. Since the oil palm in this study site is still considered early in the life cycle, C and N may be on the downward trend of the U curve described by the authors. Additionally, the lower C:N of oil palm soils would indicate faster mineralisation rates and the release of NH₄ ⁺ with subsequent uptake from organisms or leaching and loss of C through the release of CO₂ (Hossain et al. 2017).

An additional pattern that emerged was the correlation between properties in each land use, which helped to further explain these results. The difference in property interactions indicates a difference in soil functionality. For example, the negative correlation of N and C with pH in oil palm soils (positive in primary forest and logged forest) could represent a loss of the pH buffering capacity due to low SOM from removal and heavy fertilisation. The loss of the pH buffering capacity due to N fertiliser is poorly understood; however, Zhang et al. (2016) found that even on calcareous soils with a high buffering capacity, fertiliser decreased this capacity and soils acidified. Alternatively, this negative correlation could be due to the acidic soils stabilising the C content, all be it low, through chelation with sesquioxides (Sparks 1995). Therefore, as the pH decreases, SOM content could increase in these soils due to the chemical protection chelation provides.

Electrical conductivity correlated with P in primary forest soils and is best explained by the negative charge of labile orthophosphate (PO₄ ³-), which together with cations will facilitate a higher EC as it has been shown that the concentration of cations and anions positively affects EC (Carmo et al. 2016). Electrical conductivity also correlated with N and C in both primary forest and logged forest soils. These two soil properties make up a large portion of SOM with a pH-dependent charge and, therefore, can hold more cations and water, which would increase the EC as these properties increase. However, Auerswald et al. (2001) found that SOM only had a small influence on EC, which was more correlated with clay, fertilisation history and water content (near wilting point). Additionally, these results are contradictory to Lei et al. (2019), who found a negative correlation of EC with C. However, the greater SOM content can also cover P sorption sites on clay minerals and the greater SOM stock can be mineralised to release P, thereby increasing EC (Yusran 2011).

As expected, C and N were highly correlated in all land uses. However, it was surprising that N and C correlated with C:N only in logged forest soils. This was unexpected because it was thought that both N and C would be highly correlated with C:N in all land uses due to the fact the chemicals make up the ratio. Additionally, N had a positive correlation with C:N, but it was expected that N would be negatively correlated with C:N. Therefore, it is thought that if N and C were accumulating, N does so at a slower rate in logged forest soils. This hypothesis was confirmed by the coefficients of N (0.06) and C (1.16) in the linear regression. The difference in the accumulation rates could be due to the logging process and vegetation type, which can leave more lignin and, hence, more C on the soil surface.

Conclusion

This study showed the effect land use has in a tropical landscape and potential depletion of ecosystem services provided by soils. The spatial distribution of soil properties differed on each land use with logged forest showing less spatial autocorrelation and oil palm showing a different spatial structure relative to primary forest. Although the spatial distribution of soil properties changed in logged forest, the soils appear to have little to no degradation, or they have regenerated as the mean values of the soil properties were similar to those of primary forest. On the other hand, the mean soil properties changed in oil palm relative to both primary forest and logged forest, which were associated with degradation of ecosystem services. Correlations of soil properties within each land use also differed between land uses, and the results indicated a loss of pH buffering capacity in oil palm soils. These changes suggested SOM degradation from mechanical removal and management practices such as fertiliser and pesticide application in oil palm soils resulting in threats to biodiversity, soil and water quality and sustainability of oil palm agriculture. Further research should include measuring soil physical properties to better understand the dynamics of the conversion of ecosystems for different land uses.