Soil mite communities (Acari: Mesostigmata) as indicators of urban ecosystems in Bucharest, Romania

The aim of the present study was to establish the effect of management type and of environmental variables on the structure, abundance and species richness of soil mites (Acari: Mesostigmata) in twelve urban green areas in Bucharest-Romania. Three categories of ecosystem based upon management type were investigated: protected area, managed (metropolitan, municipal and district parks) and unmanaged urban areas. The environmental variables which were analysed were: soil and air temperature, soil moisture and atmospheric humidity, soil pH and soil penetration resistance. In June 2017, 480 soil samples were taken, using MacFadyen soil core. The same number of measures was made for quantification of environmental variables. Considering these, we observed that soil temperature, air temperature, air humidity and soil penetration resistance differed significantly between all three types of managed urban green area. All investigated environmental variables, especially soil pH, were significantly related to community assemblage. Analysing the entire Mesostigmata community, 68 species were identified, with 790 individuals and 49 immatures. In order to highlight the response of the soil mite communities to the urban conditions, Shannon, dominance, equitability and soil maturity indices were quantified. With one exception (numerical abundance), these indices recorded higher values in unmanaged green areas compared to managed ecosystems. The same trend was observed between different types of managed green areas, with metropolitan parks having a richer acarological fauna than the municipal or district parks.


Material and methods
Study area. The study was conducted in June 2017, in twelve urban green areas in Bucharest (Fig. 1). Bucharest is situated in Central Eastern Europe, in the lower Danube region. The mean annual temperature is about 10-11 °C, and annual precipitation is 615 mm 52 . Of the total city area (22,800 ha), urban green space occupies 2275 ha (9.97%) while urban parks represent 29.9% of Bucharest and approximately 3% of the total city area but are unequally distributed. Urban parks were classified as (a) metropolitan: MtP (average area = 52.4 ha, over 5000 visitors/weekend day), (b) municipal: MnP (average area = 40.3 ha, 2000-5000 visitors/weekend day), and (c) district: DrP (average area = 6.4 ha, under 2000 visitors/weekend day) 52,53 .
The ecosystems studied were classified in three management scenarios (MS): (1) urban protected area (PA) i.e. Văcăreşti (V), (2) managed green urban areas (MGA) i.e. metropolitan parks [Tineretului (TN) and Titan (TT)], municipal parks (Plumbuita (PL) and Carol (C) and district parks (Floreasca (FL) and Crângaşi (CR), and (3) unmanaged green urban areas (UGA) i.e. Băneasa (B), Pantelimon (P), Griviţa (G), Lacul Morii (LM) and Fundeni (F). A detailed description of these areas is given in Supplementary Appendix 1. Văcăreşti Natural Park is the biggest (183 ha) and the newest park designated in Bucharest as well as the first urban nature park in Romania. It is located at 5 km distance from the city centre and it was classified in the category V ("protected landscape") following the IUCN criteria 54 .
In all managed green areas (with one exception-Văcăreşti Natural Park), the local administrations together with the Plant Protection Centre subordinated to the General Council of Bucharest carries out treatments to combat plant pests (through biological and integrated methods): mites, aphids, thrips, moths, powdery mildew, rust, wasps, defoliating caterpillars, San Jose lice, woolly lice, cicadas, etc. Weed control is also applied in all managed green areas, using manual methods. These methods are governed by the Regulation of the public service for the administration of the public and private domain, regarding the activities of arrangement and maintenance of the parks and of the gardens in the municipality of Bucharest. All managed areas are irrigated.
Soil fauna. Soil samples were collected in June 2017, using a MacFadyen soil core (5 cm diameter) to 10 cm depth. In total, 480 soil samples were collected randomly (40 samples/urban area). In each urban ecosystem, the sampled area was 100 m 2 Due to the small number of immature individuals, they were not considering separately in our statistical analysis. The mites were extracted with a Berlese-Tullgren funnel, in ethyl alcohol, www.nature.com/scientificreports/ clarified in lactic acid and identified to species level, using published identification keys [55][56][57][58][59][60][61][62][63] . The mites extraction lasted 10 days. Some specimens were mounted on permanent slides. All species were deposited in the collection of the Institute of Biology-Bucharest, Romanian Academy-Research Station Posada.
Environmental variables. In order to establish the relationship between soil mite communities and environmental variables, we measured air and soil temperature (T-air, T-soil as °C), air and soil moisture (H-air and H-soil as %), soil pH, and soil penetration resistance (RP as MPa). Air and soil moisture and temperature were measured with a digital thermo-hygrometer PCE-310. Air temperature and moisture were measure at 5 cm above the soil level. Penetration resistance was determined with a soil penetrometer, Step System GmbH, 41,010. The pH was measured with a C532 Jasco Consort pH-meter. The average values of environmental variables are presented in Tables 1 and 2.
Data analysis. The completeness of the species inventory was examined by visual inspection of Mao Tau species accumulation curves for each management scenario and for pooled data of all sites.  www.nature.com/scientificreports/ We examined differences across management scenarios: managed green areas, unmanaged green areas and sites (Băneasa, Pantelimon, Griviţa, Lacul Morii, Fundeni, Carol, Plumbuita, Tineretului, Titan, Floreasca and Crângaşi) using non-metric multidimensional scaling (NMDS) with the Bray-Curtis (B-C) distance index.
The NMDS was conducted using the function "metaMDS" in the R package "vegan" 64. The function automatically applies a square root transformation to the data matrix of species abundance. The environmental variables were standardised to have mean zero and unit variance. We ran NMDS using 500 random starts and tested the goodness of fit of the data using the R2 value and examining the Shepard plot (i.e. the scatter around the regression of the distances between each pair of communities against their original dissimilarities). The significance of differences between communities from management scenarios and sites was assessed with an overall PER-MANOVA based on B-C dissimilarities with the function "adonis" within the vegan package and pairwise using the function "pairwise.perm.manova" within the "RVAideMemoire" package 65 . The same analyses were repeated to examine the differences among different types of MGAs, i.e. metropolitan, municipal and district parks.
To determine whether physical characteristics influenced the species communities, we conducted a canonical correspondence analysis (CCA). Before performing the CCA, we selected the significant environmental variables using the ANOVA permutation test for CCA (999 random permutation) and only introduced the variables with P > 0.05 in the CCA.
We used linear mixed models (LMMs) to test whether the main community feature (i.e. the species richness SR) was related to: (i) management scenario, (ii) environmental variables and (iii) a combination of management scenarios and environmental variables. In the LMMs, the management scenario and the environmental variables were introduced as fixed effects and sites as random effects. We assessed the relative performance of the models using the selection technique based on Akaike's information criterion corrected for sample size-AICc 66,67 . We ranked the models and the model with the lowest AICc was used as the reference for calculating the AIC difference (Δi) and the likelihood of a model given the data and model weights (wi). Models within two AIC units of the AICmin were considered competitive and more plausible than others 66 .
To identify any taxonomic group-specific pattern across management scenarios we examined the extent to which each management scenario was characterised by rare species, i.e. singletons and doubletons. All analyses were performed using R version 3.2.1.
Based on the reproduction strategy of soil mites ("k" or "R") 68 , the index of maturity was calculated for each studied urban green area (Supplementary Appendix 2). Using PAST software, we quantified the following parameters: dominance (D), Shannon index of diversity (H) and equitability (J) 69 .
To identify mite species characteristic to the three MS we applied the indicator Value method using the "indicspecies" R package 70,71 . The method assesses the specificity (uniqueness to a particular MS) and fidelity (frequency of occurrence) of mite species. A high indicator value (IndVal, expressed as percentage) indicates that a mite species can be considered characteristic to a particular or combination of MS. The IndVal for each mite species was calculated based on the matrix of mite species abundance. We used a random reallocation procedure of sites among site groups to test for the significance of IndVal measures for each mite species 70 . According to Dufrêne and Legendre (1997), a mite species is considered characteristic to a MS or combination of MS if the species IndVal is > 25% and significant at p < 0.05.

Results
Analysing the entire Mesostigmata community, we discovered 68 species, with 790 individuals and 49 immatures. If we consider the protected area, 8 species were identified, with 37 individuals and 15 immatures. In managed urban areas 33 species were identified, with 435 individuals and 25 immatures while in unmanaged green ecosystems, we recorded 49 species, with 318 individuals and 9 immatures (Tables 3, 4). None of the species accumulation curves for each management scenario or for pooled data from all sites approached an asymptote, indicating that more samples are required to detect all the species theoretically expected (Fig. 2).
Examining the mean values for environmental variables, we observed that soil temperature, air temperature, air humidity and soil penetration resistance differed significantly between all three types of managed urban green area (p < 0.05). The MGA areas were characterised by the highest mean values of air and soil moisture, RP and by the lowest mean values of air temperature and pH. The PA had the highest mean values of air and soil Table 2. Mean ± standard deviation (range) of the environmental variables for the different categories of managed green areas. The p value represents the significance level of the global one-way analysis of variance performed separately for each variable. www.nature.com/scientificreports/ www.nature.com/scientificreports/ temperature, pH and moderate mean values of the remaining abiotic factors. In UGA, the lowest mean values of soil temperature, RP air and soil moisture were recorded (Table 1). Measuring the environmental variables from the three types of MGA, we recorded a significant difference between their mean values (p < 0.05). MtP was characterised by the highest mean values of air and soil temperatures, soil moisture and the lowest mean value of pH. MnP were characterised by the highest values of atmospheric humidity (RP) and lowest air temperature. In DrP, the lowest average values of soil temperature and moisture (RP) and the highest mean value of pH were observed ( Table 2).
The Shepard plot showed that original dissimilarities were well preserved in the reduced number of dimensions and low stress was found for both NMDS analyses examining the differences across the MS (R2 = 1, stress = 0.001) and different types of MGA (R2 = 0.966, stress = 0.055) (Supplementary Appendix 3), respectively.
The NMDS analysis showed low clustering of data points by both management scenarios, types of MGA and sites (Supplementary Appendix 4).
PERMANOVA showed significant differences between MS (F [2,243] = 3.77, R2 = 0.03, p < 0.001), MGA (F [2,106] = 2.91, R2 = 0.050 = 3, p < 0.001) and sites (F [11,243] = 3.04, R2 = 0.126, p < 0.001). Pairwise comparisons of species communities identified marginally significant differences only between MGA and UGA (p = 0.05) and DrP and MtP (p = 0.05). Comparisons between species assemblages between all pairs of sites are shown in Table 5.   www.nature.com/scientificreports/ The CCA model including all the environmental variables explained a total of 52.00% of the overall variation in species composition, of which 26.01% was explained by the first axis and 25.9% by the second axis, respectively. The first canonical axis was highly correlated with pH (− 0.944) whereas the second canonical axis was correlated with T-soil (0.710), T-air (− 0.598) and H-air (0.712). The biplot classification on the two first axes showed species from the upper right quadrant were associated with H-soil, Rp, T-soil and H-air. The lower quadrant contains species associated with T-air. Species associated with pH were placed in the upper left quadrant (Fig. 3).
The ANOVA permutation test for CCA showed that all the environmental variables were significantly related to community assemblage (Supplementary Appendix 5).
Based on species richness, the model selection using AICc indicated one highly supported model that included the pH ( Table 6). The second best supported model included the MS. Species richness significantly increased with pH (F [1,231] = 4.00, p = 0.04) for all management scenarios (Fig. 4). No significant differences in species richness were found among the UGA (F [2,106] (Table 3).
Turning to the dominance index, the mite communities from TT, V and PL are represented by a few dominant species: Pergamasus quisquiliarum and Arctoseius cetratus in TT, Rhodacarellus silesiacus in V, and Hypoaspis aculeifer in all three ecosystems. In the remaining green areas studied, representation of remaining species was almost equitable. This fact was indicated by the equitability index, characteristic for each urban area (Table 3).
Comparing the taxonomic structure of the mite communities, we observed that each managed urban green area was defined by some characteristic species i.e. (a) for the natural park: Pergamasus crassipes and Parasitus fimetorum, (b) for metropolitan parks: Cherosieus bryophilus and Leioseius minusculus, (c) for municipal parks: Paragamasus similis, Protogamasellus singularis and Olopachys vysotskajae, and (d) for district parks: Pergamasus laetus, Parasitus beta, Pseudolaelaps doderoi and Laelaps astronomica.
In the different types of managed parks in Bucharest, we observed that in MtP, the number of species as well as numerical densities is higher than in DtP. In MtP the average value of soil moisture is higher than in DtP (15.63 in comparison with 10.84). The mean value of soil temperature is higher in MtP (18.06 °C) that in DtP (13.41 °C). These environmental variables could influence the mite communities from the two types of managed urban areas. Table 6. Akaike statistics for model including the species richness, the management scenario and the environmental variables, AIC (Akaike's Information Criterion) differences (ΔAICc) and Akaike weights (wi) were used to rank models relative to the best model (minimum AIC). K number of parameters, LL log likelihood.  www.nature.com/scientificreports/ Making an ecological analysis of the mite communities from UGA areas, the highest number of identified species was recorded in B park (22 species), as well as the Shannon index of diversity (H = 2.791) and equitability (J = 0.903). The highest value for numerical abundance was obtained in park F (110 individuals). In contrast, park P has lower values for all these parameters, with the exception of dominance, whose value was D = 0.294 and of maturity index (M = 0.71) ( Table 4).

Discussion
Comparing the present research with data obtained from other studies in Europe on the Mesostigmata of urban soils, the number of species identified in Bucharest is comparable with totals obtained in Latvia and Poland, in such different habitats as roadsides, greenery, parks and urban forests (25-51 species), but is higher than that for mite communities in a botanical garden in Hungary (12 species) 40,42,45,47 . However, the number of recorded species is lower than suburban and urban green areas from Poland (62-68 species), urban grasslands from Latvia or cemeteries and botanical garden in Slovakia (123 species, in three years of study) 34,35,38,39,43 .
Analysing the above presented results, we observed that management intensity in urban green areas, from Bucharest-Romania, impacts mesostigmatid mite diversity and community composition. The effect of management intensity on mesostigmatid mite diversity was highlighted by the results of pairwise comparison of species communities that revealed significant differences between MGA and UGA areas. The number of recorded species in UGA is higher than that for MGA. The Romanian results agree with those obtained in Latvia for Riga city, where there was a high diversity of soil invertebrates in the disturbed urban forest habitats, but undisturbed soils harbour a greater species richness of mites than disturbed soils 40 .
The management type of the urban green areas influenced the species composition as well. Species composition differs in the three management scenarios for urban green areas. About 27.54% of the total number of species were characteristic of MGA, with Pergamasus quisquiliarum, Dendrolealaps foveolatus and Gamasellus bicolor having the highest recorded abundance. Other research has indicated that Pergamasus quisquiliarum is a species with affinity for parks and street green areas in Warsaw, Poland. Dendrolealaps foveolatus and Gamasellus bicolor have been identified in urban forest in Riga city and in suburban green areas in Warsaw 35,45 .
In UGA, 47.83% are species identified only in these parks. In these mite communities, species from the Macrochelidae and Pachylaelapidae families were best represented, with Macrocheles recki, Aliphis halleri and Parasitus fimetorum most numerous. These species were also associated with industrial and post-industrial areas in Poland 38 .
Comparing the mite communities of PA with those from the remaining MGA, we observed not only the lowest number of species and numerical densities in PA but also the lowest value of the Shannon index of diversity (H = 1.696). The same situation is found if we compare PA mite communities with all communities from unmanaged urban green areas. These results are supported by the MANOVA pairwise comparisons of species communities among sites (Table 5). Analysing abiotic factors, we observed that in the PA, the highest mean values of soil and air temperature were recorded (16.87 °C and 30.02 °C), potentially affecting the structure of mite communities. Despite the fact that this PA (Văcăreşti) is wetland, higher temperatures during summer and autumn cause high evapotranspiration and water depletion at depth in the soil 72 .
Comparison between species communities from MGA revealed that those from DrP and MtP were significantly different (according to PERMANOVA test). In managed types of district park green areas the number of species as well as numerical densities are higher than in metropolitan parks. Some environmental variables are significantly different between parks e.g. soil moisture and temperature. There are also differences in community composition structure. All these differences are correlated with specific environmental conditions. Higher soil humidity and temperature provided more favourable environmental conditions for predatory soil mites in metropolitan parks. Drought or lower soil humidity has a negative impact on soil microarthropod fauna, both in general and on mites specifically 73 .
Using the presented information, based on the obtained results from this research, we could propose that the future new created urban parks should be designed taking into account the natural/unmanaged model. This means that, from acarological point of view, the autochthone vegetation and native soils constitute a proper habitat for Mesostigmata species. In unmanaged urban green areas, the ecological process are natural and undisturbed (not anthropized as in managed parks) and these ecosystems will follow the natural ecological succession. If we put into discussion the relation between ecological characteristics of the managed urban green areas and Mesostigmata fauna, we consider that their management could be reconsidered. This reconsideration means to use native species of plant, trees and soil, reach in organic matter (not indigenous brought from Asia or Central America, usually accompanied by the alochtonous soil and other invasive species).
Research studies from all over the world revealed that the Mesostigmata communities correlated with environmental variables could constitute a very useful biological tool to characterized the ecological quality of terrestrial ecosystems and the quality of soils 9,10,18,21,[25][26][27][29][30][31]33 . In the case of urban green areas from Bucharest, under different management scenarios, their ecological characterization was made taking into account the variation of the environmental parameters, in correlation with spatial dynamics of the communities' parameters of Mesostigmata fauna (species diversity and numerical abundance). www.nature.com/scientificreports/ Analysing the environmental parameters, the highest average value of soil moisture (13.16%) was recorded in MGA areas where the soil penetration resistance is higher (RP = 1.8). However, the differences between mean values of soil moisture for all the green areas are not significant (p = 0.334), due to the fact the soils throughout Bucharest are mainly clay ( Table 1). Because of their very small particle size and consequent large surface area, clays are able to retain greater amounts of water than sandy or most anthropogenic soils 74 . Other research studies revealed that there was a strong correlation between soil penetration resistance, pore space, soil moisture and quantity of organic matter 75,76 .
In general, an increased value of RP means a decreased quantity of organic matter, a lower value of soil moisture and reduced pore space. Although soil moisture in our study is higher in MGA than in UGA or PA, it is possible that the lack of organic matter is the main limiting factor for species diversity. However, clay particles from MGA soils had a much greater tendency to stick together than sands (which were present mainly in UGA), and not only did soil moisture have recorded high values but also the RP is higher 74 .
In the investigated plots, especially those in UGA, the dominant species of Mesostigmata are predatory mites. Their presence is linked to organic matter (correlated with pH), which represents a trophic reservoir and favourable habitat for other invertebrate groups-a relationship supported by the higher values of the soil maturity index 15,28,77,78 . The relationship between clay soil and pH is well established. Clays have a large specific surface and are predominantly negatively charged, retaining nutrients against leaching and buffering the soil against extreme pH changes 79 . It has also been discovered that on basic soil some species of bacteria develop, with these two variables being linearly correlated 80 . Bacterial community development constitutes a favourable factor for the soil invertebrates that represent the food source for predatory mites 9,77 . This may also explain the linear relation between pH value and species richness obtained in our research.
Soil structure is undisturbed in UGA compared with MGA; in MGA there is some modification every year due to the establishing of different ornamental plants. According to specialists 81 , some soils in Bucharest were identified that, although less modified, nonetheless received (and were still receiving) mainly negative impacts caused by daily household or industrial activities. Soils of this type were found in those green spaces less modified by urbanism, and soils from peripheral and suburban areas.
If we consider the influence of the environmental variables on all the mite species identified, we observed that species such as Dendrolealaps foveolatus, Pergamasus quisquiliarum, Laelaps astronomica, Hypoaspis praesternalis and Onchodellus karawaiewi were influenced by soil moisture and air humidity, soil temperature and RP. In Europe, these species have been identified mostly in grasslands, grassy arable fallows or urban habitats, their communities being correlated with soil moisture (Onchodellus karawaiewi being a euryhygrophilous species) or organic matter content. In our study, this explains the relationship between such mites and RP 20,25,44,60,79 .
Three species (Rhodacarus denticulatus, Olopachys suecicus and Glyptholaspis americana) were influenced by soil pH. Soil pH itself is influenced by different vegetation types, including those from urban green areas 20,82 . Olopachys suecicus is eurytopic, with wide ecological plasticity, able to colonise poorly shaded habitats (urban parks, scrub, grasslands). Glyptholaspis americana is cosmopolitan and mostly found in acid habitats such as dung or leaf litter 58 .
Rather surprisingly, in these urban ecosystems air temperature influenced the composition of mite communities indirectly. In such ecosystems, the most striking difference from surface habitats was that the moisture/ humidity and temperature extremes were much less pronounced in the soil. Similar patterns have been identified in other habitat types, such as screes 83 .
Spatial dynamics of the investigated parameters of soil mite communities related to environmental factors could constitute a response of this invertebrate group to the specific conditions of urban green areas under different types of management scenario.

Conclusions
The presented study established the effect of management type and of environmental variables on the structure, abundance and species richness of soil mites (Acari: Mesostigmata). For the first time in Romania, an intensive study on urban Mesostigmata fauna has been made. Twelve urban areas were characterised by specific environment conditions and mite communities. Six environmental factors were investigated (soil and air temperature, soil moisture and atmospheric humidity, soil pH, soil penetration resistance). In total 480 soil samples were analysed for the soil fauna, as well for the abiotic factors.
Examining the mean values of environmental variables, we observed that soil temperature, air temperature, air humidity and soil penetration resistance showed significant differences between all three types of managed urban green area. All environmental variables, especially soil pH, significantly influenced the mite communities. The MGA areas were characterised by the highest mean values of air and soil moisture, the PA had the highest mean values of air and soil temperature, pH and the UGA by the lowest values of these parameters. From all collected soil samples, 68 Mesostigmata species were identified, with 790 individuals and 49 immatures. In terms of the number of species and numerical abundance in relation to management regimes and to the specific investigated environmental conditions, differences in mesostigmatid mite communities were recorded between managed green areas (MGA) and unmanaged green areas (UMG). In comparison with MGA, UGA were characterised by higher values of the community parameters (Shannon diversity, dominance and equitability), as well as by the highest values of the soil maturity index. Comparative analysis between the three types of MGA revealed that the species communities from metropolitan parks were richer than those from district parks.
The investigated indices that showed different values in relation to environmental variables demonstrate important links between mite communities in specifically urban ecosystems that are under anthropogenic pressure. On the other hand, we consider that unmanaged urban green areas as "hotspots" of Mesostigmata diversity. www.nature.com/scientificreports/ This will constitute an argument of the usage of native vegetation and soil, for creation of new urban parks or for management of those existing ones.

Data availability
The data sets analysed during the current study are available from the corresponding author on reasonable request.