Tropical forest conversion to rubber plantation affects soil micro- & mesofaunal community & diversity

Tropical rainforests play important roles in carbon sequestration and are hot spots for biodiversity. Tropical forests are being replaced by rubber (Hevea brasiliensis) plantations, causing widespread concern of a crash in biodiversity. Such changes in aboveground vegetation might have stronger impacts on belowground biodiversity. We studied tropical rainforest fragments and derived rubber plantations at a network of sites in Xishuangbanna, China, hypothesizing a major decrease in diversity with conversion to plantations. We used metabarcoding of the 18S rRNA gene and recovered 2313 OTUs, with a total of 449 OTUs shared between the two land-use types. The most abundant phyla detected were Annelida (66.4% reads) followed by arthropods (15.5% reads) and nematodes (8.9% reads). Of these, only annelids were significantly more abundant in rubber plantation. Taken together, α- and β-diversity were significantly higher in forest than rubber plantation. Soil pH and spatial distance explained a significant portion of the variability in phylogenetic community structure for both land-use types. Community assembly was primarily influenced by stochastic processes. Overall it appears that forest replacement by rubber plantation results in an overall loss and extensive replacement of soil micro- and mesofaunal biodiversity, which should be regarded as an additional aspect of the impact of forest conversion.

stages of some metazoan species 23 . These organisms have been often neglected from soil biodiversity surveys even though they are considered to play an important role in land ecosystem functioning. For example, nematodes play important roles in soil ecological processes such as nutrient cycling and plant growth 24,25 and are sensitive enough to environmental changes to be used as bioindicators for ecotoxicological assessment of soil 26 .
In the past few years, DNA sequencing-based analyses of soil and sediment micro-and mesofaunal communities have revealed divergent lineages [27][28][29] , accompanied with large datasets of unknown sequences within known phyla that hint at a vast but poorly known diversity. Amongst Nematoda, for instance, it has been estimated that <4% of total species diversity is formally described 30 . This information gap is principally due to difficulty in studying these organisms because of their small size and hard-to-distinguish morphology at the species taxonomic level.
In this study, we focus on rainforests in the global biodiversity hotspot of the Xishuangbanna region of Yunnan, SW China. In Xishuangbanna, the original primary forest, as well as the derived logged forest, has been converted rapidly to rubber (Hevea brasiliensis Müll. Arg.) plantation. In 1992, rubber covered 87,000 ha, rising to 153,000 ha in 2002 and 424,000 ha in 2012 16 . The impact of forest conversion to rubber plantation in southern China has already been explored with respect to soil nematodes, which is one aspect of micro-and mesofaunal diversity. Using traditional morphological methods, Xiao et al. 31 compared nematode diversity in natural forests and rubber plantations, finding that conversion was associated with a decrease in both αand βdiversity.
The present study goes further in encompassing the soil micro-and mesofauna and by including a larger number of sites than was sampled by Xiao et al. 31 . The present study also utilizes high-throughput sequencing interpreted with the help of biologically informative taxonomic resources covering alternate loci such as the 18S rRNA gene 32 . As such the technique is capable of recognizing cryptic biodiversity even within relatively well-studied groups (e.g., Nematoda), as well as undocumented biodiversity in relatively poorly studied groups (e.g., Tardigrada).
Our main hypotheses in this study were as follows: (1) We expected that forest conversion to rubber plantation would significantly alter the soil community composition, with each land-use type harboring a distinct soil metazoan community, with the forest conversion driving changes in soil chemistry which in turn modified the soil biota. (2) We predicted lower soil micro-and mesofaunal αdiversity in rubber plantation.
(3) We also predicted lower βdiversity of micro-& mesofauna within rubber plantation compared to the native forest. (4) We predicted that the micro-& mesofaunal community assembly in both the native forest and derived rubber plantations would be strongly influenced by deterministic rather than stochastic/neutral processes, just as is the case for soil and sediment microbes from similar environments [33][34][35] .

Materials and Methods
Study location. The

Study area description.
Xishuangbanna has a seasonal outer tropical climate, with hot, rainy summers and somewhat cooler and drier winters, which remain free of frost 36 . The tropical rainforests of Xishuangbanna are characterized by high species richness with up to 4180 species of vascular plants/2500 m 2 plot and prevalence of epiphytes and epiphyllous mosses on the trees with the ubiquitous presence of woody lianas and strangling plants 37 . The tropical rainforest can be classified into two subtypes: tropical seasonal rain forest (in the lowlands; <900 m a.s.l., or meters above sea level) and tropical montane rain forest (900-1600 m a.s.l.) 38 , but with soil type differentiation within each of these (for soil types see Supplementary Table S1). For example, some areas of the lowland and montane forest are on limestone soils, with a distinctive plant community.
Field sampling. Sampling was carried out during the month of August 2013. Seven rainforest sampling sites were chosen within a 10 km radius of XTBG. Seven rubber plantation sites were chosen adjacent to each old-growth forest sampling site. In this study, there are fourteen sites in total (seven forests and seven rubber). At each sampling site, four quadrats (10 m × 10 m), were located >30 m apart along a linear transect. See SI Appendix Fig. S1 for sampling design; see Supplementary Table S1 for more details. Each 10 m × 10 m quadrat provided a single sample. This individual sample was the result of combining and mixing soil from 5 equal subsamples (each approximately 50 g) of soil from the top 10 cm, sampled using a small trowel marked to 10 cm depth starting down from the base of the litter layer. The five subsamples, one taken at each corner and one at the center, were gathered from each 0.01 ha area and mixed into a single soil sample bag 39 . Soil sample processing and DNA isolation were done one day after soil sampling. A list of the sampling sites and data on environmental parameters is provided as Supplementary Table S1. This study led to a total of 56 samples (7 sites × 4 biological replicate samples per site = 28 forest, and similarly 7 × 4 = 28 rubber plantation samples).
Soil micro-and mesofaunal DNA extraction and PCR amplification. To maximize capture of soil microfauna and mesofauna, 200 g of each soil sample was suspended in around 2 I distilled water (DW) and then gently sieved through a set of sieves. For the combined micro-& mesofauna isolation under the Baermann funnel method, we used three different sieves in a vertical series. We placed a 2 mm sieve at the top and a 150 µm sieve below it. Below this was a 20 µm www.nature.com/scientificreports www.nature.com/scientificreports/ sieve (coarse to fine from top to bottom). 2 I of soil suspension was slowly poured onto the 2 mm sieve. We gently washed down all of the soil particles that would pass through the 2 mm sieve using more distilled water, so that stones, roots, and litter, etc., remained on the 2 mm sieve. We then collected accumulated particles -containing fauna -from both the remaining sieves below (150 µm & 20 µm) using a wash bottle to wash through finer particles. Mesofauna would be selectively accumulated on top of the 150 µm sieve, and microfauna on top of the 20 µm sieve. The accumulated material from each sieve was removed by spatula and placed together onto cotton gauze in a Baermann funnel for further extraction followed by sugar flotation 40 .
As well as removing stones, roots, etc., the 2 mm sieve ensured that the macrofauna (>2 mm) did not end up in the collected material. Smaller eukaryotes (e.g., smaller protists) and clays, along with water, would be able to pass through the 20 µm holes of the finest sieve.
In this manner, we could collect both the micro-& mesofauna from the soil while keeping the contamination from macrofauna to a minimum. This collected soil was then used for both Baermann funnel and sugar flotation extraction. The Baermann Funnel method is most often used to capture Nematoda, but will also capture most other types of motile micro-& mesofaunal populations, as it depends on organismal locomotion 40,41 .
After 24-36 hours, fauna that had moved down out of the soil-derived material, to the water at the bottom of the funnel were collected by draining the fluid into a Falcon tube. Tubes were then centrifuged for 10 min approximately at 5000 rpm (settling micro-& mesofaunal populations at the bottom), and the individual sediment pellet was collected from each for DNA extraction 42 . Less active/dead components remaining in the soil material within the funnel were then captured by sugar flotation and centrifugation using a 40% (w/v) sugar solution, to yield a pellet in each tube 43 . This combined procedure was intended to increase the efficiency of capture of a range of micro-and mesofauna by capturing both active (funnel) and less active or dead components (flotation).
The extracted pellets from both methods were combined before DNA extraction at the tropical forest ecology lab in XTBG, China, using the MoBio Power Soil DNA isolation kit (MoBio Laboratories, Carlsbad, CA) according to the manufacturer's instructions.
The isolated DNA was stored at −80 °C, and was later used as a template to amplify a ~400 bp diagnostic region, defined by primers NF1 (C. elegans 1226-1250 bp position) and 18Sr2b (C. elegans 1567-1588 bp position) towards the 3′ end of the 18S rDNA with PCR reaction conditions as described by Porazinska et al. 44 . Detailed PCR procedure is described in SI Appendix, Section 1b. This primer pair was originally designed for nematodes, but we were able to test its suitability in coverage for the whole of subkingdom Metazoa. We used the Ecotaxaspecificity module (http://pythonhosted.org/OBITools/scripts/ecotaxspecificity.html) under the eco-Primers platform (http://pythonhosted.org/OBITools/scripts/ecoPrimers.html) which evaluates barcode resolution at different taxonomic ranks 45 . For the primer pair tested, ecoPrimers calculates a taxonomic coverage, which is the number of amplified target species relative to the total number of target species in the input database. Taxonomic coverage rates of 100% (SI Appendix: Section 1c,d) across all metazoan sequences for phyla in the SILVA SSU database suggest that the barcode designed for nematodes is degenerate enough to be used as a general primer for subkingdom Metazoa. All metazoan sequences across all phyla in the SILVA SSU database are predicted to amplify with this primer pair (Complete details can be found in SI Appendix: Section 1c,d).
All sequences shorter than 150 nucleotides, with homopolymers longer than eight nucleotides and all reads containing ambiguous base calls or incorrect primer sequences, were removed. Remaining sequences were aligned against the SILVA SSU eukaryotic aligned database. The sequences were further filtered to remove gaps generated after alignment, then preclustered using the Mothur implementation of pseudo-single linkage preclustering algorithm from Huse and colleagues 47 . Putative chimeric sequences were detected and removed via the Chimera Uchime algorithm contained within Mothur 48 in de novo mode, which first splits sequences into groups and then checks each sequence within a group using the more abundant groups as a reference. Taxonomic classification of Metazoa was performed at a Bayesian cut-off of 55% with 1000 iterations against SILVA-ARB (database containing aligned 18S ribosomal RNA sequences with a minimum length of 1200 bases for Eukarya~ http://www. arb-silva.de/download/arb-files/; August 23, 2013). The remaining reads were then clustered using Mothur's average algorithm. All samples were standardized by random subsampling to 1885 reads per sample (based on the smallest sample size by default), using the sub.sample command (http://www.mothur.org/wiki/Sub.sample), for calculating richness (OTUs at ≥99% sequence similarity), rarefaction curves, diversity and community compositional indices and matrices to be used later for the statistical analyses.
Statistical analyses. All data were analyzed using R 3.3.0 49 except for NMDS analysis which was performed in Primer-6 Software. To test for differences in richness (OTUs), abundance (read counts), diversity indices (Shannon, and Faith Phylogenetic diversity; PD) and soil physicochemical parameters i.e., TOC (Total Organic Carbon), TN (Total Nitrogen), pH, SX (Soil Texture), ST (Soil Temperature), Ele (Elevation), GSW (Gravimetric Soil Water) & AP (Available Phosphorus), (SI Appendix: Section 1a) between forest and rubber plantation, we used GLMM with a Poisson distribution in package lme4 & nlme and site as a random term (for integer dataset: OTUs) and glmmPQL with a Quasi-Poisson distribution in package MASS as this allows quasi distributions with cluster as a blocking term (for non-integer data-set: abundance, Shannon, and PD etc.) 50 . The functions glmm and lmer do not work with a quasi-Poisson distribution 51 .
We employed the methods of Jost 52 to multiplicatively partition γ-diversity into independent αand βcomponents for each management type. We used the "multipart" function in Vegan package and ran 1000 simulations to compare values against a null model. Following Jost 52 , we partitioned species richness using a multiplicative framework, in which γ = α × β. We assessed α-& βdiversity at three spatial scales, such that overall γdiversity within each forest was expressed as: ( where α.1 & β.1 denote species richness and βdiversity between sites, respectively; α.2 & β.2 denote species richness and βdiversity between forest and rubber plantation, respectively; α.3 & β.3, denote species richness and βdiversity within sites, respectively. We then assessed whether forest clearance and conversion to plantation had influenced the βdiversity, by partitioning γdiversity, separately for the forest samples and samples from rubber plantation, respectively. We performed Non-metric Multi Dimensional Scaling (NMDS) to explore patterns in phylogenetic community composition using both the unweighted UniFrac & β-MNTD (β-mean nearest taxon distance) distance matrix with Primer-6 software 53 . The unweighted UniFrac measures the distance between two communities by calculating the fraction of the branch length in a phylogenetic tree that leads to descendants in either, but not both, of the two communities 54 . β-MNTD is the mean phylogenetic distance to the closest relative in a paired community for all taxa and is sensitive to the changes of lineages close to the phylogenetic tips 55 . We used an analysis of similarity (ANOSIM) with 999 permutations to test if community composition was significantly different between the forest and samples from rubber plantation. To test whether the taxonomic composition results may have been influenced by pseudoreplication or spatial arrangement of study sites (Ramage et al. 2013), we used a Mantel test under the function mantel (999 permutations/analysis) in package Vegan in R 56 to compare the taxonomic composition to geographic distance between pairs of transects within a land-use type.
We performed a redundancy analysis (RDA) based variation partitioning analysis 57 to assess the relative effects of environmental and spatial variables on community composition. We used Hellinger transformed OTU abundance data as the response variable and two sets of explanatory variables which included environmental variables (pH, Ele, TOC, AP, GSW, SX, TN, and ST) and spatial variables (geographical co-ordinates for sampling sites), respectively. Before the RDA, the environmental variables with high variance inflation factor (VIF) >10 were eliminated to avoid collinearity among factors. The importance of environmental and spatial variables in explaining species composition was determined by an RDA analysis using Monte Carlo permutation tests (999 unrestricted permutations) followed by forward selection to remove the non-significant variables from each of the explanatory sets.
To further evaluate the relative importance of each environmental variable and spatial distance on the community phylogenetic dissimilarity, we used a multiple regression on matrices (MRM) approach 58 . Before applying MRM to the dataset, we looked for redundant environmental factors using the VARCLUS procedure 59 in the Hmisc R package 60 and removed TN as it was highly correlated with TOC (Spearman's ρ2 = 0.80). In MRM, non-significant factors were removed sequentially, and the analysis was repeated until only significant factors were left in the model. Significance was tested by permutations (9999), and P-values of the two-tailed tests are reported for this analysis.
Phylogenetic analysis. Phylogenetic analyses were performed to answer whether the community assembly in both the native forest and derived rubber plantations is influenced by deterministic or neutral processes (stochastic). A maximum-likelihood tree was constructed using aligned 18S rRNA gene sequences of representative OTUs in FastTree 61 . To evaluate the phylogenetic signal, an environmental optimum for each OTU was calculated for each environmental variable as in Stegen et al. 62 . Between-OTU environmental optimum differences were calculated as Euclidean distances using optima for all the environmental variables. The correlation coefficients between differences in environmental optima and phylogenetic distances were measured using Mantel correlogram with 999 random permutations 34,35 .
To calculate the turnover in phylogenetic community composition, we calculated the βMNTD 55,63 using 'comdistnt' function (abundance.weighted = TRUE) of Picante R package 64 . Further, to investigate the influence of deterministic and stochastic processes on community assembly of soil micro-and mesofauna, a null modeling approach was implemented 34,35 . First, we calculated the β-nearest taxon index or βNTI, which is the difference between observed βMNTD and the null distribution (999 times) of βMNTD measured in units of its standard deviation. Values of βNTI < −2 and >+2 indicate significantly less and more than expected phylogenetic turnover (deterministic assembly), respectively. However, comparisons falling within the null distribution of βMNTD (|βNTI|<+2) indicate that the observed difference in community composition is not the result of deterministic selection and are instead attributable to stochastic assembly 63 .

Results
The 18S rRNA gene metagenomic data (SI Appendix Fig. S2) gives details of both taxonomic identity and relative abundances of reads within each taxonomic category. It would be unrealistic to try to assign relative abundances of individuals to each taxonomic category since even closely related species of Metazoa from the same class or phylum can differ dramatically in size -and the samples pick up a great diversity of Metazoan taxa. For multicellular organisms, the number of cells and thus the number of copies of phylogenetic markers (such as 18S rRNA here) per individual are influenced by their size and their biomass. Consequently, under the premise that there was no technical bias, the number of reads obtained after PCR and NGS will probably be more closely correlated to their fractional biomass than abundance in terms of individuals.
The NMDS results revealed that phylogenetic community composition of both land-use types, based on both unweighted UniFrac and β-MNTD, was driven by forest conversion to rubber plantation, with samples from rubber plantation forming separate clusters apart from forest samples ( Fig. 1; ANOSIM results, R = 0.3, ρ = 0.001). These results show clearly that the micro-and mesofaunal community composition of the rubber plantation differed from the original native forest. We did not find any significant effect of distance on taxonomic composition, by testing for phylogenetic signal using Mantel correlograms of soil microbiome within land-use types for forest samples (Mantel test, P > 0.05). However, we did find a significant effect of distance on composition matrices for samples from rubber plantation (Mantel test, R = 0.22, P < 0.05). As the correlation was weak, most of the variation observed in community composition may be explicable by land use change with a small fraction as a result of pseudoreplication. According to Oksanen 65 , true treatment replication is sometimes impossible to achieve in a study, and thus its absence need not weaken the value of a study, provided that the results are judged properly.
Testing for phylogenetic signal using Mantel correlograms showed significant correlations between differences in OTU environmental optima and OTU phylogenetic distances but only across relatively short phylogenetic distances (ρ < 0.05, SI Appendix Fig. S4). We examined the relationship between βNTI and land-use type to infer changes in the relative influences of deterministic and stochastic assembly processes on the micro-and mesofaunal community. The pairwise comparisons of βNTI values within each land-use indicated similar patterns (Fig. 2). In forest and rubber plantation, the mean of pairwise βNTI values fell between the null distribution of βMNTD (|βNTI| <+2) indicating that community assembly is primarily influenced by stochastic processes.
Land-use and diversity. The rarefaction analysis showed that both richness and sample heterogeneity of OTUs were higher in the forest than in rubber plantation (Fig. 3). The α-diversity of 18S rRNA gene-based OTUs was significantly altered by forest conversion to rubber plantation. At a depth of 1885 reads, micro-and mesofaunal richness (OTUs) was significantly higher in the forest than in rubber plantation (GLMM; t = −2.28, p < 0.05). When measured separately for each of the three most abundant phyla (Table 1), only arthropods yielded significant results with OTU richness being higher in native forest (GLMM; t = −2.62, p < 0.01).  www.nature.com/scientificreports www.nature.com/scientificreports/ Other α-diversity indices also showed that diversity was significantly higher in native forests (glmm PQL results: Shannon; t = −4.7, p < 0.001, Faith's PD; t = −3.5, p < 0.01), compared to the rubber plantations. Similar results were found when tested separately for each of the three major phyla present except for nematodes which did not yield any significant difference between the two land use types (Table 1).
Multiplicative diversity partitioning indicated that the native forest had significantly higher (p < 0.001) metazoan αand β-diversity than rubber plantation at all levels in the hierarchy ( Table 2). Multiplicative partitioning yielded similar results for analyses done separately for each of the three most abundant phyla (Table 2). Annelids, when taken separately, showed exactly the same patterns of diversity as observed for the total metazoan community. Nematodes and arthropods, however, showed a different pattern, with higher diversity in rubber plantation at some (but not all) hierarchal scales. For arthropods, α-diversity when measured between forest sites and for forest overall emerged as significantly higher than for rubber plantation. By contrast, beta diversity when calculated between rubber plantation sites, and for the total rubber plantation dataset, was significantly higher than for forest samples. For nematodes, at all levels in the taxonomic hierarchy, the diversity was higher in samples from rubber plantation -except for total beta diversity.
Variation partitioning & RDA results revealed that in forest samples, environmental variables (i.e., pH, TOC & GSW) and Ele explained 25.3% of the total variation (ρ < 0.005) of metazoan communities, of which 21.9% of the total variation (ρ < 0.005) could be explained by soil variables alone ( Fig. 4; top). In isolation, spatial variables were insignificant and could not explain any of the variation observed (ρ > 0.5). Likewise, within the rubber plantation metazoan community, environmental variables (pH & ST) and Ele explained the largest proportion (9.9%) of the total variation (ρ < 0.01) (Fig. 4; bottom). Unlike the forest, around 3.3% of total metazoan community variation in the rubber plantation could be explained by spatial variables (ρ < 0.05), while abiotic and spatial variables together explained around 14.1% of the variation (ρ < 0.01). Similar results were found for each of the three most abundant phyla (SI Appendix Fig. S5).   Table 2. Multiplicative diversity partitioning for the Metazoan community; γ-diversity was multiplicatively partitioned into independent αand β-components following Jost 52 . (a) Separate runs for Native Forest and Rubber Plantation; (b) Single run for the overall dataset. *Alpha.1: between sites, Alpha.2: total, Beta1: between sites, Beta2: total. *Alpha.1: between sites, Alpha.2: between forest and rubber plantation, Alpha.3: within sites, Beta1: between sites, Beta2: between forest and rubber plantation, Beta.3: within sites. Significance level: ***p < 0.001, **p < 0.01, *p < 0.05. www.nature.com/scientificreports www.nature.com/scientificreports/ The MRM model explained a significant proportion of the variability in total phylogenetic community structure, with most of the environmental parameters in both land-use types ( Table 3). Among the environmental variables measured, soil pH and TOC were correlated to UniFrac distance and β-MNTD for the total phylogenetic community for the forest whereas only soil pH was correlated to β-MNTD for the samples from rubber plantation. Spatial distance was significantly correlated to the unweighted UniFrac for both the forest and rubber plantation samples, but to β-MNTD only for samples from rubber plantation (Table 3).
When checked for differences between the two land-use types for major soil variables, a significant difference in TOC (glmm PQL; t = −2.8, p < 0.05) and a marginally significant difference in TN (glmm PQL; t = −2.0, p = 0.051), was present. Soil C: N ratio was also lower in samples from rubber plantation (p < 0.05).

Discussion
Impacts of conversion to plantation on community composition. Conversion of the forest to rubber plantation results in consistent changes in soil micro-& mesofaunal community composition. In comparison to rubber plantation, a larger proportion of variation for forest can be explained using abiotic parameters (Fig. 4, Variation Partitioning Results). Similarly, MRM results (Table 3) suggest that soil pH and TOC delimit the community composition to a certain extent, mostly for forest sites. It seems that for forest samples, environmental variables are a major force in structuring communities. There may be some analogy with a recent soil microfaunal study from Antarctica where microfaunal distribution was controlled by nutrients such as pH, organic matter and soil moisture 66 .
The phylogenetic signal 55,67 , detected across relatively short phylogenetic distances, indicates that closely related micro-and mesofaunal OTUs tend to occupy similar niches (ρ < 0.05, SI Appendix Fig. S4). Similar patterns have been observed in other studies of both prokaryotic and eukaryotic microbial communities [62,68,69 . From phylogenetic null modeling 34 we found a strong impact of stochastic processes on the assembly of micro-and mesofaunal communities in soils of both forest and rubber plantation, suggesting that despite its other evident effects, the shift to a highly altered rubber plantation environment has not altered the role of stochasticity in community structure.
According to other studies, conversion from forest to rubber plantation is also accompanied by a progressive increase in soil bulk density and continuous decrease in pore space (soil compaction) in the 0-10 cm soil layer 70,71 , which may decrease the niche availability for many small soil animals. There are also significant changes in the other soil properties 70,71 , which may affect the microbial community which is, in turn, grazed upon by many soil metazoans. Physical disturbance in rubber plantation sites at XTBG also results from the prevalent management practices, which might suppress metazoan abundances. Reports have shown positive effects of limiting physical disturbance on metazoan abundance, especially the nematode fauna 72 .
Reduction in soil carbon on conversion to rubber plantation may also be an important factor. Soil carbon provides a resource for fungi and bacteria, on which many soil micro-and mesofaunal community feeds, and thus changes in this food source could adversely affect their diversity. Also possibly important in producing lower soil micro-and mesofaunal diversity is the dominance of a non-native tree species (rubber) which fewer soil microand mesofauna may be adapted to feeding off.
Considering the phyla individually, annelids were significantly more abundant in the rubber plantation (p < 0.05, Table 1). We did not find any indication of invasive earthworm species such as Pontoscolex corethrurus 73 , and most reads belonged to highly diverse small annelid families (Enchytraeidae or potworms, Naididae & Tubificidae). Arthropoda, which comprises the vast majority of the described tropical rainforest fauna and is estimated in total at around 2.5-3.7 million species 74,75 was represented by more OTUs in the forest (p < 0.01, Table 1). Soil calcium concentration is a limiting factor for terrestrial ostracods (Crustacea) 76 and millipedes (Diplopoda) 77 .
Nematoda showed no significant changes in relative abundance with land-use change, although the family level feeding guild structure of the nematode community showed striking differences in relation to the land-use type (SI Appendix Fig. S6). An increase in Bacteria Feeding (BF) along with Fungi Feeding (FF) guilds in rubber plantation could be because the Omnivore Predator (OP) guild are sensitive to land use change and need more time to establish compared to more rapidly growing BF and FF nematodes 78,79 . The Plant-Feeding (PF) guild was significantly reduced in the rubber plantation soils, which is unsurprising as their abundance has been found to be directly correlated to plant community richness and soil C: N ratio 80 .  Table 3. Multiple Regression on Matrices (MRM) analysis of metazoan communities using unweighted UniFrac and β-MNTD distance in forest and rubber plantation soils. The variation (R 2 ) is explained by the remaining variables, and the partial regression coefficients of the final model are reported. Significance level: ***p < 0.001, **p < 0.01, *p < 0.05, NS = non-significant.
www.nature.com/scientificreports www.nature.com/scientificreports/ Impacts of plantations on alpha-diversity. We predicted a substantial impact of forest conversion on soil micro-and mesofaunal α-diversity, and the results of this study suggest that indeed the conversion of the tropical forest of southern Yunnan to rubber plantation has resulted in a consistent decrease in α-diversity for soil micro-and mesofaunal community, except for nematodes, when these are considered separately. Possible reasons why α-diversity is lower in rubber plantations than in forest environments include: lack of small-scale heterogeneity in terms of various sizes/types of trees and shrubs, lower diversity of plants giving fewer organic substrates for niche differentiation in feeding, and a relative lack of carbon in litter and coarse debris flux resulting in fewer niches 81,82 .
Pesticides are generally used on rubber plantations, for example, to kill weeds, and are a possible factor producing the observed differences in soil metazoan diversity, as they are toxic to non-target species including soil animals 83,84 . A prevalent management practice in rubber plantations is repeated application of glyphosate, (alone or in a mixture with other pesticides) 31 which increases the presence of glyphosate residues in soils [85][86][87][88] and can be detected in top soils even after two years since the last spraying 86 . Glyphosate is reportedly a potent microbiocide 89,90 and has been reported to severely affect the casting activities and reproduction of earthworms 91 .
Another conventional control practice used for powdery mildew and anthracnose diseases of rubber plantations is annual spraying of sulfur powder during the leaf expansion period of the rubber tree. This application over a period of time has been reported to reduce soil pH 92 , and it could be one of the reasons behind the lower pH observed in our samples from rubber plantation. However, in a study Li et al. 92 sulfur spraying did not significantly alter the soil microbial community composition and function over the long term (48 years) some specific microorganisms, such as acidophilic bacteria, in the soil can use sulfur as an energy source and oxidize sulfur to SO 4 2− accompanied with a decrease in soil pH which will persist over the years 93 . Changes in soil physicochemical characteristics and microbial community with the use of pesticides -as shown by these two examples -may eventually have effects on the micro-and mesofaunal populations of the rubber plantation soils.
Ecosystem engineers, such as earthworms, have been found to show a drastic decrease in abundance under rubber monoculture (possibly due to decline in soil exchangeable Ca 94 and glyphosate herbicide 91 in soil from rubber plantations). This may influence the diversity of other soil animals by regulating the availability, diversity, and spatial distribution of resources available to them 95,96 , and may help to explain the lower α-diversity in samples from rubber plantation.
Impact of conversion to plantations on beta-diversity. As anticipated, and consistent with earlier studies of larger organisms, a decrease in micro-and mesofaunal β-diversity with the conversion from forest to rubber plantation was observed 97,98 except for arthropods and nematodes. Still, total β-diversity for nematodes in samples from rubber plantation was marginally lower than forest samples ( Table 2).
The uncoupling of αand β-diversity after ecosystem conversion that is observed here for arthropods is generally not seen in examples from plants and larger animals 98 : conversion usually results in decreases in both αand β-diversity in parallel as observed above for the total population. Although rarely found in higher organisms 99 the converse situation -a decrease in β-diversity along with an increase in α-diversity in response to human disturbances has been reported elsewhere 100,101 . Similar uncoupling between αand β-diversity has also been previously reported for bacterial communities, where the natural vegetation was either converted to agriculture 20 or plantation 21 .
Various ecological mechanisms could be behind the greater arthropod β-diversity we observed in rubber plantations. With a random loss of the species present in the original rainforest (associated here with deforestation and initial and continuous human disturbance) and limited dispersal, species richness at the α-diversity level will be low and few species will be shared between disturbed plots. This will lead to an increased β-diversity between rubber plantation plots. Apart from the effects of dispersal limitation 102,103 , the increased arthropod β-diversity in rubber plantations might be a consequence of some unknown aspect of habitat heterogeneity that leads to species sorting by environmental selection 103,104 .
The diversity of a plant community is likely to be one of the major environmental factors that influences the diversity of the soil micro-and mesofaunal community, through shaping the resource availability to the soil community 105 . Although soil micro-and mesofauna need very small patch sizes and comparatively smaller sources of energy and nutrients to maintain their diversity, soil animal diversity is expected to increase with greater plant diversity due to enhanced opportunities for niche differentiation with respect to habitat use and sources of energy and nutrients [106][107][108] . However, the patterns are not consistent, and it has been often shown that the effect of plant species identity is greater than that of plant species richness 107,[109][110][111] . The effect of plants on soil animal diversity is generally regulated in two different ways, i.e., availability of a resource and diversity of a resource 107,112 . Plants not only play a role in the availability of metabolic resources, but also have effects on other ecological factors that further shape soil habitats and niche space, including the presence and abundance of ecosystem engineers 95,96 , microclimatic conditions, and general soil properties.

Conclusion
It is clear that conversion of forest to rubber plantation involves a major turnover of lower level taxa (expressed as OTUs) of micro-and mesofauna, and an overall loss of diversity. Hypothetically, this loss of micro-and mesofaunal diversity could result in lower soil and ecosystem resilience, by analogy with other systems where removal of diversity has resulted in decreased stability 112 . It is also relevant to consider whether the loss of the original rainforest soil micro-and mesofauna might make it more difficult to re-establish tropical rainforest on former rubber plantations, as these organisms have important and specialized roles in nutrient cycling in tropical forests. Such questions should be investigated in further studies of this and other analogous systems.