Soil properties of North Iberian wet heathlands in relation to climate, management and plant community

Heathlands are a broad vegetation type characterized by the dominance of evergreen shrub species that thrive on nutrient-poor soils, thus sensitive to small changes in soil conditions. Here we aimed to identify soil gradients related to climate, management and main species in wet heathlands. Soil nutrient levels, organic matter and acidity were studied at two soil depths on ninety plots from 18 sites in Erica mackayana wet heaths of NW Iberian Peninsula, in relation to presence and cover of structural plant species (shrubs and graminoids) at two scales: plot and site (landscape) scale. We identified one main soil gradient explained by soil organic matter (SOM), the effective cation exchange complex (eCEC), available phosphorus (P), and Aluminium-Calcium ratio (Al:Ca). Cattle density had a positive correlation with the main gradient in the surface layer, all other climate and management factors were unrelated to soil conditions. Molinia caerulea had a positive relation with SOM, eCEC, basic cations and low Al:Ca ratio. Erica cinerea showed the opposite reaction at both scales. Ulex gallii showed a negative correlation with C:N ratio at the plot scale and deeper layer. SOM accumulation, low nutrient levels and Al toxicity explain the uniqueness of E. mackayana heathland vegetation and soils. Main indicator species react to soil conditions at plot and landscape scales. Cattle density correlated positively with soil nutrient levels, but density of wild ponies is unrelated to soil conditions. Large herbivores, especially ponies, are needed for conservation management of these heathlands.


Abstract
Background and aims Heathlands are a broad vegetation type characterized by the dominance of evergreen shrub species that thrive on nutrient-poor soils, thus sensitive to small changes in soil conditions. Here we aimed to identify soil gradients related to climate, management and main species in wet heathlands. Methods Soil nutrient levels, organic matter and acidity were studied at two soil depths on ninety plots from 18 sites in Erica mackayana wet heaths of NW Iberian Peninsula, in relation to presence and cover of structural plant species (shrubs and graminoids) at two scales: plot and site (landscape) scale.
to small changes in soil acidity and the concentration and availability of nutrients (Braakhekke and Hooftman 1999;De Graaf et al. 2009;Field et al. 2017;van der Wal et al. 2003). These changes may drive important shifts in vegetation structure, community composition and species richness and diversity (Bobbink et al. 1998;Fagúndez 2013;Tibbett et al. 2019).
Heathland soils are highly organic, acidic and with limited nutrient loads (Gimingham 1981;Zas and Alonso 2002). These are soils of an aluminous, podzolic or gleyic type, which can develop as a catena along the slope. In the highest or steepest parts, where erosive processes and extractive dynamics dominate, with intense drainage and nutrient desaturation, skeletal soils develop, little evolved, with polycyclic features and A umbric or H hystic horizons. In the less energetic sectors of the landscape, the accumulation of a large amount of organic matter favors the slowing down of drainage and the development of more or less intense hydromorphic conditions, a reduction in microbial activity and the appearance of peaty soils.
An important underlying gradient in wet heaths is the amount of Soil Organic Matter (SOM), that functions as a stock of nutrients for plants, providing the largest pool of macronutrients like Nitrogen (N), Sulphur (S) and Phosphorus (P) (Trivedi et al. 2018). SOM is directly related to the effective Cation Exchange Complex (eCEC) that retains basic cations needed for plants. However, the acid soils are characterized by nutrient deficiency and toxicity of metals such as manganese (Mn), iron (Fe) and aluminium (Al); with Al toxicity being the main factor limiting plant growth in this type of soil (Bose et al. 2015;Gupta et al. 2013;Kochian et al. 2004). The eCEC is partly covered by protons (H + ), therefore Base Saturation (BS) of other cations like Calcium (Ca) is low. In addition, part of the eCEC is occupied by Al limiting basic cations availability. Aluminium is mobilized at low pH values (<5.5) increasing its solubility in suboxic environments (Takeno 2005), a highly toxic element-species that constrains plant development in acid soils (De Graaf et al. 2009;Houdijk et al. 1993;Kleijn et al. 2008). A high tolerance to Al toxicity in heathland species like Calluna vulgaris is an advantage for competition in such limiting environments (De Graaf et al. 1997;Vogels et al. 2013). Al toxicity is related to high levels of acidification, when other cations such as Ca and magnesium (Mg) are mobilised, allowing Al to replace them in the eCEC (De Graaf et al. 1997;Vogels et al. 2020). The high levels of Al:Ca ratio in heathland soils is a good indicator of this gradient (De Graaf et al. 2009). Iron (Fe) is mobilized in highly acid soils, and may become a toxic element at high concentrations. Ericoid mycorrhizal endophytes aid in regulating Fe availability for heathers such as C. vulgaris ). In general, ericoid mycorrhiza buffers metallic elements intakes, an advantage over other species with different strategies (Read 1996).
Soil acidification and presence and availability of nutrients have shown to play a major role in community composition and habitat quality of northern European wet heaths (Damgaard et al. 2014;Nielsen et al. 2017;van Voorn et al. 2016). A major shift of heathlands transformation into grasslands has been described in northern European heaths since the second half of the twentieth century, changes that have been linked to high levels of atmospheric deposition of different forms of Nitrogen (NHx, NOx) and Sulphur (SOx) (Bobbink et al. 1998;Heil and Diemont 1983;Hofland-Zijlstra and Berendse 2010). This chronic supply causes acidification in the system because N saturation causes basic cations leaching and increases free protons concentration in the soil solution (Grennfelt and Hultberg 1986). It also increases N soil availability compared to P or K, and subsequently stoichometric unbalance of N:P and N:K ratios on plant tissues (Bobbink et al. 1998;Roem and Berendse 2000;de Vries et al. 2009). These changes can trigger processes of habitat degradation such as higher herbivory rates and disruption of interspecific competition (Vogels et al. 2013). In heathlands, grasses, ericoid shrubs and other species like gorse compete for scarce resources. Under the restrictive conditions of soil acidity and low nutrient levels, evergreen woody species with ericoid leaf types commonly associated to ericoid mycorrhiza become dominant, but these species have a low competitive ability in richer environments (Aerts 1990). In general, higher nutrient concentrations favours grasses like Deschampsia flexuosa and Molinia caerulea over shrubs like C. vulgaris or Erica tetralix (Aerts 1989;Alonso et al. 2001;Damgaard et al. 2017). For example, small variations in pH and basic cations (Ca, Mg) availability favours M. caerulea in wet heaths (Riesch et al. 2018).
In the northwest of the Iberian Peninsula, heathland occupies an area of at least 600 k ha as dominant formations, together with a significant extension of underwood formations that have not been inventoried. The specific Erica mackayana Atlantic wet heathland is a habitat of priority conservation in the European Union (European Comission 2007). It covers large areas in the mountains of northern Galicia, in the north west of the Iberian Peninsula (Díaz-Varela et al. 2018;Fagúndez 2016Fagúndez , 2018Fraga et al. 1990). This is a unique habitat that support high levels of biodiversity of different groups such as plants, insects, amphibians, birds and soil fauna, and provide many ecosystem services such as water supply and carbon storage. In northwest Iberian, the community develops under a range of soil types, from mineral to peaty soils, developed over granite, schists, or even serpentine or limestone (Fagúndez 2006(Fagúndez , 2016García-Arrese et al. 2009;Loidi et al. 2010;Pillon et al. 2019 (Fagúndez 2016) and interspecific competition (Fagúndez 2018), along other topographic, climatic and soil factors. Traditionally, these heathlands have been used as rangelands for livestock, with associated management practices to prevent shrub encroachment. The abandonment of the traditional management practices has led to a continuous loss of the area covered by heaths since the 1960s, changed into afforested areas of pine (Pinus spp.) and, more recently, eucalyptus (Eucalyptus spp.), improved pastures, or successional Ulex-Cytisus scrub followed by the formation of deciduous forests of birch (Betula alba) and oak Trophic interactions and ecosystem functioning in wet heaths are modulated by the effects of large herbivores, mainly wild ponies. Albeit it constitutes the largest European populations of free living horses with over ten thousand animals, the Galician wild ponies have been poorly studied (Nuñez et al. 2016;Fagúndez et al. 2021). Wild ponies have a positive impact on vegetation structure, species diversity, and other ecosystem services, but this traditional low intensity management is being abandon due to depopulation, constraining policies and other factors (Fagúndez 2016;López-Bao et al. 2013). On one hand, absence of ponies causes shrub encroachment, particularly by gorse (Ulex spp.), what constrains cattle grazing and increases the risk of wildfires (Muñoz-Barcia et al. 2019). On the other hand, over grazing by cattle and related activities such as manure addition enhances grasses and result in a homogeneous herbaceous community with a lower diversity of rare species (Fagúndez 2016;González-Hernández et al. 2020;López-Bao et al. 2013;Silva-Pando et al. 2002).
In contrast to northern countries, soil properties and nutrient fluxes have scarcely been investigated in southern European Atlantic wet heaths. In this study we provide the first regional approach to species responses to the main parameters, including soil organic matter (SOM), soil acidity (pH), effective cation exchange complex (eCEC), main nutrients represented by phosphorus (P) and basic cations (Ca, Mg, K, Na), toxic elements such as aluminium (Al), and two ratios (C:N and Al:Ca) of the Erica mackayana wet heaths in northern Galicia. Analyses were performed at two surface soil depths (between 0 and 20 cm depth) on ninety plots from 18 sites throughout the region, and were related to plant species presence and abundance at two scales: plot and site (landscape) scale. We aimed to describe the soil condition of the community, and to evaluate vegetation trends in relation to management. The general drivers of habitat quality loss such as gorse expansion and dominance, and Molinia caerulea encroachment, were particularly examined in the light of acidification and P availability, and tolerance to Al in the soil. We aimed to i) describe the main soil gradients for the habitat, ii) evaluate the relation between environmental constrains (climate, herbivory) and such gradients, and iii) study the specific response of the main species to identify reliable indicators of soil conditions. These results will provide the knowledge base for the establishment of regional patterns of soil nutrients in wet heaths, and to evaluate potential changes in its floristic composition under higher densities of cattle livestock and fewer wild ponies, or potential threats such as increasing atmospheric pollution. Appropriate management for conservation of the habitat must consider the plant-soil system to avoid impoverishment of the community and maintain habitat quality.

Study area
Eighteen heathland sites were prospected in July-August 2012 and 2013 in the north of Galicia, covering as much as possible the regional variation range in location throughout the region, elevation, lithology and management ( Fig. 1). Details on the study sites can be found in Fagúndez (2016) and supplementary  Table S1.
Climate data were obtained from the digital cartography services of the University of Extremadura (http:// secad. unex. es; Felicísimo 2011), with a resolution of 1 km × 1 km. Values at the centre of each field studied were extracted for climate and elevation. Seasonal temperature ranges and precipitation are the main climatic constraints for the different heathland types (Loidi et al. 2007), therefore we included mean temperature of the coldest month (January), mean temperature of the hottest month (July), summer precipitation (June+July+August) and annual precipitation (Fagúndez 2016). We also included altitude for each site. Site management was classified as grazed or ungrazed (binary). The main livestock species were cattle and/or ponies. For each species, we defined an index of animal density based on dung counts in a fixed area of 50 m 2 , described in detail in Fagúndez (2016). Only one site was grazed by sheep, therefore it was not analyzed independently.

Floristic data
The floristic assessment of plant frequency and cover was calculated in two complementary surveys with two datasets: i) a presence-absence matrix from five 2 × 2 m plots at each study site (hereafter plot dataset), and ii) a species cover matrix from three transects at each study site (hereafter site dataset). The five plots were set randomly at a minimum distance of 25 m from field edge and from each other. Species occurrence was recorded for all vascular plants. From this data, we created a presence/absence matrix for each species. For the site scale, ten-meter transects were randomly located in each study site, perpendicularly to the main slope. In these transects, a vertical marked bar of 2.5 cm width was used to calculate vegetation cover and height. Every 25 cm, the bar was vertically lowered and every species hit on vegetative organs (i.e. stems or leaves) was recorded together with height to the nearest cm. If a species were hit more than once, only the first (i.e. tallest) hit value was recorded. Species cover was calculated as the number of positive hits for the species divided by the total number of measures in the transect (40). To counteract the neighbouring effect (i.e. high probability of hitting the same species in proximity after a positive hit), we performed three independent transects, at a distance of 25 m or more, and then calculated the arithmetic mean cover value for each species and site (Fagúndez 2016(Fagúndez , 2018. Data were collected for all structural species: shrubs, grasses and sedges. Only those species recorded in at least four sites were included in the analyses. Subsequently, six recorded species were not included: Erica tetralix (one site), Erica vagans (2), Ulex europaeus (1), Thymelaea coridifolia (3), Carex panicea (1) and Luzula multiflora (1). The complex group of stoloniferous Agrostis species which includes A. canina, A. hesperica and A. capillaris were pooled due to difficulties on species identification under A. canina s.l. A final set of seven shrubs and nine "graminoids" (seven grasses and two species of Carex) was obtained.

Soil sampling and analyses
In the approximate centre point of each plot (18 sites, five plots per site), we collected topsoil samples of approximately 20 cm depth, and placed them in a plastic bag labelled and stored at 4 °C until processing. Soils were taken to the lab, and allowed to air-dry for 48 h. Soils showed a general profile of a dark organic or mineral very rich in organic material layer of 3-6 cm (horizon H/O/A), hereafter named Surface Horizon (SH) followed, with variable depth, by a more mineral, but still rich in organic matter, layer (horizon A) and, in a few places, hydromorphic organic layers (horizon H), hereafter named Sub-Surface Horizon (SSH). These two horizons were identified following the general guidelines for soil description (FAO 2006), including boundaries topography and transition of boundaries, presence, nature and abundance of sands and gravels, features of aeromorphic organic layers, color, hydromorphism, presence and color of mottles, structure and consistency of aggregates, porosity, abundance and size of roots.
A wide range of preliminary soil analyses for general description included texture, colour, bulk density or water content and a range of chemical parameters (supplementary data Table S2).
Analytical determinations were performed on airdried soil samples sieved through a 2 mm mesh. Only for those determinations that required it, the samples were dried at 105 °C and/or ground in an agate mill to a size of less than 5 μm. All soil samples were analyzed in duplicate and repeated when the coefficient of variation exceeded 10%. Only those selected for subsequent analyses are described.
The soil pH in water (pHw) and 0.1 M KCl solution (pHK) were determined using a 1:2.5 soil:solution ratio (Guitián and Carballas 1976). The pH in water reflects the actual acidity and the pH in KCl the potential acidity changeable in the presence of ions.
Soil Organic Matter (SOM), in %, was obtained by the Loss-On-Ignition (LOI) procedure, in which dried soil fractions (105 °C to constant weight) were burned in the oven at 550 °C. In the absence of carbonates or hydrated clay minerals, the lost weight was identified as the organic matter content of the soil.
Total C and N contents were analyzed in finely milled samples using an elemental autoanalyzer FISONS EA 1108 and infrared detection, with an absolute error of less than 0.3%. The detection limit for C is 50 mg kg −1 and for N is 100 mg kg −1 . The equipment was calibrated and verified using different reference materials (EDTA-502-092, 502-309-SOIL, SOIL-502-308). Total C is considered organic C because of the lack of carbonates in these soils.
Plant available phosphorus (P) extracted with a 0.03 M NH 4 F plus 0.1 M HCl solution (Bray and Kurtz 1945) was determined colorimetrically as ortophosphate with an ascorbic acid reduced molybdate-blue complex.
Effective cation exchange capacity (eCEC) was estimated (Gillman and Sumpier 1985) as the sum of the exchangeable basic cations (K, Na, Ca, Mg) displaced with 1 M NH 4 Cl (Peech et al. 1947) and exchangeable Al, Fe and Mn extracted with 1 M KCl (Lin and Coleman 1960). In extracts with pH less than 4, the H + content was determined by titration with 0.01 M NaOH at end point 4.0. Base saturation is the sum of the extractable basic cations including Ca, Mg, K, Na, expressed as a percentage of eCEC.

Data analysis
Many soil parameters are highly correlated, which causes multicolinearity issues in standard analyses (Sena et al. 2002). We performed bivariate correlation analyses by calculating Spearman's rho correlation index, and Principal Component Analyses (PCA) including different combinations of parameters for better resolution. Based on this preliminary analyses and general assumptions on key elements, we included the following parameters in the analyses: SOM, pH w (in distilled water), C:N ratio, extractable P, Al:Ca ratio, eCEC, Basic cations Saturation (BasS: Ca, Mg, Na, K in eCEC).
The Kruskall-Wallis H-test (H k-w ), a non-parametric procedure for several independent samples, was also applied, taking as modulators of soil properties the environmental variables orientation and slope and the type of horizon.
Principal component analysis (PCA) was used to evaluate the relevance of climate and herbivory in the soil. Principal components (PCs) with eigenvalues higher than one were retained and subjected to varimax rotation with Kaiser normalization to obtain the proportion of the variance of each soil parameter explained by PCs. The correlation of the management (grazing and herbivore species and densities) and environmental factors (climate and elevation) with the main PCs (eigenvalues >1) was measured using Spearman rho.
At the site scale we used frequency of each species as the mean cover of the three transects and measured correlation with above-mentioned PCs using Spearman rho. Since some species were absent from a large number of areas, correlation was also calculated for presence-only to address the weight of low-frequency species (not shown). At the plot scale we evaluated specific response of each species to each soil parameter independently. We used the presence/absence matrix and performed non-parametric tests (Man-Witney-Wilcoxon test) to compare values of each soil parameter in plots with and without the species under analysis.

Soil description
The studied soils are dark, well-structured with high total porosity, with high total values of SOM and low nutrient concentrations. Mean values and variability for each parameter, and contrasting values between SH and SSH are shown in Fig. 2 and supplementary material Table S2.
Soil organic matter (SOM), determined from LOI, presented higher values in SH (H k-w ; p < 0.001). It showed a strong positive correlation with C and N (r 2 : 0.936 and 0.885; p < 0.001, n: 102). The organic fraction was poorly mineralized, as depicted by the high C:N ratio (Fig. 2). The C:N ratio values determined in the studied soils are moderately high, and are similar to those found in other studies on heathland soils (Berendse 1990;Certini et al. 2015;Thaysen et al. 2017). These C:N ratios are identified with a soil organic matter with many fresh remains, where the plant source can be recognized, and a smaller proportion of more evolved material, represented by a mor a moder humus. However, their real nature will require future characterization with more specific techniques. The soils analyzed showed an acidic reaction with pH mean values in water close to 4, and 3.3 in KCl. This parameter is also controlled by SOM, depicted by a negative correlation with pH (r 2 : 0.616, p < 0.001, n: 101). The mean P concentration was 35.9 mg kg −1 , although a strong variability was detected among samples, with minimum concentrations lower than 2 mg kg −1 (Fig. 2). Available P maintains significant and positive correlations (p < 0.001) with LOI, C, N, and particularly with the C:N ratio. The eCEC was low, with an average of 7.5 cmol c kg −1 and a marked variability among samples clearly controlled by pH (r 2 : 0.503; p < 0.001; n: 101). Statistically significant positive correlations were also observed between pH and basic cations, in particular Ca (r 2 : 0.462; p < 0.001; n: 100) and negative with Al (r 2 : −0.619; p < 0.001; n: 102). The eCEC is higher when the exchange complex is dominated by basic cations, represented by BasS (r 2 : 0.947; p < 0.001; n: 101) and reaches minimum values when SAl rises (r 2 : −0.825; p < 0.001; n: 101).
The PCA showed different results for each horizon. The first three PCs in SH and the first two in SSH accounted for 85% and 69.5% of the variation respectively (Table 1). In SH, most parameters were clearly assigned to a single PC: PC1, which accounts for 46.45% of the variation, included LOI (exclusive for this PC), P, BasS and negative Al:Ca. Fe was mainly reflected in PC2 (20.78%), and pH was exclusively in PC3 (17.83%). eCEC and C:N were not clearly assigned to a single PC. In SSH, the two extracted PCs showed less clear relation to the parameters that contributed to both PCs. LOI was tentatively assigned to PC1 (37.08%), and pH and eCEC to PC2 (32.48%).

Soil parameters in relation to climate and herbivory at the site scale
Climatic and management factors in relation to the soil were evaluated at the site scale independently for each horizon. Soil factors were unrelated to any climatic variables and altitude for the 18 studied areas. Cattle density had a statistically significant positive correlation at SH with PC1 (Spearman rho = 0.555; p = 0.017), which mainly accounts for variation in LOI, and concentrations of available P and basic cations (Ca, Mg, Na, K), and negative Al:Ca ratio. Grazing (binary) and ponies' densities showed no correlation with any soil parameter at both depths (p > 0.1; results not shown). Soil parameters and species frequency at the site scale The seventeen selected species (seven shrubs and ten graminoids) showed markedly different frequency values at the site scale (Fig. 3a). Erica mackayana, C. vulgaris and U. gallii were nearly present in all areas but with variable frequencies. Erica umbellata and E. ciliaris were only present in four sites, although the latter had prominent cover values above 20% in three of them. In contrast, there were

Soil parameters and species frequency at the plot scale
Species frequency at the plot scale showed comparatively similar results to the site scale (Fig. 3b). E. ciliaris only occurred in three sites but was highly frequent in them, while H. marginata depicted an opposite trend with low frequencies at many sites (15 out of 18).
At the plot scale each parameter and species were analyzed independently. Erica mackayana occurred in all plots but one, therefore was not included in the analyses. Agrostis canina/ capillaris were excluded due to identification uncertainties (see methods). Calluna vulgaris and E. umbellata showed no statistically significant differences for any soil parameter. All other species showed differences for one or more soil parameters (Table 2).

Discussion
The Erica mackayana wet heathland soils in the Atlantic context The Erica mackayana wet heaths occur under highly organic and acidic soils, as has been shown for other heathlands across the Atlantic region (Field et al. 2017). We found a strong variation among sites, samples and soil depths for several parameters, for example soil organic matter (SOM) varied from c. 6% to 90%, available phosphorus (P) from 2 to 210 mg kg −1 , and pH from 3.4 to 5.6 (Fig. 2). However, our results showed that the E. mackayana wet heaths follow the general patterns described for heathland plant-soil systems for the C:N ratio values in SOM (Certini et al. 2015;Richardson and Simpson 2011) and aluminium (Al) saturation (De Graaf et al. 1997). Considering the Al:Ca ratio values obtained, the process of Al saturation of the effective cation exchange complex (eCEC) was constrained by soil depth, and correlated with SOM. A more intense recycling kinetics seem to be favored in the topsoil, where the dominant cation of eCEC is Ca, while in deeper layers the Al saturation values of eCEC can be higher than the phytotoxic threshold, affecting plant development. Similarly, P maintains a positive correlation with SOM, and in particular with the C:N ratio, which is associated to microbial mineralization-immobilization balances the promoters of ion exchange at the rhizosphere level (Richardson and Simpson 2011). According to the differences found between soil depths (Fig. 2), and positive correlation with SOM, a P recycling mechanism seems to be operating in the heathland topsoil.

Soil parameters in relation to climate and management
The Erica mackayana wet heathlands occurs in the northern most of the Iberian Atlantic biogeographic region, restricted to a narrow Cantabrian belt characterised by a mild hyper-oceanic climate (Díaz-Varela et al. 2018;Fagúndez 2016;Loidi et al. 2010;Ojeda et al. 1998). Our results showed that the community is less constrained by soil condition and can occur in a wide range of soil types, depicted by the high variation in SOM and nutrients, within its geographic distribution area. Moreover, we did not find a correlation between climatic factors (total and summer precipitation, mean and summer temperature) or altitude, and soil parameters. Other environmental factors like bed rock type can influence soils and thus be responsible for part of this variation (García-Arrese et al. 2009) but this has not been investigated here.
Grazing is the main land use for wet heathlands in the region, mainly by wild ponies and cattle and rarely by sheep (Fagúndez 2016;Muñoz-Barcia et al. 2019). Presence of ponies, but not density, have a positive correlation with heterogeneity in vegetation structure and presence and frequency of rare herbs with conservation value (Fagúndez 2016). In turn, when cattle density increases, rare species diversity is negatively affected. We found that cattle density was positively correlated with the main gradient of SOM and soil nutrients in SH. Soil nutrients in heavily grazed areas may rise due to direct deposition of cattle dung and urine (Kohler et al. 2004), but other sources of nutrient supply such as manure addition by local farmers to improve pastures cannot be ruled out. Selective grazing pressure of domestic large herbivores can promote heterogeneity by causing a net nutrient input in richer and a net loss in poorer communities (Augustine et al. 2003), and this can be the case in the study area, where cattle are often moved to improved pastures in the winter. In turn, wild ponies seem to have a neutral effect on the soil, but positive on plant diversity, adding to the evidence of a positive influence of ponies in the habitat and low rates of disturbance from their performance (Fagúndez 2016;González-Hernández et al. 2020). The only area grazed by sheep had similar cover values as those grazed by ponies or cattle compared to ungrazed areas. Sheep may have a positive impact on the community, and has been recommended for heathland conservation in other areas (e.g. Gallet and Roze 2001).

Searching for indicators: Species individual response to soil conditions
The floristic composition of the E. mackayana wet heaths includes some widespread species such as Calluna vulgaris, Deschampsia flexuosa and Molinia caerulea, but also Iberian endemic species such as Erica umbellata. Among these, M. caerulea, Erica cinerea and Ulex gallii are the species more strongly linked to soil conditions at both scales (Fig. 4). At the site scale, M. caerulea cover correlated positively with PC1, i.e. higher SOM and available P, and lower Al:Ca ratio. At the plot scale, statistically significant higher P and basic cations saturation (BasS) values were found in plots with M. caerulea. Results for M. caerulea are in accordance with other studies from western Atlantic Europe, as the species performs better in organic than in mineral soils, responds positively to N and P addition, and to Ca and Mg availability in less acid soils (Hayati and Proctor 1990;Rowland et al. 2000;Loach 1968;Roem and Berendse 2000;Taylor et al. 2001;Friedrich et al. 2011;Riesch et al. 2018). The species has increased its dominance at a broad scale (Chambers et al. 1999) and different studies have claimed that high loads of N deposition that cause eutrophication of the soil promotes a switch from heathland to grassland, because species like M. caerulea out-compete C. vulgaris and other heathers (Alonso et al. 2001;Friedrich et al.

3
Vol.: (0123456789) 2011). In the studied heaths, the species is frequent (Fig. 3) but there is no evidence of competition with shrubs (Fagúndez 2018). Molinia caerulea is thus a faithful indicator of relatively high SOM and nutrient levels, including available P, in the E. mackayana wet heaths at all scales, and general variations in its cover should be evaluated as a potential response to eutrophication.
Among the shrub species, Erica cinerea and Ulex gallii showed the highest correlation with soil conditions at the site scale, and statistically significant differences of occurrences at the plot scale for most parameters (Fig. 4). Frequency of E. cinerea correlated negatively with SOM and levels of P and other nutrients, and high Al:Ca ratio at the site scale. The species is dominant in dry heaths, but wet heath is generally considered an accessory habitat. E. cinerea generally occurs in drier and more mineral soils (Bannister 1965;Davies 1984), therefore our results could depict a limitation due to strong ecological filters for its occurrence and abundance. We expected a similar trend for E. umbellata, a species with high frequencies in dry sub-mediterranean heaths, high slopes and shallow soils Díaz-Vizcaíno et al. 1989;García-Arrese et al. 2009), but there was no evidence of any relation to soil parameters for the species. Occurrence and abundance of E. umbellata is therefore not related to particular drier conditions  (1) and absence (0) values of Erica cinerea, Ulex gallii and Molinia caerulea at two soil depths (SH: surface horizon, SSH: sub-surface horizon) for soil organic matter (represented by lost-on-ignition in %), available phosphorus (P) and the aluminium (Al) -calcium (Ca) ratio at the plot scale. Statistically significant (Man-W-Wilcoxon pair test, p < 0.05) differences in plots with presence of species positive in blue, negative in red. Light grey means non-significant differences. The box-plot represents median, 25-75% interval (box), 5-95% (bars). The last column represents correlation between PC1 of SH and frequency (0-1) of each species at the site scale, adjusted lineal model and 95% confidence band in the E. mackayana wet heath, but can be a result of opportunistic colonization through seed dispersal. Seeds of E. umbellata are much lighter than those of E. cinerea (Fagúndez et al. 2010), an advantage for effective dispersal and gap-filling at longer distances.
Gorse (Ulex gallii) is a key species of the E. mackayana wet heaths (Fagúndez 2016;Aldezabal et al. 2013;Silva-Pando et al. 2002). Encroachment of the species denotes a low habitat quality linked to lower species diversity, constrains pasture availability for livestock and increases wild fires risk (Muñoz- Barcia et al. 2019;Gómez et al. 2019). Similarly to E. cinerea, U. gallii correlated negatively with SOM, P, BasS and eCEC, and Al:Ca ratio. However, in contrast to E. cinerea, U. gallii had a negative correlation with Fe and C:N ratio at the deeper layer (SSH), therefore in the root absorption zone. Gorse is a vigorous N-fixer in which nitrogen uptake is mediated by root association with bacteria. Total N pool and nitrification rates can increase under gorse compared to heathers like Erica umbellata (Rodríguez et al. 2007). This can improve the heterogeneity of soil nutrients and diversify the community, but gorse encroachment can indicate an increase in nitrogen related to eutrophication, which poses a threat to the habitat.
High soil acidity is a major constrain for plant development, but it may enhance the ericoid species that perform well in very acid soils thanks to mycorrhizal associations (Mitchell and Gibson 2006). Ericoid species release leaves with acidic organic matter with high concentrations of phenolics and C:N ratio. Thus, pH in heathlands is mediated by soil, but with a positive return from the dominant plant species (Eskelinen et al. 2009). Similarly, ericoid species capture low concentration metals through mycorrhiza Mitchell and Gibson 2006;Read 1996), and this may explain the positive correlation with Fe for heathers (Erica cinerea, E. ciliaris) at the deep root layer of the plot scale, and C. vulgaris cover at the site scale.
Deschampsia flexuosa showed a negative correlation with pH at SSH at the plot scale, but this variation was not observed at the site scale. According to the literature, D. flexuosa occurs in highly acid soils (Scurfield 1954;Rowland et al. 2000;Van Den Berg et al. 2005), and grows better on peaty substrates (Britton et al. 2003). The species shows null or weak reaction to changes in Ca, Mg, Fe or P (Hackett 1965), and is tolerant to high Al concentrations (Pegtel 1987). As M. caerulea, D. flexuosa can be a strong competitor with heathers, but we have no evidence of that in this particular habitat (Fagúndez 2018). Other graminoids like C. binervis, H. marginata and D. decumbens occurred in less acid soils (higher pH values) at the plot scale, but this trend was not observed for species cover at the site scale. These species occupy gaps among shrubs at the fine scale (Fagúndez 2018), where the soil could be less affected by acidification from the ericoid necromass, which is highly acid (Certini et al. 2015).

Conclusion
Conservation of Atlantic heathland remnants requires maintaining the ecological filters that prevent encroachment and dominance of Molinia caerulea or Ulex gallii. Like other heathland types across western Europe, the E. mackayana wet heaths rely on the ecological filters for plant development related to soil condition. The process of net reduction of eCEC, the increase of Al saturation and the loss of basic cations could favor the exclusion of more nutritionally demanding species in these heathlands. We found no evidence of a major change in the community due to increasing nutrient availability as has been described for northern countries (e.g. Alonso et al. 2001;Haerdtle et al. 2006), but high pressure from farming, particularly higher cattle livestock densities, may affect soils and impoverish its floristic composition (Fagúndez 2016). A positive correlation between cattle density and the main gradient of nutrient availability and SOM can be interpreted as a result of the impact of these herbivores through faeces or conditional management. Conversely, wild ponies control gorse encroachment, enhances micro-scale heterogeneous structure, and maintains high plant diversity values, without affecting soil condition. Conservation measures should consider the management-soil-vegetation system to avoid biodiversity loss and associated ecosystem services. We claim that, in this particular heathland type, wild ponies best play the role of large herbivores through biomass recycling and selective diet, and that there is a potential loss of ecosystem services if ponies are removed. The wild ponies' population of the Galician mountains is declining, therefore active measures for preserving them should be applied (Fagúndez 2016;Fagúndez et al. 2021; 1 3 Vol.: (0123456789) González-Hernández et al. 2020;López-Bao et al. 2013).
Molinia caerulea, E. cinerea and U. gallii were the three species with the strongest correlation with soil parameters such as SOM, P and interchangeable basic cations (Fig. 4). Molinia caerulea occurs under higher nutrient conditions and SOM at the plot scale, and increases its frequency at the site scale. In contrast, E. cinerea negatively correlates with those parameters at both scales. U. gallii followed a similar pattern to E. cinerea but it showed lower values of C:N at deeper layers probably due to N enrichment by atmospheric N-fixing. Conversely, N inputs through fixation reduces the advantage of heathers associated to mycorrhizal fungi in such nutrient-limited soils (Read 1996). Deschampsia flexuosa occurred in highly acid soils, and other graminoids showed the opposite, although this trend was only observed at the plot scale. This set of species can be considered a first reference list of indicators for evaluating soil conditions through floristic assessment in the E. mackayana heathlands. In addition, it adds to a better understanding of some widespread dominant plants like M. caerulea and D. flexuosa, which respond similarly to northern latitudes to soil condition.
Funding Open Access funding provided thanks to the CRUE-CSIC agreement with Springer Nature.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.