Next Article in Journal
Bird Assemblage Recovery in a Chronosequence of Tropical Dry Forests in Costa Rica
Next Article in Special Issue
Evaluation of Soil Biodiversity in Alpine Habitats through eDNA Metabarcoding and Relationships with Environmental Features
Previous Article in Journal
Implications of Reduced Stand Density on Tree Growth and Drought Susceptibility: A Study of Three Species under Varying Climate
Previous Article in Special Issue
Responses of Soil Microbial Community Composition and Enzyme Activities to Land-Use Change in the Eastern Tibetan Plateau, China
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Effects of Livestock Pressure and Vegetation Cover on the Spatial and Temporal Structure of Soil Microarthropod Communities in Iberian Rangelands

by
Carlos Lozano Fondón
1,*,
Jesús Barrena González
2,
Manuel Pulido Fernández
2,
Sara Remelli
1,
Javier Lozano-Parra
3 and
Cristina Menta
1
1
Department of Chemistry, Life Sciences and Environmental Sustainability, University of Parma, Viale delle Scienze 11/A, 43124 Parma, Italy
2
GeoEnvironmental Research Group, University of Extremadura, Avda. Universidad s/n, 10071 Cáceres, Spain
3
Instituto de Geografía, Pontificia Universidad Católica de Chile, Avda. Vicuña Mackenna 4860, Santiago de Chile 7820436, Chile
*
Author to whom correspondence should be addressed.
Forests 2020, 11(6), 628; https://doi.org/10.3390/f11060628
Submission received: 27 March 2020 / Revised: 22 May 2020 / Accepted: 27 May 2020 / Published: 2 June 2020
(This article belongs to the Special Issue Relationship between Forest Biodiversity and Soil Functions)

Abstract

:
Forests, including their soils, play an important role since they represent a large reservoir of biodiversity. Current studies show that the diversity of soil fauna provides multiple ecosystem functions and services across biomes. However, anthropogenic practices often pose a threat to soil fauna because of changes in land use and soil mismanagement. In these terms, rangelands in the southwest of Spain present several problems of soil degradation related to livestock activity and soil erosion, the intensity of which compromises the soil fauna’s functions in the ecosystem. Therefore, the aim of this study is to evaluate the response of community metrics and the spatial distribution of soil microarthropods to livestock activity and vegetation in such ecosystems. A photo interpretation analysis of an experimental catchment used as a study area was developed to identify and classify the intensity of livestock pressure. A total of 150 soil samples were collected throughout 2018. Soil biological (CO2 efflux) and physical-chemical parameters (pH, bulk density, organic matter, and water contents), and such meteorological variables as precipitation, temperature, and evapotranspiration were considered as variables affecting the composition of microarthropod communities in terms of taxa diversity, abundances, and their adaptation to soil environment (evaluated by QBS-ar index). Results showed higher abundance of microarthropods and higher adaptation to soil environment outside the influence of trees rather than beneath tree canopies. Moreover, the classification of livestock pressure revealed by the photo interpretation analysis showed low correlations with community structure, as well as with the occurrence of well-adapted microarthropod groups that were found less frequently in areas with evidence of intense livestock activity. Furthermore, abundances and adaptations followed different spatial patterns. Due to future climate changes and increasing anthropogenic pressure, it is necessary to continue the study of soil fauna communities to determine their degree of sensitivity to such changes.

1. Introduction

Soils are one of the most important reservoirs of biodiversity in the world [1,2,3,4,5,6]. Today, problems such as climate change, an increasing of human population, changes in land use, and land abandonment compromise the diversity and functions of soil biota [7,8] and, subsequently, the provision of ecosystem services [2,4,8,9,10] by the soil complex. Moreover, soil functions and healthy soil communities are closely correlated, and, together, they are essential for safe and sustainable food production [9], and they also maintain ecosystem stability and resilience [8]. Furthermore, the diversity of the soil community is often used to provide soil quality indicators, such as the composition and abundance of microarthropod communities [11,12]. Indeed, it is widely accepted that soil microarthropods are very sensitive to disturbances because of their adaptation to a soil environment [12,13]. Although the contribution of microarthropods to the total amount of energy fluxes and biogeochemical transformations occurring in the soil is relatively low [14], they are a key component in enhancing the resilience and resistance of the soil food web by supporting structural stability [15] since they link microorganisms to macrofauna in the context of an interconnected network [16]. Such close correlations between bacterial and fungal channels to mesofauna [9] also determine top-down and bottom-up forces that modify the structure of the entire community and, therefore, the efflux of CO2 produced by the soil food web during its metabolic activity [3].
Many ecological functions have been attributed to soil microarthropod communities [17]. However, the functions they perform can be compromised via the reduction of biodiversity caused by disturbances in the soil environment [2]. Therefore, the loss of functional groups of microarthropods, such as detritivores, which are related to the soil carbon cycle could determine the interruption of several steps in the organic matter degradation chain [18,19]. In such a context, Mediterranean bioclimatic areas with semi-arid conditions, such as rangelands in the Southwestern Iberian Peninsula, are susceptible to this fact since they have been catalogued as ecosystems under risk due in future climate scenarios to land mismanagement and livestock intensification [20].
Dehesas are traditionally-managed rangelands commonly characterized by a two-layered vegetation structure: a savanna-like open tree layer (15–40 trees/ha) with an understorey pasture in the same land unit [21,22,23,24]. Moreover, pools of soil nutrients are frequently limited due to poor parent material and extremely arid conditions during the Mediterranean summer [21,24]. It is a system particularly subject to abandonment [22,25], soil degradation [22,24,26], and subsequent loss of soil biodiversity because of the increase of livestock density and the progressive abandonment of land by farmers. However, patches of vegetation are important for dehesas to maintain biodiversity associated with spatial heterogeneity [27]. In this context, trees play an important role in regulating environmental features such as soil temperature [28,29] and moisture [23]; the modification of chemical characteristics, such as availability of nutrients [30], and the direct promotion of the development of detritivorous microarthropod communities via the reduction of sunlight availability and litter inputs [31]. Moreover, such habitat heterogeneity at multiple spatial scales [22,30] could represent areas for the conservation of biodiversity in farmlands, as indicated by Moreno et al. [27] in which the authors used the term “habitat condition” to refer to areas that sustain certain levels of aboveground biodiversity in rangelands [32]. We adapted this concept to our study area in trying to define combinations of environmental features and elements of the landscape (mostly in reference to vegetation and livestock pressure) that could also drive the spatial distribution, structure, abundances, and adaptation to soil of microarthropod communities. In order to clearly define such combinations of factors in this work, we used the term “soil habitat condition” (SHC) [32].
Thus, the central questions of our study are: do different intensities of livestock activity induce changes in soil microarthropod communities? Is the structure of a microarthropod community affected by niche environmental factors associated with the presence of the tree? Do the adaptations and abundances of microarthropods follow different spatial patterns? With regard to these questions, we defined mainly three aims: (i) determining changes on microarthropod communities associated to seasonality, proximity to trees, and intensity livestock pressure; (ii) identifying the most sensitive biological forms of the microarthropod group to livestock pressure; and (iii) exploring the spatial patterns of microarthropod abundances and the occurrence of morphological traits that indicate high adaptation to the soil environment.

2. Study Area

Research was conducted on a farmland with agro-silvo-pastoral land use located in the province of Cáceres, in the southwest of Spain, where an experimental catchment was delimited (Figure 1). The study area (151.6 ha) is representative of a traditionally-managed system, commonly known as a dehesa, which is dominated by several vegetation layers including scattered oak trees (Quercus ilex L.), a shrub layer (Retama sphaerocarpa L.), and a herbaceous layer composed of annual species (grasses such as Vulpia bromoides L. [Gray], Bromus sp., Aira caryophyllea L., and legumes such as Ornithopus compressus L., Lathyrus angulatus L., and several species of Trifolium) [21,23]. Climate (Table A1 in Appendix A) is typical of the Mediterranean area, with semi-arid conditions characterized by cold winters and a period of hydric stress during the summer. Mean annual precipitation is about 524.2 mm. Rainfall events are common in autumn and spring; however, dry seasons and longer dry periods are frequent. Mean annual temperatures oscillates from 14° to 16 °C.
Geomorphologically, the study area is in old erosion surfaces (Figure 1A), which are formed by schist and greywacke of the Precambrian age [26]. Soils are shallow with a thickness of usually less than 50 cm [24,26]; soil textures are sandy-loam in low-slope areas and silty-loam in areas with a higher slope. Soils reactions oscillate from 4.3 to 7.3, and they are poor in organic matter (mean values are about 3% in the A horizon) [26]. They are classified as Luvisols and Cambisols [33].
Farm management is conventional: livestock walk freely inside the farm, which means that livestock charges per hectare inside the study area are not equally distributed. Moreover, the presence of several “points of reunion,” such as eating zones and water reservoirs, influence the frequency of trampling and grazing of surrounding areas close to them. In 2018, the livestock at the farm comprised 1200 sheep and goats (southeast area), 50 pigs (northwest area), 37 cows, and one bull (southwest and central areas of the farm).

3. Materials and Methods

3.1. Determination of the Intensity of Livestock Pressure

The study involved a description of the farm management by interpreting orthoimages (0.5 MP size) taken in 2016 by the Spanish National Information Center [34]. Parameters such as the density of the vegetation cover and the bare soil area were identified and related to livestock activity (mostly trampling and grazing) [30]. For the identification of zones with different grazing and trampling intensities, a supervised object-based image analysis (OBIA) classification [35,36] was used. The procedure was developed in the eCognition Developer 9 software (Trimble Germany Gmbh, Munich, Germany), avoiding the “salt and pepper effect” that occurs with pixel-oriented classifications [37,38].
Broadly, three categories were defined by OBIA based on the effects of livestock activity and the characteristics of the vegetation cover (Figure 1B). We then confirmed the field classifications: (1) SHClow: characterized by a shrub-encroached herbaceous layer, typically 40%–70% Retama sphaerocarpa L. cover with a dense tree layer, absence of bare soil and no signs of livestock pressure (i.e., defecation, trampling, or grazed vegetation); (2) SHCmedium: herbaceous layer, mostly 10%–40% of R. sphaerocarpa L. cover with a sparse tree layer, <10% of bare soil, and slight indicators of livestock presence, and (3) SHChigh: herbaceous layer with a sparse tree layer but no shrub cover, 50% or more bare soil, and evident signs of livestock pressure. See Figure 2 for an example of the general characteristics of each SHC on the field.

3.2. Soil Sampling

Two sampling campaigns were carried out in 2018. A total of 60 points were sampled in spring (April), and 90 in the late fall (December). In both campaigns, points were equally distributed among the three SHC categories previously identified by OBIA (Figure 1B). Inside each SHC, half of the points were established beneath the oak canopies (in-C) on the northern cardinal point in relation to the stem. The other half were located outside the canopy (out-C), at least 8 m far away from every tree (stem) [18]. Both in-C and out-C points were established considering the presence of herbaceous vegetation cover and avoiding bare soil because it is widely accepted that low densities of soil fauna occur in absence of vegetation [13,39]. We made this choice to accurately compare the different SHCs because bare ground was only present in SHChigh.
For the spring sampling, 60 PVC cylinders (21 cm ø and 15 cm height) were embedded at each point the day before the beginning of the sampling campaign, as recommended by the LI-COR 8100A protocol of soil CO2 efflux measurement [40]. The cylinders were installed up to 10 cm deep, leaving 5 cm of the total height free in order to carry out the measurement correctly. After the stabilization of the microbial soil community [40], sampling was conducted as follows at each point: three repeated measurements of soil CO2 efflux using a LI-COR 8100A survey device were executed; then three measurements of soil moisture were taken in three points beside the cylinder using a TDR device and, finally, three undisturbed soil cores were collected using a steel cylinder at a known volume (100 cm3) and a soil sample extractor. The soil volumes collected from inside the cylinders used for the soil CO2 efflux measurement (approx. 3.5 dm3) were taken to the laboratory, where microarthropod extraction was carried out. Once organisms were collected, soil samples were meshed at 2 mm ø. Then, six replicates were picked up to determine the pH and soil organic matter (SOM) content of each sample (3 for pH and 3 for SOM). Due to logistical reasons, measurements of soil CO2 efflux, pH, water content, and bulk density could not be revealed in the autumn campaign. At this time, only SOM has been measured in the 90 points sampled.
In the laboratory, pH was determined by dissolving 1 g of soil in 3 g of H2O using a pH meter. SOM was revealed using the loss-on-ignition method. Soil cores inside the steel cylinders were used to calculate soil bulk density and gravimetric moisture by the wet-minus-dry weights of the samples in relation to the volume of the undisturbed soil cores (100 cm3).

3.3. Analysis of the Microarthropod Communities

The analysis of the microarthropod community was based on the QBS-ar methodology [11,12]. The QBS-ar index (i.e., biological soil quality based on arthropods) evaluates the capacity of a soil to harbor animals that are sensitive to disturbances because of their morphological characteristics. Therefore, based on the number of well-adapted microarthropods to the soil environment at a given time, it is possible to make a judgment about the quality of the soil in a given area (i.e., the higher the number of such organisms, the higher the soil quality). The QBS-ar of a soil sample is calculated as the sum of the ecomorphological indices (EMIs) of each biological form. The EMI is a dimensionless score that varies between 1 and 20, and it evaluates the degree of adaptation of the morphological traits that soil animals share by evolutionary convergence. For more details on QBS-ar application, see [13,17].
In this study, soil microarthropods were extracted from the 150 samples (both spring and autumn) using Berlese-Tullgreen funnels (2 mm mesh size) and conserved in 70% ethanol solutions. The extraction time was about eight days, depending on the humidity of samples. Then, the extracted microarthropods were observed under stereomicroscope (40×), counted, identified, and classified as indicated in Table A2 from Appendix A. Once analysis of the microarthropod community was completed and raw data were collected, taxa diversity was evaluated via the Shannon’s index; a value of QBS-ar was associated with each soil sample, and the numbers of individuals per taxon for each soil sample were obtained.

3.4. Statistical Analyses

Three response variables were considered for this work: abundances, taxa diversity (defined by Shannon’s index), and community adaptation to soil environment (defined by the QBS-ar index). A first screening of the data was carried out following the protocol proposed by Zuur et al. [41]. Three categorical predictors were considered: (i) the season in which samplings were accomplished; (ii) the location of each sampling point (in-C: beneath tree canopy, out-C: outside the canopy); and (iii) the SHC representing the intensity of livestock activity surrounding each sampling point.
Initially, a three-way ANOVA test was performed to evince statistical differences among categories. Abundances were log-transformed based on the range of the data and significance level was established at p = 0.05. A Tukey pairwise test was applied post hoc to highlight the significant differences between pairs of categories. Once seasonal variation of community metrics was statistically confirmed, the rest of the analyses were carried considering seasons separately.
Prior to analysis, collinearity was tested with Pearson’s correlation coefficient in order to eliminate variables with identical trends [41]. When Pearson’s correlation was found to be higher than 0.4999, covariates were considered as collinear and, subsequently, one was excluded from the analysis [42]. Methods such as, non-metric multidimensional scaling (NMDS) and non-parametric permutational multivariate ANOVA (PERMANOVA) were chosen to study dissimilarities in microarthropod communities. NMDS based on Bray–Curtis distances was used to order the relationships among communities’ compositions in a specified number of axes [42]. A stress level score of ≤0.2 was used to account for goodness of fit. PERMANOVA, also based on Bray–Curtis distances, was then used to study environmental variables causing dissimilarity in the community structure [43,44]. In order to identify the sensitivity of each biological form, NMDS was also applied on each taxon based on the Bray-Curtis dissimilarity index. In order to accomplish these aims, community matrices were split into (1) two log-transformed abundance matrices (60 soil samples (rows) × 27 taxa (columns) in spring; 90 soil samples (rows) × 27 taxa (columns) in autumn); and (2) two EMI value matrices (60 soil samples (rows) × 27 taxa (columns) in spring, and 90 soil samples (rows) × 27 taxa (columns) in autumn) representing the morphological adaptation of biological forms to the soil environment. Environmental factors were summarized in a matrix presenting soil parameters, categorical predictors (SHC and out-C/in-C locations), geospatial characteristics (UTM coordinates, slope, and altitude of each sampling point), and meteorological variables such as maximum, minimum and mean temperature of the sampling day; average of maximums, minimums and mean temperatures of the 20 days prior to the sampling day, effective precipitation of the sampling day, effective cumulative precipitation of the 20 days prior to the sampling day, evapotranspiration of the sampling day, average evapotranspiration of the 20 days prior to the sampling day, and average hydrological balance of the 20 days prior to the sampling day (22 columns in spring and 18 in autumn). A stepwise model selection based on the significance criterion was used to choose the best combination of variables explaining the variance of the data. These analyses were carried out with the “vegan” package [45] from RStudio.
Generalized additive models (GAMs) were applied to investigate the effects of spatial distribution of niche-environmental factors upon spatial distribution of QBS-ar and total log-transformed abundances. In order to model the dispersion of community metrics across the space, two protocols were executed to run the models: the first, a random effect on “pure” spatial coordinates was used in order to seek spatial dependence of the response variables; on the second, a random effect on ordination coordinates extracted from NMDS replaced the spatial coordinates in order to find the best descriptor of the community metrics variation [41]. Stepwise model selection was based on Akaike’s information criterion (AIC) [46]. The R package “mgcv” [47] was used to perform this analysis.
Graphics were generated using the “ggplot2” [48] package from RStudio.

4. Results

4.1. Soil Parameters

Organic matter and water contents, bulk density and pH were measured in the three SHC in both in out-C and in-C locations (Table 1). Averages slightly differed based on SHC categorization. However, decreasing values of organic matter content were found from SHChigh to SHClow in out-C locations; differences were less evident when pH, bulk density, and water content were compared. Soil pH was found to be acidic (5.61 to 5.96) inside the study area with no broad variations both in either out-C or in-C. Moreover, bulk density averages in out-C varied slightly around 1.5 g cm−3. However, values of the soil parameters on in-C were smaller in the case of bulk density (≈1.2 g cm−3), but far higher in organic matter (≈10% compared to ≈5% for in-C and out-C, respectively). Instead, water content and soil CO2 efflux were the less-variable parameters considering both out-C and in-T locations (≈20% and ≈5 µmol m−2 s−1, respectively).

4.2. Seasonallity, Environmental Characteristics, and Community Structure

Generally, 113,579 organisms belonging to 27 taxa were individually identified and counted. Collembola, Acari, larvae of dipterans, and larvae of coleopterans were the most frequent taxa (58%, 35%, 3%, and 1%, respectively) representing 97% of the total abundances. In detail, frequencies varied significantly from spring to autumn (p < 0.01), as evidenced by results of the three-way ANOVA showed in Table 2. In spring, collembolans represented 50% of total abundances, which increased to 62% in autumn. Mite populations were higher in spring (44%) than in autumn (32%), and larvae of dipterans and larvae of coleopterans maintained their populations varying from 2% to 3% in the case of dipterans, remaining at 1% in the case of coleopterans. Other taxa, such as coleopteran adults, hemipterans, pauropods, thysanopterans, and ants, were frequently found but their abundances were lesser. Moreover, diplopods were found only in spring, and isopods conversely only in autumn, as well as individuals belonging to Zygentoma taxon. However, the most important source of variation was the location of each sampling point (in-C or out-C), showing a strong influence on the response variables. Although SHC was not revealed as a significant source of variation by the ANOVA, some taxa, such as diplopods, were found only in SHClow when trees were present. Embiopterans were found in the three SHCs, both in out-C and in-C, during spring campaign, but only in in-C locations during autumn. Proturans and pseudoscorpions followed similar patterns, as they were found in the same categories, but pseudoscorpions lacked in autumn in-C-SHClow. Such variations were also reflected by the QBS-ar and H′ indices, the values of which differed not only according to seasonality, but also to the location of sampling points (F = 4.490 and F = 6.232, both significant, respectively) as showed in Figure 3 and Table 2. The highest mean value of QBS-ar was detected in out-C-SHCmedium, followed by out-C-SHChigh, both in autumn. Generally, higher values of QBS-ar were found in out-C. Instead, mean values of H′ were closer to 1 in in-C locations in spring, but not in autumn.
See details in Appendix A: results of the post hoc Tukey test are shown in Table A3, mean values of QBS-ar and H’ are in Table A4 and Table A5, and total abundances are in Table A6 and Table A7.
Results of NMDS and PERMANOVA are reported in Figure 4 and Table 3, respectively. The ordination of sampling points based on community abundances showed stress values below the 0.15 threshold (Figure 4A,B). However, the occurrence of high EMI scores did not show convergence, as evidenced by stress values above 0.20 (Figure 4C,D). The scaling of community abundances was characterized by high overlapping based on categorical predictors (locations and SHC). Nevertheless, more dissimilarity among communities was attributed to locations (p < 0.05 and p < 0.001 in spring and autumn, respectively) than to SHC, which was significant in autumn (p < 0.05) but not significant at all in spring (Table 3). Representation of significant taxa abundances in spring was lower than in autumn, as shown in Figure 4A,B. In spring, such groups as Acari and Coleoptera (adults and larvae) avoided points where values of organic matter and pH were higher, but they were also positively related to high bulk density. Otherwise, Collembola, Thysanoptera, larvae of Diptera, and larvae of Coleoptera were positively related to points where soil water content was higher. In autumn (Figure 4B), Collembola, Chilopoda, and Pauropoda avoided points where SOM was higher. Moreover, these groups were significantly related to out-C-SHCmedium and out-C-SHClow in autumn. Several environmental variables were significant in both seasons, but after the model selection performed on PERMANOVA, only water content (p < 0.01), pH (p < 0.05), and soil CO2 efflux (p < 0.05) were deemed significant (R2 = 17.4%). In the autumn model, environmental variables, such as slope (p < 0.05), SOM (p < 0.01), effective precipitation (p < 0.05), and mean temperatures of the 20 days prior to sampling day (p < 0.05), explained a wider percent of the total variance of PERMANOVA (R2 = 22.1%) as compared to the spring model. It is noteworthy that the location of sampling points was found significant in both abundance models (p < 0.05 and p < 0.001 in spring and autumn, respectively).
As the NMDS ordinations based on EMI matrices were almost random in both seasons, a clear effect of SHC categories causing differences among sampling points was not determined. Nonetheless, segregation of communities by locations was quite evident in autumn as indicated by PERMANOVA (p < 0.001) in Table 3 and Figure 4D. Moreover, the number of significant taxa that fit the scaling based on EMI scores was higher when compared to taxa abundance ordination in both seasons. Once again, variables fitting communities was higher on NMDS than in PERMANOVA, but after the model selection, soil water content was the only significant variable related to spring communities (p < 0.05). In contrast, the autumn model was related to the variation of slope (p < 0.01) and mean temperature of the 20 days prior to sampling day (p < 0.01), including the location factor. Despite this, and similarly to abundances, total variances explained by EMI spring and autumn models were low (R2 = 8.4 and R2 = 19.7%, respectively).

4.3. Spatial and Temporal Patterns of Abundances and QBS-ar

GAMs demonstrated how both response variables (total abundances and QBS-ar) changed in spring and autumn in relation to the spatial structure of environmental variables (Table 4 and Figure 5). In general, variance explained by the models was very high (R2 = 87.4% and R2 = 89.1% for abundances; R2 = 85.5% and R2 = 91.8% for QBS-ar in spring and autumn, respectively), the largest effect of which was the contribution of the random effect in NDMS coordinates (F = 21.750 and F = 23.880 for abundances in spring and autumn; F = 12.870 and F = 37.430 for QBS-ar in spring and autumn, respectively). Moreover, the smoothness of the models was conditioned by location factor, where out-C significantly explained a wider proportion of the variance over the unexplained (F = 7.920 and F = 8.923 in spring and autumn, respectively) rather than in-C. Also, effective precipitation of 20 days prior to sampling day resulted statistically significant in autumn (F = 2.539). Otherwise in autumn, QBS-ar scores followed similar spatial patterns to SOM content in autumn (F = -2.214). It is important to note that smoothing was better when several variables such as pH, SOM and bulk density (spring abundances), temperature (autumn abundances), and effective precipitation (autumn QBS-ar) were available. Despite their insignificant effects on smoothing, AIC obtained better scores when they were present.
Figure 5A,B shows the spatial patterns of smooth isolines for total abundances. The maximum order of magnitude in autumn (e4) indicates a higher variation than in spring (e2). Instead, smooth isolines for QBS-ar reached variations from 0 to 80 in both seasons (Figure 5C,D).

5. Discussion

Understanding the ecosystem processes governing reservoirs of soil biodiversity, and the practices threatening them (e.g., anthropogenic mismanagement), would strongly benefit from characterizing the microarthropod community composition associated with traditionally managed rangelands. Moreover, the use of morphological traits in identifying spatial patterns and diversity of biological forms and relationships with above- and belowground environmental characteristics is a helpful method to detect areas under risk in terms of loss of soil multifunctionality [4,12]. This need becomes even more urgent since Mediterranean areas are especially sensitive to new climate change scenarios [7]. Features of Iberian rangeland ecosystems, such as the patchy distribution of vegetation and the unequal pressure exerted by livestock, may be a major structuring force of soil microarthropod communities at local scales. Nevertheless, several stochastic events (e.g., colonization and extinction) usually more related to regional scales can also take place locally [49].

5.1. The Response of Community Abundances

Our analyses suggest that abundances of microarthropod communities differed in response to the presence of trees and livestock pressure. These differences among community and taxa abundances can be simply explained by the environmental characteristics of each sampling point [50]. In this context, the analysis based on NMDS allowed us to identify dissimilarities among communities’ compositions, and sensitivity of taxa populations in relation to spatial and temporal dynamics of niche-environmental parameters. Our results suggest that, in general, areas outside the influence of trees harbor higher abundance of microarthropods rather than areas beneath tree canopies. NMDS also reveals a great degree of overlap between locations × SHC categories that could be representative of a geographic dispersal of taxa (at the study area scale) due to the absence of physical barriers and proximity of categories [51]. It is noteworthy that approximately 80% of variation was undetermined by the multivariate analysis (PERMANOVA), thus, it could be due to other non-spatially structured environmental factors that were not measured in the field [52], or even because of biotic interactions within the microarthropod community due to its spatial aggregation [49,53]. Results of GAM, based on abundances in both seasons, support the results of PERMANOVA. GAM models reflected a clear spatial structure of total abundances, as indicated by approximately 88% of variance explained in both seasons. This indicates that spatiotemporal processes (e.g., dispersal) in relation to local environmental factors drive total abundances of microarthropod communities [52,54,55,56]. Therefore, soil environmental variables, characteristics of vegetation cover, and livestock pressure alone could not explain by themselves the community aggregation phenomena but a dispersion of total abundances. This fact coincides with [53,54], who concluded that the spatial structure of variables plus the spatial and temporal structure of total abundances reveals that overall abundances are mediated both by dispersal and environmental factors, where the effects of the latter are weaker.
Several relationships inferred from NMDS analysis are relevant. Significant taxa, such as Acari, Collembola, Pauropoda, Araneae, and Chilopoda, were related to out-C-SHCmedium and out-C-SHClow in autumn, while Acari, Collembola, and Thysanoptera were positively associated with out-C-SHCmedium and out-C-SHClow in spring. Such convergence of positive correlations in the same direction could indicate that lower livestock pressure allows the development of populations and the co-occurrence of microarthropod taxa in the absence of trees, which accords with Mulder et al. [57]. However, an ensemble of biotic and abiotic interactions can also cause such relationships: for instance, soil parameters analyses indicate that trees contribute to rising SOM values, mostly due to litter inputs [21], which provide habitat and energy budget to the detrital community. In addition, significant negative relationships between pH, SOM, and abundances of such groups as Acari and Coleoptera in spring, as well as between SOM, and Collembola and Chilopoda (and Coleoptera) in autumn, have been identified. Except for Chilopoda (mainly predators), the other taxa show a wide variety of feeding habits. Indeed, the spectrum of feeding habits within the microarthropod community is wide [5]. However, our study is limited since the assemblage of communities was performed at a low taxonomical level. Thus, in order to shed light on relationships between feeding habits, environmental parameters, and community structure, an analysis of the functional roles of microarthropods at higher taxonomical level should be performed.
Another inference suggested by results based on NMDS analysis is that Acari and Collembola abundances avoid each other, as mites were more related to out-C, whereas collembolans were more related to in-C areas. However, this affirmation must be put in context: mites and collembolans are dominant and normally coexist in soil. Locally, competition for resources between both taxa (e.g., those with detritivorous feeding habits) is expected to benefit some functional groups over others due to negative interactions, as suggested by Caruso et al. [54] in a study of Antarctic microarthropod communities. As our study lacked a functional characterization of taxa, this affirmation is not sustained by our results. Otherwise, several authors [24,29] indicate that trees (and litter) in semiarid wood pastures decrease soil evapotranspiration rates by sunlight interception [31], and consequently decrease wider fluctuations of soil moisture when compared to open spaces. This causes a response of the entire microarthropod community to light availability as demonstrated by Jiménez-Chacón et al. [31], who concluded that detritivores preferred darker microsites (e.g., beneath tree canopies). Hence, a higher accumulation of litter and SOM beneath tree canopies, associated with lower rates of evapotranspiration, promoted collembolan abundances and likely inhibited the development of detritivorous mites by competition for resources. The opposite occurs in localities where mite populations are greater than collembolans, but our results do not allow for speculation about such dynamics.

5.2. The Response of Biological Forms’ Evolutive Adaptation

In the context of this study, rangeland mismanagement leads to a promotion of undesirable vegetation in terms of livestock profit. However, from the perspective of biodiversity, a greater number of well-adapted microarthropod communities live in such areas. The measurement of biological forms’ adaptation to soil environment involves the evaluation of such traits as depigmentation; reduction of appendages, such as antennae and legs; presence, absence, or reduction of the visual apparatus; presence, absence, or reduction of wings; dimensions; and body shape [11]. The higher the score attributed to each trait, the better adapted to a soil environment the organism is. The sums of all these scores serves as the EMI of each biological form. In the same way, the greater the number of biological forms with high scores that lives in the soil, the better quality and stability of the soil (i.e., the less disturbance) [58]. This is the main concept upon which QBS-ar relies and, by definition, it is based on environmental filtering theory. Overall, our analyses suggest that microarthropod communities’ evolutive adaptation to their soil environment differed mostly in response to the presence of trees. These differences reinforce the hypothesis that vegetation cover and environmental characteristics (i.e., habitat degradation caused by livestock pressure) are major forces that structure microarthropod communities even at evolutive level, which accords with the main basis of QBS-ar and environmental filtering theory.
Stress values of the NMDS analysis based on EMI communities suggest that the goodness of fit on morphological traits was poor. Nevertheless, a clear and significant differentiation between out-C and in-C was identified (stronger in autumn than in spring). This supports our hypothesis that vegetation cover shapes the microarthropod adaptation, but it contradicts the result obtained by Meloni et al. [39], who studied the community composition of ground arthropods, in terms of abundances and richness, in patches of vegetation and interpatch in Mediterranean drylands. In that study, patches of vegetation harbored higher richness and abundances than in interpatch areas when compared to the open space values of our study. However, in Meloni et al.’s study it is noteworthy that vegetation cover on interpatch areas was absent (bare soil).
NMDS analysis also showed a great degree of overlap between SHC categories, and insignificant effects on EMI communities. Overall, a greater number of taxa showed significant fits to SHC areas and locations. Unfortunately, it is difficult to extract clear inferences in relation to SHC based on these results due to stress values over 0.20. However, disturbances driven by livestock could explain it since several authors consider that well-managed silvopastoral systems, for instance with livestock charges at 1 AU ha−1 or below [30], could enhance resource allocation within soil food webs [59], by, for example, altering the C:N ratio [60]. This fact supports the niche-environmental hypothesis, and it could explain why SHCmedium showed similar patterns to SHClow on QBS-ar values.
Finally, poor values of variance explanation resulted from PERMANOVA. Approximately 87% of the total variance remained undetermined. This implies that the occurrence of morphological adaptation may be related to spatially structured variables or biotic interactions (or both) that were not measured in the field. Results of GAM analyses based on QBS-ar confirm the hypothesis that morphological adaptation also follows a spatially structured distribution. Moreover, it is even stronger than total abundances models. Smoothing patterns also differed from spring to autumn, which were negatively related to the spatial position of trees and SOM in autumn, and only to trees in spring. Therefore, the response to the third question of this work is that morphological adaptation and abundances did not follow identical, but similar, spatial patterns as confirmed by NMDS, PERMANOVA and, finally, GAM analysis.

5.3. Object-Based Image Analysis and SHC Classification

Correlations between microarthropod communities’ structures, metrics and evolutive adaptations and SHC classification using OBIA were not as high as expected. This might have been due to the fine scale at which microarthropod populations develop themselves, or to the relatively rapid dynamics of annual grasses. OBIA was chosen as the best candidate to remotely classify objects on the ground when compared to pure pixel classification techniques, but it obviously presents problems regarding the pixel’s dimensions of the image. However, OBIA turned out to be a useful technique to identify livestock effects, which would be a useful analysis when performed at larger scales. That being said, we were able to confirm that the results of the analysis corresponded with the areas in which livestock spend more time, as demonstrated in precedent studies about physical-chemical indicators of soil quality [25], impacts of livestock [30], and soil erosion studies [30,61] realized within the study area.

6. Conclusions

Results of this study suggest that there is a clear effect of spatial heterogeneity and spatially distributed variables (measured and unmeasured in the field) on structuring community metrics and community composition. This study demonstrates several facts: (1) landscape characteristics play a crucial role on the occurrence of evolutive adaptation of microarthropod biological forms; (2) abundances and the occurrence of morphological adaptation did not follow identical, but similar spatial patterns; (3) and the effect of environmental characteristics, such as the patchy distribution of vegetation being high on abundances and taxa diversity, which is likely due to environmental filtering, but to stochastic dispersal as well. In contrast, environmental filtering better explains the spatial distribution of QBS-ar and community composition based on EMI scores. (4) Higher abundances and adaptation to soil environment were related to open spaces rather than areas under arboreal influence. Smoothing of GAM models responded to the spatial positions of trees in terms of overall abundances and QBS-ar in both seasons. Moreover, the contribution to GAM models of the spatial structure of soil parameters and livestock pressure was unexpectedly low in abundances, and almost absent in QBS-ar patterns (with the exception of SOM content in autumn). This indicates that stochastic dispersal in relation to local environmental factors (e.g., non-spatially structured abiotic factors, as well as biotic interactions) drive abundances and adaptation of microarthropod communities. (5) High livestock pressure influenced microarthropod communities’ composition and metrics. Better values of community metrics were reached in medium and low areas, which indicates that lower livestock activity trends to enhance microarthropod zoocenoses. (6) The delimitation of SHC areas via OBIA technique showed unexpectedly lower correlations with microarthropod communities’ composition. For future studies, we strongly recommend the use of UAVs and multispectral images in order to reduce pixel dimensions, but also to determine the physiological state of vegetation. In this context, identifying the links between vegetation and belowground communities might be crucial to accurately quantify the resilience of ecosystems, and the consequences of climate change for humankind.

Author Contributions

Conceptualization, development of the idea, as well as experimental design and sampling: C.L.F., J.B.G., and M.P.F.; laboratory activities: C.L.F. and J.B.G.; data provision for completion of the database: J.L.-P.; statistical analyses: C.L.F. and S.R.; writing—original draft preparation: C.L.F.; supervision of the whole work: C.M.; review and editing: C.M. and J.L.-P. All authors have read and agreed to the published version of the manuscript.

Funding

This research received no external funding.

Acknowledgments

We are very grateful to the GeoEnvironmental Research Group (GIGA) from the University of Extremadura who provide us with technical support, instruments, and facilities to carry out our analysis. Moreover, we are also grateful to undergraduate students who help us in laboratory works.

Conflicts of Interest

The authors declare no conflict of interest.

Appendix A

Table A1. Means ± SD of meteorological parameters calculated for the period 01/01/2013 to 31/12/2018. From: Redarex; Meteorological station: Valdesalor (CC18) which, is approximately 31 km away from the study area; altitude: 382 m.; coordinates UTM H30 X: 730,101, Y: 4,361,000; Extracted from [62].
Table A1. Means ± SD of meteorological parameters calculated for the period 01/01/2013 to 31/12/2018. From: Redarex; Meteorological station: Valdesalor (CC18) which, is approximately 31 km away from the study area; altitude: 382 m.; coordinates UTM H30 X: 730,101, Y: 4,361,000; Extracted from [62].
Meteorological VariableUnitsValue
Annual solar radiationW/m216.59±6.76
Net solar radiationW/m27.52±4.22
Mean annual temperature°C15.02±6.30
Maximum mean temperature of the coldest month°C15.84±0.62
Minimum mean temperature of the coldest month°C2.41±0.27
Maximum mean temperature of the warmest month°C30.60±6.46
Minimum mean temperature of the warmest month°C13.61±4.22
Mean annual rainfallmm524.2±28.4
Mean annual effective precipitationmm249.8±14.4
Mean annual evapotranspirationmm1363.1±75.5
Table A2. Eco-morphological Indexes (EMIs) score for each microarthropod taxa. Groups shown are those found in this work.
Table A2. Eco-morphological Indexes (EMIs) score for each microarthropod taxa. Groups shown are those found in this work.
TaxaEMI Score
Pseudoscorpiones20
Opiliones10
Araneae1–5
Acari20
Isopoda10
Diplopoda10–20
Pauropoda20
Symphyla20
Chilopoda10–20
Protura20
Diplura20
Collembola1–20
Psocoptera1
Hemiptera1
Thysanoptera1
Zigentomi10
Embioptera10
Orthroptera1–20
Coleoptera1–20
Hymenoptera1–5
Diptera1
Lepidoptera1
Coleoptera (larvae)10
Diptera (larvae)10
Hymenoptera (larvae)10
Lepidoptera (larvae)10
Holometabolans (adults)1
Table A3. Significant comparisons from post hoc Tukey tests are shown.
Table A3. Significant comparisons from post hoc Tukey tests are shown.
MetricsFactorsPairs ComparisonDifferencep
Ln abundancesLivestock pressureLow − High0.5680.045
LocationIn-C − Out-C–0.07<0.001
SeasonSpring − Autumn–0.0510.007
Livestock pressure × LocationHigh × In-C − Low × Out-C–1.1080.009
High × In-C − Medium × Out-C–1.1470.004
Medium × In-C − Medium × Out-C–0.90.03
Livestock pressure × SeasonLow × Autumn − High × Autumn0.8820.038
High × Spring − Low × Autumn–1.0930.032
Medium × Spring − Low × Autumn–0.9240.036
Location × SeasonIn-C × Autumn − Out-C × Autumn–0.6510.036
In-C × Spring − Out-C × Autumn–1.199<0.001
Livestock pressure × Location × SeasonHigh × In-C × Spring − Low × Out-C × Autumn–1.5260.048
High × In-C × Autumn − Medium × Out-C × Autumn–1.4270.043
High × In-C × Spring − Medium × Out-C × Autumn–1.7330.019
High × In-C × Spring − Low × In-C × Autumn–1.7810.035
QBS-arLocationIn-C − Out-C–19.693<0.001
Livestock pressure × LocationHigh × In-C − High × Out-C–24.2220.037
High × In-C − Low × Out-C–28.1010.009
High × In-C − Medium × Out-C–33.791<0.001
Medium × In-C − Medium × Out-C–23.9220.027
Livestock pressure × SeasonHigh × Spring − Low × Autumn–23.2950.09
Location × SeasonIn-C × Autumn − Out-C × Autumn–27.661<0.001
In-C × Spring − Out-C × Autumn–24.2430.002
Livestock pressure × Location × SeasonHigh × In-C × Autumn − High × Out-C × Autumn–39.4190.08
High × In-C × Autumn − Low × Out-C × Autumn–34.0040.051
H′Location × SeasonIn-C × Spring − Out-C × Spring0.1970.063
Table A4. Average QBS-ar values ± SD found at each SHC (high, medium and low livestock pressure), in each Location (beneath = in-C and outside the canopy = out-C) during both sampling campaigns.
Table A4. Average QBS-ar values ± SD found at each SHC (high, medium and low livestock pressure), in each Location (beneath = in-C and outside the canopy = out-C) during both sampling campaigns.
SpringAutumn
Livestock PressureOut-CIn-COut-CIn-C
SHChigh66.9 ± 28.372.1 ± 27.2104.4 ± 34.264.9 ± 23.2
SHCmedium97.1 ± 39.480.3 ± 22.2105.9 ± 27.876.8 ± 32.3
SHClow91.3 ± 22.484.8 ± 22.798.9 ± 18.489.0 ± 31.5
Table A5. Average Shannon’s diversity values ± SD found at each SHC (high, medium and low livestock pressure), in each Location (beneath = in-C and outside the canopy = out-C) during both sampling campaigns.
Table A5. Average Shannon’s diversity values ± SD found at each SHC (high, medium and low livestock pressure), in each Location (beneath = in-C and outside the canopy = out-C) during both sampling campaigns.
SpringAutumn
Livestock PressureOut-CIn-COut-CIn-C
SHChigh0.79 ± 0.360.78 ± 0.330.97 ± 0.270.89 ± 0.23
SHCmedium0.76 ± 0.201.04 ± 0.410.86 ± 0.340.83 ± 0.29
SHClow0.69 ± 0.291.00 ± 0.480.76 ± 0.260.68 ± 0.24
Table A6. Absolute numbers of microarthropod found at each SHC (high, medium and low livestock pressure), in each Location (beneath tree canopy = In-C or outside the canopy = Out-C) during the first sampling campaign.
Table A6. Absolute numbers of microarthropod found at each SHC (high, medium and low livestock pressure), in each Location (beneath tree canopy = In-C or outside the canopy = Out-C) during the first sampling campaign.
Spring
Out-CIn-C
TaxaSHChighSHCmediumSHClowSHChighSHCmediumSHClow
Pseudoscorpiones-3---3
Opiliones------
Araneae1-3-22
Acari34103928288290211162030
Isopoda------
Diplopoda-----2
Pauropoda1813--
Symphyla-103-3-
Chilopoda217184198
Protura-1---3
Diplura30292437-
Collembola9954368104634245584936
Psocoptera---527
Hemiptera15130143413
Thysanoptera12476208
Zigentomi------
Embioptera12-642
Orthroptera---1--
Coleoptera4618028251812
Hymenoptera526641715331
Diptera-12276
Lepidoptera--111-
Coleoptera (larvae)363130143134
Diptera (larvae)292186364120135
Hymenoptera (larvae)---4--
Lepidoptera (larvae)131-22
Holometabolans------
Total455989004203449571233234
Table A7. Absolute numbers of microarthropods found at each SHC (high, medium and low livestock pressure), in each Location (beneath tree canopy = In-C or outside the canopy = Out-C) during the second sampling campaign.
Table A7. Absolute numbers of microarthropods found at each SHC (high, medium and low livestock pressure), in each Location (beneath tree canopy = In-C or outside the canopy = Out-C) during the second sampling campaign.
Autumn
Out-CIn-C
TaxaSHChighSHCmediumSHClowSHChighSHCmediumSHClow
Pseudoscorpiones161-1-
Opiliones------
Araneae1512431
Acari4003537510370132125992498
Isopoda----1-
Diplopoda------
Pauropoda163835612997
Symphyla18920712
Chilopoda5112-111
Protura227-11
Diplura81--31
Collembola70799787134724708489210,574
Psocoptera-----2
Hemiptera10222182431217
Thysanoptera34291173613
Zigentomi--31--
Embioptera---725
Orthroptera------
Coleoptera968529103632
Hymenoptera31717805810
Diptera31151672513
Lepidoptera------
Coleoptera (larvae)152101921375269
Diptera (larvae)13747122262388471
Hymenoptera (larvae)------
Lepidoptera (larvae)7016212142
Holometabolans361121
Total12,42416,38224,3106700783913,410

References

  1. Jeffery, S.; Gardi, C.; Jones, A.; Montanarella, L.; Marmo, L.; Miko, L.; Ritz, K.; Peres, G.; Römbke, J.; Van der Putten, W.H. European Atlas of Soil Biodiversity; Publications Office of the European Union: Luxembourg, 2010. [Google Scholar]
  2. Adhikari, K.; Hartemink, A.E. Linking soils to ecosystem services—A global review. Geoderma 2016, 262, 101–111. [Google Scholar] [CrossRef]
  3. Lal, R. Soil carbon sequestration to mitigate climate change. Geoderma 2004, 123, 1–22. [Google Scholar] [CrossRef]
  4. Wagg, C.; Bender, S.F.; Widmer, F.; van der Heijden, M.G.A. Soil biodiversity and soil community composition determine ecosystem multifunctionality. Proc. Natl. Acad. Sci. USA 2014, 111, 5266–5270. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  5. Coleman, D.C.; Wall, D.H. Chapter 5: Soil Fauna: Occurrence, Biodiversity, and Roles in Ecosystem Function. In Soil Microbiology Ecology and Biochemistry; Eldor, P.A., Ed.; Elsevier Inc.: Amsterdam, The Netherlands, 2015; pp. 111–149. ISBN 9780124159556. [Google Scholar]
  6. Bardgett, R.D. Causes and consequences of biological diversity in soil. Zoology 2006, 105, 367–374. [Google Scholar] [CrossRef] [PubMed]
  7. IPCC. 2014 Climate Change 2014: Mitigation of Climate Change. Contribution of Working Group III to the Fifth Assessment Report of the Intergovernmental Panel on Climate Change; Edenhofer, O., Pichs-Madruga, R., Sokona, Y., Farahani, E., Kadner, S., Seyboth, K., Adler, A., Eds.; Cambridge University Press: Cambridge, UK, 2014. [Google Scholar]
  8. Smith, P.; Nkem, J.; Calvin, K.; Campbell, D.; Cherubini, F.; Grassi, G.; Korotkov, V.; Le Hoang, A.; Lwasa, S.; McElwee, P.; et al. Interlinkages between Desertification, Land Degradation, Food Security and GHG fluxes: Synergies, trade-offs and Integrated Response Options. In Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security, and Greenhouse Gas Fluxes in Terrestrial Ecosystems; Shukla, P.R., Skea, J., Calvo Buendía, E., Masson-Delmotte, V., Pörtner, H.O., Roberts, D.C., Zhai, P., Slade, R., Connors, S., van Diemen, R., Eds.; IPCC: Geneva, Switzerland, 2019; in press. [Google Scholar]
  9. Mulder, C.; Boit, A.; Bonkowski, M.; De Ruiter, P.C.; Mancinelli, G.; Van der Heijden, M.G.A.; Van Wijnen, H.J.; Vonk, J.A.; Rutgers, M. A Belowground Perspective on Dutch Agroecosystems: How Soil Organisms Interact to Support Ecosystem Services. In Advances in Ecological Research; Elsevier Ltd.: Amsterdam, The Netherlands, 2011; Volume 44, pp. 277–357. ISBN 9780123747945. [Google Scholar]
  10. IPCC. Climate Change and Land: An IPCC Special Report on Climate Change, Desertification, Land Degradation, Sustainable Land Management, Food Security and Greenhouse Gas Fluxes in Terrestrial Ecosystems; Shukla, P.R., Skea, J., Calvo Buendía, E., Masson-Delmotte, V., Pörtner, H.O., Roberts, D.C., Zhai, P., Slade, R., Connors, S., van Diemen, R., Eds.; IPCC: Geneva, Switzerland, 2019; in press. [Google Scholar]
  11. Parisi, V.; Menta, C.; Gardi, C.; Jacomini, C.; Mozzanica, E. Microarthropod communities as a tool to assess soil quality and biodiversity: A new approach in Italy. Agric. Ecosyst. Environ. 2005, 105, 323–333. [Google Scholar] [CrossRef]
  12. Menta, C.; Conti, F.D.; Pinto, S.; Bodini, A. Soil Biological Quality index (QBS-ar): 15 years of application at global scale. Ecol. Indic. 2018, 85, 773–780. [Google Scholar] [CrossRef]
  13. Menta, C.; Conti, F.D.; Pinto, S.; Leoni, A.; Lozano-Fondón, C. Monitoring soil restoration in an open-pit mine in northern Italy. Appl. Soil Ecol. 2014, 83, 22–29. [Google Scholar] [CrossRef]
  14. Soong, J.L.; Nielsen, U.N. The role of microarthropods in emerging models of soil organic matter. Soil Biol. Biochem. 2016, 102, 37–39. [Google Scholar] [CrossRef] [Green Version]
  15. Andrés, P.; Moore, J.C.; Simpson, R.T.; Selby, G.; Cotrufo, F.; Denef, K.; Haddix, M.L.; Shaw, E.A.; de Tomasel, C.M.; Molowny-Horas, R.; et al. Soil food web stability in response to grazing in a semi-arid prairie: The importance of soil textural heterogeneity. Soil Biol. Biochem. 2016, 97, 131–143. [Google Scholar] [CrossRef]
  16. Van Dam, N.M.; Heil, M. Multitrophic interactions below and above ground: En route to the next level. J. Ecol. 2011, 99, 77–88. [Google Scholar] [CrossRef]
  17. Menta, C.; Leoni, A.; Gardi, C.; Delia Conti, F. Are grasslands important habitats for soil microarthropod conservation? Biodivers. Conserv. 2011, 20, 1073–1087. [Google Scholar] [CrossRef]
  18. Wardle, D.A.; Yeates, G.W.; Williamson, W.; Bonner, K.I. The response of a three trophic level soil food web to the identity and diversity of plant species and functional groups. Oikos 2003, 102, 45–56. [Google Scholar] [CrossRef]
  19. Lindo, Z.; Whiteley, J.; Gonzalez, A. Traits explain community disassembly and trophic contraction following experimental environmental change. Glob. Chang. Biol. 2012, 18, 2448–2457. [Google Scholar] [CrossRef]
  20. Lacetera, N.; Segnalini, M.; Bernabucci, U.; Ronchi, B.; Vitali, A.; Tran, A.; Guis, H.; Caminade, C.; Calvete, C.; Morse, A.; et al. Climate Induced Effects on Livestock Population and Productivity in the Mediterranean Area. In Regional Assessment of Climate Change in the Mediterranean. Volume 2: Agriculture, Forests and Ecosystem Services and People; Navarra, A., Tubiana, L., Eds.; Advances in Global Change Research; Springer: Dordrecht, The Netherlands, 2013; pp. 135–156. [Google Scholar]
  21. Moreno, G.; Pulido, F.J. The Functioning, Management and Persistence of Dehesas. In Agroforestry in Europe; Springer: Dordrecht, The Netherlands, 2009; Volume 10600, pp. 127–160. ISBN 9781402082719. [Google Scholar]
  22. Fernández, M.P.; Contador, J.F.L.; Schnabel, S.; Gutiérrez, Á.G.; Lozano-Parra, J. Changes in Land Management of Iberian Rangelands and Grasslands in the Last 60 Years and their Effect on Vegetation. In Vegetation; Sebata, A., Ed.; IntechOpen: London, UK, 2018. [Google Scholar]
  23. Lozano-Parra, J.; Maneta, M.P.; Schnabel, S. Climate and topographic controls on simulated pasture production in a semiarid Mediterranean watershed with scattered tree cover. Hydrol. Earth Syst. Sci. 2014, 18, 1439–1456. [Google Scholar] [CrossRef] [Green Version]
  24. Schanbel, S.; Dahlgren, R.A.; Moreno-Marcos, G. Soil and water dymamics. In Mediiterranean Oak Woodland Working Landscapes: Dehesas of Spain and Ranchlands of California; Campos, P., Starrs, P.F., Huntsinger, L., Montero, G., Díaz, M., Eds.; Springer-Verlag: New York, NY, USA, 2013; pp. 91–121. ISBN 978-94-007-6706-5. [Google Scholar]
  25. Pulido, M.; Schnabel, S.; Lavado-Contador, F.; Lozano-Parra, J.; Gómez-Gutiérrez, Á. Selecting indicators for assessing soil quality and degradation in rangelands of Extremadura (SW Spain). Ecol. Indic. 2017, 74, 49–61. [Google Scholar] [CrossRef]
  26. Lozano-Parra, J.; Schnabel, S.; Ceballos-Barbancho, A. The role of vegetation covers on soil wetting processes at rainfall event scale in scattered tree woodland of Mediterranean climate. J. Hydrol. 2015, 529, 951–961. [Google Scholar] [CrossRef]
  27. Moreno, G.; Gonzalez-Bornay, G.; Pulido, F.; Lopez-Diaz, M.L.; Bertomeu, M.; Juárez, E.; Diaz, M. Exploring the causes of high biodiversity of Iberian dehesas: The importance of wood pastures and marginal habitats. Agrofor. Syst. 2015, 90, 87–105. [Google Scholar] [CrossRef]
  28. Pregitzer, K.S.; King, J.S.; Burton, A.J.; Brown, S.E. Responses of tree fine roots to temperature.pdf. New Phytol. 2000, 147, 105–115. [Google Scholar] [CrossRef] [Green Version]
  29. Lozano-Parra, J.; Pulido, M.; Lozano-Fondón, C.; Schnabel, S. How do soil moisture and vegetation covers influence soil temperature in drylands of Mediterranean regions? Water (Switzerland) 2018, 10, 1747. [Google Scholar] [CrossRef] [Green Version]
  30. Pulido, M.; Schnabel, S.; Francisco, J.; Contador, L.; Lozano-parra, J.; González, F. The impact of heavy grazing on soil quality and pasture production in rangelands of SW Spain. Land Degrad. Dev. 2016, 29, 219–230. [Google Scholar] [CrossRef]
  31. Jiménez-Chacón, A.; Homet, P.; Matías, L.; Gómez-Aparicio, L.; Godoy, O. Fine Scale Determinants of Soil Litter Fauna on a Mediterranean Mixed Oak Forest Invaded by the Exotic Soil-Borne Pathogen Phytophthora cinnamomi. Forests 2018, 9, 218. [Google Scholar] [CrossRef] [Green Version]
  32. Fuls, E.R. A technique for objective habitat condition assessments in rangelands. J. Arid Environ. 1992, 22, 195–198. [Google Scholar] [CrossRef]
  33. FAO. FAO-UNESCO Soil Map of the World; FAO: Rome, Italy, 2006. [Google Scholar]
  34. Ministerio de Fomento Centro Nacional de Información Geográfica. Available online: https://www.cnig.es/home (accessed on 31 May 2020).
  35. Drǎguţ, L.; Tiede, D.; Levick, S.R. ESP: A tool to estimate scale parameter for multiresolution image segmentation of remotely sensed data. Int. J. Geogr. Inf. Sci. 2010, 24, 859–871. [Google Scholar] [CrossRef]
  36. Ma, L.; Li, M.; Ma, X.; Cheng, L.; Du, P.; Liu, Y. A review of supervised object-based land-cover image classification. ISPRS J. Photogramm. Remote Sens. 2017, 130, 277–293. [Google Scholar] [CrossRef]
  37. Blaschke, T.; Lang, S.; Lorup, E.; Strobl, J.; Zeil, P. Object-oriented image processing in an integrated GIS/remote sensing environment and perspectives for environmental applications. In Environmental Information for Planning. Politics and the Public, Vol. 2; Cremers, A., Greve, K., Eds.; dc.Publisher: Marburg, Germany, 2000; pp. 555–570. [Google Scholar]
  38. Myint, S.W.; Gober, P.; Brazel, A.; Grossman-Clarke, S.; Weng, Q. Per-pixel vs. object-based classification of urban land cover extraction using high spatial resolution imagery. Remote Sens. Environ. 2011, 115, 1145–1161. [Google Scholar] [CrossRef]
  39. Meloni, F.; Civieta, B.F.; Zaragoza, J.A.; Bautista, S. Vegetation Pattern Modulates Ground Arthropod Diversity in Semi-Arid Mediterranean Steppes. Insects 2020, 11, 59. [Google Scholar] [CrossRef] [Green Version]
  40. LI-COR Inc. LI-8100A Automated Soil CO2 Flux System Instruction Manual; LI-COR Inc.: Lincoln, NE, USA, 2010. [Google Scholar]
  41. Zuur, A.F.; Ieno, E.N.; Elphick, C.S. A protocol for data exploration to avoid common statistical problems. Methods Ecol. Evol. 2010, 1, 3–14. [Google Scholar] [CrossRef]
  42. Borcard, D.; Gillet, F.; Legrende, P. Numerical Ecology with R.; Springer Science & Business Media: New York, NY, USA, 2011; ISBN 978-0-387-78170-9. [Google Scholar]
  43. Anderson, M.J. A new method for non-parametric multivariate analysis of variance. Austral Ecol. 2001, 26, 32–46. [Google Scholar]
  44. McArdle, B.; Anderson, M.J. Fitting Multivariate Models To Community Data . Ecology 2001, 82, 290–297. [Google Scholar] [CrossRef]
  45. Oksanen, J.; Guillaume Blanchet, F.; Kindt, R.; Legendre, P.; Minchin, P.R.; O’Hara, R.B.; Simpson, G.L.; Solymos, P.; Stevens, H.H.; Wagner, H. Vegan: Community Ecology Package, R Package Version 1.15-1; R Foundation for Statistical Computing: Vienna, Austria, 2019.
  46. Zuur, A.F.; Ieno, E.N.; Walker, N.J.; Saveliev, A.A.; Smith, G.M. Mixed Effects Models and Extensions in Ecology with R; Springer Science & Business Media: New York, NY, USA, 2009; ISBN 9781351414234. [Google Scholar]
  47. Wood, S.N. Fast stable restricted maximum likelihood and marginal likelihood estimation of semiparametric generalized linear models. J. R. Stat. Soc. (B) 2011, 73, 3–36. [Google Scholar] [CrossRef] [Green Version]
  48. Wickham, H. Ggplot2: Elegant Graphics for Data Analysis; Springer-Verlag: New York, NY, USA, 2016; ISBN 978-3-319-24277-4. [Google Scholar]
  49. Gao, M.; Liu, D.; Lin, L.; Wu, D. The small-scale structure of a soil mite metacommunity. Eur. J. Soil Biol. 2016, 74, 69–75. [Google Scholar] [CrossRef]
  50. Tilman, D. Niche tradeoffs, neutrality, and community structure: A stochastic theory of resource competition, invasion, and community assembly. Proc. Natl. Acad. Sci. USA 2004, 101, 10854–10861. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  51. Gao, M.; He, P.; Zhang, X.; Liu, D.; Wu, D. Relative roles of spatial factors, environmental filtering and biotic interactions in fine-scale structuring of a soil mite community. Soil Biol. Biochem. 2014, 79, 68–77. [Google Scholar] [CrossRef]
  52. Sha, D.; Gao, M.; Sun, X.; Wu, D.; Zhang, X. Relative Contributions of Spatial and Environmental Processes and Biotic Interactions in a Soil Collembolan Community. Chin. Geogr. Sci. 2015, 25, 582–590. [Google Scholar] [CrossRef]
  53. Gao, M.; Sun, X.; Qiao, Z.; Hou, H.; Lu, T.; Wu, D. Distinct patterns suggest that assembly processes di ff er for dominant arthropods in above-ground and below-ground ecosystems. Pedobiol. J. Soil Ecol. 2018, 69, 17–28. [Google Scholar] [CrossRef]
  54. Caruso, T.; Trokhymets, V.; Bargagli, R.; Convey, P. Biotic interactions as a structuring force in soil communities: Evidence from the micro-arthropods of an Antarctic moss model system. Oecologia 2013, 172, 495–503. [Google Scholar] [CrossRef]
  55. Caruso, T.; Taormina, M.; Migliorini, M. Relative role of deterministic and stochastic determinants of soil animal community: A spatially explicit analysis of oribatid mites. J. Anim. Ecol. 2012, 81, 214–221. [Google Scholar] [CrossRef]
  56. Dong, C.; Gao, M.; Guo, C.; Lin, L.; Wu, D. The underlying processes of a soil mite metacommunity on a small scale. PLoS ONE 2017, 12, e0176828. [Google Scholar] [CrossRef]
  57. Mulder, C.; Den Hollander, H.A.; Hendriks, A.J. Aboveground Herbivory Shapes the Biomass Distribution and Flux of Soil Invertebrates. PLoS ONE 2008, 3, e3573. [Google Scholar] [CrossRef]
  58. Parisi, V.; Menta, C. Microarthropods of the soil: Convergence phenomena and evaluation of soil quality using QBS-ar and QBS-C. Fresenius Environ. Bull. 2008, 17, 1170–1174. [Google Scholar]
  59. Schon, N.L.; Mackay, A.D.; Minor, M.A.; Yeates, G.W.; Hedley, M.J. Soil fauna in grazed New Zealand hill country pastures at two management intensities. Appl. Soil Ecol. 2008, 40, 218–228. [Google Scholar] [CrossRef]
  60. Peco, B.; Navarro, E.; Carmona, C.P.; Medina, N.G.; Marques, M.J. Effects of grazing abandonment on soil multifunctionality: The role of plant functional traits. Agric. Ecosyst. Environ. 2017, 249, 215–225. [Google Scholar] [CrossRef]
  61. Schnabel, S.; Gómez Gutiérrez, A.; Lavado Contador, J.F. Grazing and soil erosion in dehesas of SW Spain. In Advances in Studies on Desertification; International Nuclear Information System: Vienna, Austria, 2009; pp. 725–728, 732. ISBN 978-84-8371-888-9. [Google Scholar]
  62. Junta de Extremadura REDAREX: Red de Asesoramiento al Regante de Extremadura. Available online: http://redarexplus.gobex.es/RedarexPlus/ (accessed on 15 May 2019).
Figure 1. (A) Study area within the Extremadura region (blue boundaries) in the Iberian Peninsula and (B) the result of the object-based image analysis (OBIA) classification for the entire experimental catchment. Red indicates the ensemble of characteristics defining the SHChigh area, yellow represents SHCmedium, and SHClow is indicated by green. Sampling point distribution across the SHC categories is also shown: points sampled in spring are black and points added in autumn are indicated by dark red. A tree symbol (in-C; black or dark red) indicates the geographic location of a sampling point beneath the tree canopy; a grass symbol (out-C; black or dark red) indicates the geographic position of sampling points outside the canopy.
Figure 1. (A) Study area within the Extremadura region (blue boundaries) in the Iberian Peninsula and (B) the result of the object-based image analysis (OBIA) classification for the entire experimental catchment. Red indicates the ensemble of characteristics defining the SHChigh area, yellow represents SHCmedium, and SHClow is indicated by green. Sampling point distribution across the SHC categories is also shown: points sampled in spring are black and points added in autumn are indicated by dark red. A tree symbol (in-C; black or dark red) indicates the geographic location of a sampling point beneath the tree canopy; a grass symbol (out-C; black or dark red) indicates the geographic position of sampling points outside the canopy.
Forests 11 00628 g001
Figure 2. A picture taken in the study area showing the general characteristics of the environments classified as SHChigh, SHCmedium, and SHClow.
Figure 2. A picture taken in the study area showing the general characteristics of the environments classified as SHChigh, SHCmedium, and SHClow.
Forests 11 00628 g002
Figure 3. Boxplots showing the distribution of data for each response variable by soil habitat condition (SHC) and location factor (outside the canopy = Out-C or beneath tree canopy = In-C) in both seasons. (A,B) plots show log-transformed abundances; (C,D) plots show QBS-ar values; (E,F) plots show H′ index values.
Figure 3. Boxplots showing the distribution of data for each response variable by soil habitat condition (SHC) and location factor (outside the canopy = Out-C or beneath tree canopy = In-C) in both seasons. (A,B) plots show log-transformed abundances; (C,D) plots show QBS-ar values; (E,F) plots show H′ index values.
Forests 11 00628 g003aForests 11 00628 g003b
Figure 4. Community composition fits categorical predictors (SHC and location factor), environmental parameters and taxa (abundances and EMIs). Log-transformed abundance matrices for spring and autumn samplings are shown in graphics (A,B), respectively; EMI-values matrices for spring and autumn samplings are shown in panels (C,D), respectively. Black arrows show the fitting of significant taxa, whereas red arrows show the fitting of significant environmental variables. Location indicates whether the sampling point was established beneath the tree canopy or outside the canopy; SHC indicates the characteristics of the surrounding environment, as well as the pressure of the livestock where sampling points were placed; bd = bulk density; ET_avg20 = average evapotranspiration of the 20 days prior to sampling; mo = soil organic matter content; T_0 = mean temperature of the sampling day; T_avg20 = average temperature of the 20 days prior to sampling; T_max0 = maximum temperature of the sampling day; T_min0 = minimum temperatures of the sampling day; T_min.avg20 = average minimum temperature of the 20 days prior to sampling; wc = soil water content.
Figure 4. Community composition fits categorical predictors (SHC and location factor), environmental parameters and taxa (abundances and EMIs). Log-transformed abundance matrices for spring and autumn samplings are shown in graphics (A,B), respectively; EMI-values matrices for spring and autumn samplings are shown in panels (C,D), respectively. Black arrows show the fitting of significant taxa, whereas red arrows show the fitting of significant environmental variables. Location indicates whether the sampling point was established beneath the tree canopy or outside the canopy; SHC indicates the characteristics of the surrounding environment, as well as the pressure of the livestock where sampling points were placed; bd = bulk density; ET_avg20 = average evapotranspiration of the 20 days prior to sampling; mo = soil organic matter content; T_0 = mean temperature of the sampling day; T_avg20 = average temperature of the 20 days prior to sampling; T_max0 = maximum temperature of the sampling day; T_min0 = minimum temperatures of the sampling day; T_min.avg20 = average minimum temperature of the 20 days prior to sampling; wc = soil water content.
Forests 11 00628 g004
Figure 5. GAM plots representing spatial smoothing of the response variables. (A) GAM model for spring total abundances; (B) GAM model for autumn total abundances; (C) GAM model for spring QBS-ar; (D) GAM model for autumn QBS-ar. Codes represent the quantity by which each response variable varies. Black-solid isolines represent the spatial smoothing that belonged to a determined interval of variation. Red-dashed isolines represent the upper variation of that interval associated with each black solid isoline sharing the same code. Green-dashed isolines represent the lower variation of that interval associated to each black solid isoline sharing the same code.
Figure 5. GAM plots representing spatial smoothing of the response variables. (A) GAM model for spring total abundances; (B) GAM model for autumn total abundances; (C) GAM model for spring QBS-ar; (D) GAM model for autumn QBS-ar. Codes represent the quantity by which each response variable varies. Black-solid isolines represent the spatial smoothing that belonged to a determined interval of variation. Red-dashed isolines represent the upper variation of that interval associated with each black solid isoline sharing the same code. Green-dashed isolines represent the lower variation of that interval associated to each black solid isoline sharing the same code.
Forests 11 00628 g005
Table 1. Mean ± standard deviations of soil parameters. SHC indicates the “soil habitat condition” categories; In-C and Out-C indicate whether values were detected beneath the tree canopy or outside the canopy, respectively.
Table 1. Mean ± standard deviations of soil parameters. SHC indicates the “soil habitat condition” categories; In-C and Out-C indicate whether values were detected beneath the tree canopy or outside the canopy, respectively.
SHChighSHCmediumSHClow
UnitsOut-CIn-COut-CIn-COut-CIn-C
Bulk densityg cm−11.5 ± 0.11.2 ± 0.21.5 ± 0.11.2 ± 0.21.5 ± 0.21.2 ± 0.1
Organic matter%5.3 ± 2.510.3 ± 4.61.9 ± 1.89.1 ± 4.63.8 ± 1.49.5 ± 4.1
pH-5.8 ± 0.76.0 ± 0.85.6 ± 0.45.9 ± 0.45.8 ± 0.25.7 ± 0.7
Soil CO2 effluxµmol m−2 s−14.6 ± 2.84.9 ± 2.74.5 ± 1.45.2 ± 1.95.2 ± 1.95.0 ± 2.2
Water content%18.4 ± 5.920.8 ± 8.322.2 ± 8.823.2 ± 8.021.8 ± 8.423.7 ± 10.5
Table 2. Three-way ANOVA on log-transformed abundances, QBS-ar and Shannon’s index (H′). Asterisks indicate levels of significance (*) = p < 0.05; (**) = p < 0.01; (***) = p < 0.001. “Location” indicates whether the sampling point was located beneath tree canopies (in-C) and outside the canopy (out-C).
Table 2. Three-way ANOVA on log-transformed abundances, QBS-ar and Shannon’s index (H′). Asterisks indicate levels of significance (*) = p < 0.05; (**) = p < 0.01; (***) = p < 0.001. “Location” indicates whether the sampling point was located beneath tree canopies (in-C) and outside the canopy (out-C).
Ln AbundancesQBS-arH′
Suorce of VariationF Testp-ValueF Testp-ValueF Testp-Value
Livestock pressure2.9110.0582.4510.0901.5320.220
Location14.655<0.001 ***17.464<0.001 ***0.9130.341
Season7.6440.007 **2.0570.1540.0070.932
Livestock pressure × Location1.2330.2951.0590.3501.2090.302
Livestock pressure × Season1.3550.2620.5110.6012.2470.110
Location × Season0.0150.9024.4900.036 *6.2320.014 *
Livestock pressure × Location × Season1.7180.1841.5810.2100.6340.532
Table 3. PERMANOVA results for matrices of log-transformed abundances and eco-morphological index score matrices. Significant results for environmental parameters causing dissimilarity are marked with asterisk: (*) = p < 0.05; (**) = p < 0.01; (***) = p < 0.001. Location indicates if the sampling point was established beneath the tree canopy (In-C) or outside the canopy (Out-C); SHC indicates the characteristics of the surrounding environment and pressure of the livestock where sampling points were placed; EP (–20) = effective precipitation of the 20 days prior to the sampling day; T (-20) = average temperature of the 20 days prior to the sampling day.
Table 3. PERMANOVA results for matrices of log-transformed abundances and eco-morphological index score matrices. Significant results for environmental parameters causing dissimilarity are marked with asterisk: (*) = p < 0.05; (**) = p < 0.01; (***) = p < 0.001. Location indicates if the sampling point was established beneath the tree canopy (In-C) or outside the canopy (Out-C); SHC indicates the characteristics of the surrounding environment and pressure of the livestock where sampling points were placed; EP (–20) = effective precipitation of the 20 days prior to the sampling day; T (-20) = average temperature of the 20 days prior to the sampling day.
Community MatrixSeasonSource of DissimilarityDfFR2
Log-transformedSpringLocation12.6740.041 *
Abundances Water content13.4440.052 **
pH12.4800.037 *
Soil CO2 efflux12.2630.034 *
Residuals54 0.836
Log-transformedAutumnLocation16.2170.062 ***
Abundances SHC22.1840.044 *
Slope13.1730.032 *
OM content13.5300.035 **
T (−20)12.7200.027 *
EP (−20)12.1560.022 *
Residuals78 0.779
EMIsSpringLocation12.0970.034
Water content13.1550.051 *
Residuals57 0.916
EMIsAutumnLocation111.3290.111 ***
Slope15.0550.050 **
T (−20)13.4560.036 **
Residuals82 0.803
Table 4. GAMs results for community metrics in each season. Significant results are shown in bold. Location indicates where sampling stations were located: outside the canopy (Out-C) or beneath tree canopy (In-C); SHC-low/medium/high indicates the characteristics of environment and pressure of the livestock in which points were located; T (–20) = average temperature of the 20 days prior to sampling day; EP (–20) = effective cumulative precipitation of the 20 days prior to sampling day.
Table 4. GAMs results for community metrics in each season. Significant results are shown in bold. Location indicates where sampling stations were located: outside the canopy (Out-C) or beneath tree canopy (In-C); SHC-low/medium/high indicates the characteristics of environment and pressure of the livestock in which points were located; T (–20) = average temperature of the 20 days prior to sampling day; EP (–20) = effective cumulative precipitation of the 20 days prior to sampling day.
MetricsSeasonParameterFpR2
Log-transformedSprings(NMDS1, NDMS2)21.750<0.0010.874
Abundances Location-out-C7.920<0.001
Location-In-C−0.9480.348
SHClow−1.7050.095
SHCmedium−1.4440.156
pH−1.4160.164
OM content−1.5290.134
Bulk density−1.5360.132
Autumns(NMDS1, NDMS2)23.880<0.0010.891
Location-out-C8.923<0.001
Location-In-C1.7420.087
SHClow1.8990.062
SHCmedium1.0660.291
T (−20)−1.6050.114
EP (−20)2.5390.014
QBS-arSprings(NMDS1, NDMS2)12.870<0.0010.855
Location-out-C35.507<0.001
Location-In-C−0.8730.389
Autumns(NMDS1, NDMS2)37.430<0.0010.918
Location-out-C15.892<0.001
Location-In-C0.1070.915
OM content−2.2140.030
EP (−20)−1.8880.063

Share and Cite

MDPI and ACS Style

Lozano Fondón, C.; González, J.B.; Pulido Fernández, M.; Remelli, S.; Lozano-Parra, J.; Menta, C. Effects of Livestock Pressure and Vegetation Cover on the Spatial and Temporal Structure of Soil Microarthropod Communities in Iberian Rangelands. Forests 2020, 11, 628. https://doi.org/10.3390/f11060628

AMA Style

Lozano Fondón C, González JB, Pulido Fernández M, Remelli S, Lozano-Parra J, Menta C. Effects of Livestock Pressure and Vegetation Cover on the Spatial and Temporal Structure of Soil Microarthropod Communities in Iberian Rangelands. Forests. 2020; 11(6):628. https://doi.org/10.3390/f11060628

Chicago/Turabian Style

Lozano Fondón, Carlos, Jesús Barrena González, Manuel Pulido Fernández, Sara Remelli, Javier Lozano-Parra, and Cristina Menta. 2020. "Effects of Livestock Pressure and Vegetation Cover on the Spatial and Temporal Structure of Soil Microarthropod Communities in Iberian Rangelands" Forests 11, no. 6: 628. https://doi.org/10.3390/f11060628

APA Style

Lozano Fondón, C., González, J. B., Pulido Fernández, M., Remelli, S., Lozano-Parra, J., & Menta, C. (2020). Effects of Livestock Pressure and Vegetation Cover on the Spatial and Temporal Structure of Soil Microarthropod Communities in Iberian Rangelands. Forests, 11(6), 628. https://doi.org/10.3390/f11060628

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop