Biogeodiversity and pedodiversity islands in arid lands of Europe (Almería Province, Spain)

Plant and soil landscapes across bioclimatic belts and drainage basins were studied using georeferenced databases in arid lands of SE Spain, the driest area of Europe. The syntaxonomic system was used to analyze phytocenoses and bioclimatic belts, as well as the concept of potential natural vegetation (PNV), a common approach in many countries of continental Europe. Soil types included in pedological databases were classified using the World Reference Base for Soil Resources international system (FAO 1998). Both bioclimatic belts and drainage basins effectively discriminate soil and plant assemblages in the study area of the Almeria province. The syntaxonomic perspective permits distinguishing between PNV dependent on (i) climate (climatophylous), (ii) climate and lithology, and (iii) soils (edaphophylous). Richness-area relationships of plant and soil assemblages fit well to power law distributions, showing few idiosyncratic differences. PNV, lithological associations, and soil richness are clearly correlated with the area of each climatic belt and watershed. PNV and pedotaxa richness (understood as a number of taxa at a given hierarchical level) increases from the mountain tops to the coastal lands. Around 59% of the PNV units are edaphophylous and 87% of these are edaphohygrophylous that require water supply or tolerate water excess in riverbed ramblas (dry watercourses). Edaphohygrophylous PNV are distributed in small patches within a very arid matrix. They can be considered as plant “biodiversity islands”, a concept different from that of “fertility islands” used by ecologists in arid land studies. The spatial dispersion of these phytocenoses prevents adequate preservation in the frame of conservation biology policies. At landscape level, the extent of plant communities is as follows: PNV climate dependent > PNV climate-lithology dependent > PNV soil dependent. The diversity of plant communities follows an opposite trend: PNV soil dependent > PNV climate-lithology dependent > PNV climate dependent. The PNV most conditioned by soil properties are located along the streambeds of ramblas. These fluvial sediments are not reported as soil materials in soil maps. PNV, soils and lithological associations by drainage basins conform to the predictions of the statistical tool termed nested subsets theory. However, lithological associations by climatic belts depart from this spatial pattern.


Introduction
The province of Almería houses the driest landscapes in Europe, including desert areas (Puigdefábregas and Mendizábal 1998). Many studies on soil erosion and desertification have been carried out in the SE of Spain (Brandt and Thornes 1996;Cantón et al. 2011;Barbero-Sierra et al. 2015). Less attention has been paid to the analysis of the physical geography of this region, including lithology, hydrology, soils, climate, and vegetation. A geographical study is a basic tool that should be considered before addressing more detailed research on specific aspects of structures, processes, and landscape history. In fact, whether the province of Almería and southern Murcia have suffered and/or are presently suffering intense desertification (Martínez-Fernández and Esteve 2005), or if they are already in a metastable equilibrium with the nowadays environmental conditions (Ibáñez et al. 2015) remains a matter of debate. The purpose of this study is to analyse the spatial patterns of vegetation, soils, climate (bioclimatic units), lithology, and the relations among them, implementing quantitative tools used in biodiversity and pedodiversity studies, in such a way that the results obtained can be compared with those detected in other geographical areas and environments (Ibáñez and Effland 2011). In any case, this study mainly pays attention to the understanding of soilscapes, considering several of their forming factors including vegetation, lithology, climate, as well as surrogate indicators of landforms and human impacts. Surrogate indicators of physiography, such as drainage basins (DB), are used due to lack of maps and databases on landforms for the whole of the study area. This paper does not intend to discuss processes of erosion and desertification, but rather to provide a territorial analysis that may be useful for more specific purposes, such as land planning policy as well as biodiversity and pedodiversity preservation (Magurran 2004;Ibáñez et al. 1990Ibáñez et al. , 1995Ibáñez et al. , 2012. For this purpose, previous georeferenced inventories (map layers) and scientific literature, published mainly in Spanish, were used. Likewise, different mathematical tools were applied to establish a quantitative soil and plant geography ).

Brief description of the study area
The study area covers the greastest part of the Almería Province in the south-east of the Iberian Peninsula (Figure 1). The area of Almería province is 8 769.05 km 2 (used for the bioclimatic belts analysis), while the study area covered by the drainage basins is 7 726.75 km 2 , as is explained in section 2.2.1. The tables and text showed in this paper are for the whole province unless otherwise specified.
The physiography of Almería consists of several very steep mountain ranges, rising from sea level to altitudes around 2 600 m asl in the Sierra Nevada, over a distance as short as 30-40 km (Simón 2005;Aguilar-Ruiz et al. 2004) following E-WS and SW-NE orientations. Lithology is very diverse with rocks ranging from Palaeozoic to Holocene (Villalobos-Megía 2003).
As Almería is located in the eastern end of Andalusia, its rainfall regime is isolated from the Atlantic moisture (Aguilar-Ruiz et al. 2004). Scarce rainfall, coming from the eastern Mediterranean circulation, occurs mainly in autumn. Huge evaporation causes water deficits throughout most of the year, except in the highest inland mountain ranges. The climate becomes typically xeric Mediterranean above 800-1 200 m amsl (Aguilar-Ruiz et al. 2004;Simón 2005). Low rainfall and high temperatures produce scarce and open vegetation cover that favours soil erosion during storm events.
The river channels are called "ramblas" that are similar to the ephemeral watercourses called "wadi" in Africa, "arroyos" in the western USA and "dry rivers" in other countries (Aliat et al. 2016) and their ecosystem dynamics are not well understood (Marchamalo et al. 2016;Dahlberg 2000). Watercourses are generally dry throughout the year, often for several consecutive years or decades. Locally, "ramblas" are used as arable land, roads, and for housing.
Intense stormy periods occur once every few years to decades (there are no gauging stations). These events produce flash floods typical of desert and arid environments, which cause severe damage to industrial, urban and agricultural infrastructures and erode river landforms (Villalobos-Megía 2003;Alonso-Sarriá et al. 2016;Ibáñez et al. 2015). However, flash floods also recharge aquifers along permeable streambeds and feed local groundwater that is exploited for irrigation (Shentsis and Rosenthal 2003;Lorenzo et al. 2007;Rodríguez-Lloveras et al. 2015).
Coastal lands are characterized by the presence of sand beaches, dunes, and wetlands such as saltmarshes, marine terraces, as well as small ramblas (Villalobos-Megía 2003). Badlands, alluvial fans, and endorheic basins with playas and ramblas are frequent in arid inland depressions and on gently sloping surfaces (Cooker and Warren 1973;Ollier 1976;Villalobos-Megía 2003). Typical disperse vegetation patterns termed "fertility islands" appear everywhere in the region, as also occur worldwide in arid environments (Cambardella et al. 1994;Schlesinger et al. 1996;Klausmeier 1999;Schlesinger and Pilmanis 1998;Pugnaire et al. 2004). Bare soils surrounding the fertility islands are often covered by biocrusts that increase water infiltration and soil moisture available to plants and soil fertility (Chamizo et al. 2016). Near the coast, plant communities are shrubs and bushes typical of arid/desertic regions (Brandt and Thornes 1996;Mairota et al. 1997;Puigdefábregas and Mendizábal 1998;Martínez-Fernández and Esteve 2005). Aridity results from climatic and physiographic conditions, together with the human impact since at least the last 6 000 ± 1 000 years, although Neolithic settlements could be traced back to 7 500 years (Simón 2005;Carrión et al. 2007;Cortés-Sánchez et al. 2012;García Alix et al. 2013).
Few spots of Mediterranean forests have survived long-standing land use above midaltitudes (Aguilar-Ruiz et al. 2004). At mountain tops, plant communities receive higher amounts of precipitation, favouring humid vegetation of pastures, meadows and shrubs (Blanca et al. 1998).
In arid and semiarid areas, climate and topography are the driving forces of soil formation (Ollier 1976), whereas many soil properties are inherited from parent materials (Buol 1965). Young soils such as Regosols and Leptosols are mainly located on slopes exposed to continuous erosion (Simón 2005), whereas more evolved soils including Calcisols, Solonchaks and Gypsisols are frequent on the low-lying arid flatlands (Cooker and Warren 1973;Aguilar-Ruiz et al. 2004).
Soil water and temperature regimes are dominantly aridic and thermic along the coast and in inland endorheic basins, while they are xeric and mesic above 800-1 000 m elevation (Aguilar-Ruiz et al. 2004).
Although water culture in the SE of Spain has been of paramount importance for millennia, traced back to prehistoric times, it was only after the Arab invasion of the Iberian Peninsula that water management methods were generalized through a variety of techniques. Much of the agricultural crops and pastures in the arid areas of Almería have depended on these ancient hydraulic technologies (Muñoz-Muñoz 2000).

Materials and methods
Description and statistical analysis of the Almería landscapes are based on previous georeferenced inventories (map layers) using Geographical Information Systems. Lithology, soil and vegetation layers were analysed by overlaying the georeferenced datasets by (i) drainage basins (DB) and (ii) bioclimatic belts (BB). Spatial pattern analysis was carried out along two perpendicular transects: DB land segmentation occurs perpendicular to the coastal line and BB segmentation divides the study area into altitudinal belts more or less at right angles to DB (Figure 2). The smallest DB of Hortonian rank 1 starting at low elevations drain to the sea and are usually dry.  Ibáñez et al. (2015, 2016a) carried out different analyses of soil assemblages of the Almería province by DB but not by BB using the same mathematical formalism and pedological classification.

Data layers
Bioclimatic belts (BB) are at a scale of 1:400,000 (CMAYOT 2004), and potential natural vegetation units (PNV) at a scale 1:10,000 (CMAYOT 2005). Georeferenced datasets are freely available from the spatial data infrastructure REDIAM, and the lithological map at a scale 1:400,000 is available in IGME-Junta de Andalucía (2011). PNV maps were established using the so-called syntaxonomic system (Mirkin 1989;Rivas-Martinez 2005). This approach, based on the International Code of Phytosociological Nomenclature (CPN) (Weber et al. 2000), permits the identification of landscape units in which plant communities are subject to similar environmental conditions (climatic, pedologic, geomorphic and lithological features), regardless of their stage in an ecological succession (Rivas-Martínez 2005). PNV are named mixing the CPN terminology with the main environmental factor affecting them. For example, edaphophylous phytocenoses are conditioned by soil and/ or parent material properties. Suffixes can be added such as basophilous (rocks or soils rich in nutrients), acidophilous (rocks or soils poor in nutrients), edaphohygrophylous (plants exposed to permanent or seasonal waterlogging), edaphoxerophylous (plants adapted to growing on soils with restricted water supply or xericity), etc. (Rivas-Martinez 2005). This geobotanical school has been widely used in many countries of continental Europe but not in the Anglo-Saxon world. Similarly, the European Habitat Directive strongly relies on phytosociological and syntaxonomic systems to carry out the inventory of the natural habitats to be preserved at continental scale. This approach is also used in many countries for instance in territorial planning and restoration ecology (Loidi and Fernández-González 2012).
Lithological units of the geological map were identified following the norms established by the IGME, the official Geological Service of Spain. Geological maps have many shortcomings. There is no universal classification accepted worldwide by geologists and geomorphologists. Likewise, geological maps usually do not show the distribution of all rocks exposed, as would be desirable to analyse the relationships between lithology, soils and vegetation.
The soil dataset used was the digital soil map of Almería at the scale of 1:100 000 (Aguilar- Ruiz et al. 2004). Soil types were described using the World Reference Base for Soil Resources or WRB classification (FAO 1998), and mapped at the scale of 1:50 000. The legend of the soil map provides information on the surface area covered by the soil types included in each soil association. FAO rules were used for the purpose of segregating associated soil types.
The DB basins layer was created from Terrain Digital Models of the Spanish National Geographic Institute (IGN 2009). This layer was matched with the Almería province boundary polygon to get the DB used in this study. Each DB was assigned a rank according to the Horton-Strahler method (Strahler 1957) which represents its stream order or branching complexity. The final layer consisted of 19 DB, with ranks between 1 and 4, plus two small portions of the Guadalquivir (rank 6) and Segura (rank 5) rivers which were not considered relevant for the previous analysis of DB (Ibáñez et al. 2016a). Soils and PNV were overlaid with the DB and BB georeferenced layers. Therefore the figures provided by geographic information systems differ between DB analysis and BB (the whole of the Almería Province).

Data processing
Statistical analysis was conducted on the digitized soil map using GIS tools. All data were introduced in a spatial geodatabase using ArcGIS 10.1 software. Afterwards, the polygons of each lithotype, PNV type and soil association were exported into Excel spreadsheets for statistical analyses, computing taxa richness and Shannon's index values, linear regression models, taxa-area relationships, and a nested subset analysis. Because taxa diversity richness increases with the area, a nested subset analysis informs if PNV, soils and lithological units are distributed randomly or show an idiosyncratic spatial pattern by DB and BB. The objective of the nested subset analysis is to determine if taxa present in the DB or BB with smaller areas are also present in larger ones but not the reverse. Larger DB and BB host singular taxa that do not appear in those that have less extent in the territory studied, as biodiversity and pedodiversity studies predict (Ibáñez et al. 2005b). The nested subset theory was proposed by Patterson and Atmar (1986).
After this process, BB and DB were classified with their numerical lithotypes (only for DB analysis), plant and soil diversities. Data were sorted by decreasing number or spatial extent of the polygons, obtaining the so-termed rankabundance plots (RAP). The fragmentation index was also calculated, as the inverse value of the average polygon size.
Lithology, PNV, soils, and BB datasets were analyzed using statistical procedures common in biodiversity and pedodiversity studies that are well documented in the scientific literature (Magurran 2004;Ibáñez et al. 1995;Ibáñez and Bockheim 2013).
Results were represented on single RAP, a representation preferred by experts to display relative taxa abundance in which taxa are ranged on the X axis according to their decreasing abundance, while their relative abundance is shown on the Y axis. Richness is the number of taxa present in a given assemblage. The popular Shannon Entropy index was used to estimate the diversity. Richness-area relationships (SAR) were established to test if datasets conformed to a power law as is usual in biodiversity and pedodiversity analysis. The mathematical expression of the Shannon index is as follows: ( 1) where p i is estimated by means of n i /N (n i is the area covered by the i-th taxa and N is the total area under study).
The SAR is estimated with a power law using the following formula: where S is the number of taxa encountered, A is the area estimated of a given geographical unit (e.g. a BB), C is an empirical constant, and z is the slope of the log of the number of taxa in the geographical unit (GU) considered against the log of the GU area. Thus, C and z are empirical constants which would vary from one taxon to another and from a GU to another.
Pedological and biological assemblages have been usually fitted to one of the following abundance distribution models (ADM): geometric series, logarithmic series, lognormal distributions and the broken stick models (Magurran 2004). However, more recently power laws have also been proposed to describe ADM, as well as a novel diversity index (Pueyo 2006) that was used in this paper to fit taxa assemblages.
Conditions for the development of a nested subset structure (Patterson and Brown 1991) include: (i) compared taxa and sites must have a common geographic history and similar contemporary ecological conditions; and (ii) the possibility of the existence of some kind of hierarchical pattern among taxa or sites. However, this statistical tool has raised some controversy and debates, mainly on the methods to quantify nestedness as well as the best null hypothesis (Ulrich et al. 2009). In this study, we make use of the classical approach and software termed "Temperature Calculator" Patterson 1993, 1995) in spite of its potential limitations. The purpose is to detect if PNV and soil taxa follow similar mathematical regularities. Nestedness is a property of assemblages, not of individual taxa, so it is reasonable to ask whether nestedness reveals anything that cannot be determined from the ordered distributions of constituent taxa (Simberloff and Martin 1991). This mathematical tool can be applied to datasets of soils and PNV, but unfortunately not to the dataset provided by the geological map that lacks the level of information required by this type of analysis.
Database scales, number of polygons and taxa of the databases are showed in Table 1.

Regional and ad hoc soil classifications
Lithological, vegetation, and soil taxa were classified according to their respective presence in DB of different hortonian ranks (unpublished data) and BB respectively. This ad hoc classification in eight classes followed the guidelines established by Ibáñez et al. ( , 2015 for the whole of Europe. The classes are the following: (1) Dominant taxa are those that are absent in less than two DB or one BB respectively; (2) Abundant taxa are those that occur in more than five DB or two BB respectively; (3) Common taxa are those that occur in more than half of the DB or three BB belts respectively; (4) Frequent taxa are those that only appear in three BB (only in BB analysis); (5) Rare taxa or taxa minorities are those that only appear in two DB or BB respectively; (6) Endemic or unique taxa are those that only occur in a single DB or BB respectively; (7) Idiosyncratic taxa are those that are only present in the DB of greater extent or hortonian rank (only in DB basins analysis). These are excluded of the Table 3 and 4 but specified in their respective legends; (8) Other taxa are those that cannot be included in the preceding categories (only in DB basins analysis) PNV units and pedotaxa were grouped for certain purposes in "functional" ad hoc groups to describe the structure of plant and soil landscapes as well as to avoid the use of complex and unfamiliar syntaxonomic nomenclature and facilitate the understanding of the WRB soil types.

Results and Discussion
3.1. A geographical overview of soils, vegetation and lithology using regional classification.
Tables 3, 4, and 5 show PNV units and pedotaxa according to their respective geographical classification as described in section 2.2.3. Table 3 shows soil associations of lithotypes by DB that were not used in the BB analysis. In contrast, Table 4 and Table 5 show the classification of PNV and pedotaxa by BB, but     Tables 4 and 5, are higher than for lithologies ( Table 3). This shows that fluvial landforms in large watersheds (including complex floodplains, fluvial terraces, marshes) and other PNV common of rivers that carry more water and sediments induce the appearance of idiosyncratic soils and vegetation in contrast to small DB with low sediment load and water discharge, as well as with poorly developed fluvial landforms. The larger the number of PNV units, the smaller their area of distribution in terms of BB (Table  5). Thus, the highest PNV richness corresponds to the endemic class, indicating that each BB is characterized by several endemic PNV units. Endemic and minority PNV richness is higher than soil type richness. This partially results from the large areas covered by shallow soils (Regosols and Leptosols) as is show in Section 3.2.  (Ibáñez et al. 2015). Figure 4 show the rank-abundance plots (RAP) of Reference Soil Groups and soil types at the second level of the WRB (FAO 1998) in the Almería province. As is common in this type of analysis, few taxa or units are very abundant whereas most of them are scarce (Magurran 2004   Legend: Codes suggested for WRB Reference soil Groups and qualifiers (FAO 1998  of the most populated and disturbed coastland areas. It has been shown in previous studies that human impact causes at the beginning a decrease of pedodiversity values with extinguishing natural soils, and after that manmade pedotaxa begin to appear (Lo Papa et al. 2011;Zhang 2013).

Figure 3 and
The study area covers only a small area of the Oro-Mediterranean belt, typical of the mountain tops, that extends to the west in the neighboring Granada province. For this reason, several pedotaxa typical of the high Mediterranean mountain ranges of Andalusia, rich in organic matter and poor in nutrients (Haplic Umbrisols and Leptic Umbrisols), and sometimes waterlogged (Gleyic Umbrisols, Gleyic Phaeozems, Gleysols, and Histosols) are not included.
Pedotaxa can be grouped in an "ad hoc" classification according to the soil properties that strongly influence the type of vegetation that can grow on them, as it is shown in Table 6. We term these clusters as "functional/ecological pedotaxa".   lowland soils are replaced by pedotaxa where salts are leached to the lower soil horizons as the climate becomes wetter (e.g. Cambisols and Luvisols). Well-developed Luvisols, common in mountainous areas, require stable landforms and/or a dense and permanent plant cover; consequently they are very scarce in the study area that is largely deforested and strongly eroded. Thus the Almería landscape has a sort of altitudinal catena structure. The Low-Meso-Mediterranean BB seems therefore to be a proxy for the delimitation of the arid lands.
Aridic Solonchaks occur in the most arid areas, usually in sites where also subsaline marls appear, but this parent material is not the only determinant of their presence/absence. Gleyic Solonchaks are restricted to the shoreline where saline and salt marshes are also common. Gypsisols occur in the most arid conditions of the Thermo and Meso-Mediterranean BB, but not in gypsum deposits. The area of Sorbas is known for having important gypsum outcrops. But highly soluble gypsum is easily exported by water, leaving impressive endokarstic features. Thus Aridic and Hypergypsic Gypsisols are not representative of the gypsum formations, where Gypsiric Regosols are common instead. Table 5 shows that only about 20% of the edaphophylous PNV group is edaphoxerophylous, growing on soils with a more arid pedoclimate than that of the surrounding areas (e.g. on gypsum or parent material rich in sodium chloride, also on Leptosols and associated rock outcrop cracks). The remaining 80% are edaphohygrophylous, appearing in many sites irrespective of other environmental conditions. Edaphoxerophylous PNV are more or less homogeneously distributed in most BB, while edaphohygrophylous PNV are more abundant in arid territories at medium and lower altitudes, in small physiographic patches that favor the storage of water in the solum. They appear as scattered spots in the surrounding soilscape matrix. Most of the edaphohygrophylous PNV are endemic or are soil minorities, are frequently associated with larger basins (idiosyncratic pedotaxa) and have no clear relationship with Fluvisols. This spatial pattern resembles that of the fertility islands identified by various authors (Schlesinger and Pilmanis 1998) in arid environments such as Almería, at more detailed scales (Pugnaire et al. 2004). Thus, about 40% of PNV diversity in Almería is located in small geoforms subject to permanent or seasonal waterlogging without correspondence with the typical soils of fluvial landforms. Table 7 shows a summary of the PNV and soil richness as well as their respective Shannon Index values by the BB belts of the Almería province.

Figures 5 show taxa-area curves for PNV (a)
and soil types in the study area (b). In both cases, the fit to a power law distribution is clear. This pattern also occurs when PNV and soilscapes are analyzed by DB (Ibáñez et al. 2015) in agreement with most pedodiversity and biodiversity literature. However, the relationship between richness and area of BB belts is better for PNV than for pedotaxa. The pedotaxa number increases from mountain tops to coastlands. Taking into account the area of each BB (Table 7) , it could be suggested that factors other than the area are needed to explain the high richness of the coastal soilscapes, such as the abundance of fluvial landforms and manmade soils. In contrast, PNV-area relationships are clearer. In both cases diversity index and richness values follow the same trend.  [ IBÁÑEZ J. J., PÉREZ-GÓMEZ R., OYONARTE C. & ZINCK A. ] In fact most of the variables in the study area fit well to power laws even when soil associations (SMU) are used instead of soil types: (i) SMU area-PNV richness; (ii) SMU area -number of climatic belts; (iii) Soil types area -SMU richness; (iv) PNV area -SMU richness; (v) PNV area-soil types richness; (vi) PNV area -number of climatic belts; (vii) Climatic belts area -SMU richness; (viii) Climatic belts areasoil type richness; (ix) Climatic belts area-PNV richness, being the exception SMU area -soil type richness that does not conform to a power law. The statistical significance of the fit was tested using the classical R 2 (Feoli et al. 2015).
The results of the nested subset analysis of soils and PNV by BB do not show a clear nested structure, in contrast to the one obtained with the same datasets by DB (Ibañez et al. 2016a) or in previous studies carried out using archipelagoes with different island sizes or fragmenting continental landscapes by DB of different sizes (Ibáñez et al. 2005a,b) or the whole of the pedosphere by administrative boundaries  (Ibáñez and Feoli 2013). PNV and pedotaxa of small DB are nested subsets of the larger ones with a very high degree of statistical significance (Ibáñez et al. 2016a). Likewise, the idiosyncratic pedotaxa of the larger DB are effected by (i) the capture of larger areas from the shoreline to the mountain tops (and thus including many different BB) and (ii) the presence of more developed fluvial landforms and sediments with idiosyncratic soils and associated PNV, as was described in the geographical classification in section 2.2.2. Therefore, different perspectives lead to the detection of distinct but complementary landscape patterns. In other words, spatial patterns identified for PNV and soil depend on the units resulting from the land fragmentation (DB or BB in this paper). The reason is that, in contrast to the DB segmentation, BB violates one of the conditions necessary to apply this type of statistical analysis: the comparison requires that taxa and sites have similar environmental conditions, and that does not occur in the analysis by bioclimatic belt (Patterson and Brown 1991). It could be concluded that there is no universal pattern independent of landscape fragmentation decisions used by experts.
3.5. Soil-plant community relations Table 8 shows a general picture of the PNV assemblages in the study area according to the criteria of the geographical classification proposed in section 2.2.3. To provide more complete information on plant landscape richness and diversity the data should be supplemented by their respective coverages, as is shown in Table 9.
The study area is mainly covered by climatophylous and lithological-climatically dependent PNV, in contrast to edaphophylous PNV. There are no dominant PNV (Table 8), while Table 9 shows that lithological-climatic and climatophylous PNV follow opposite trends. Lithological dependent PNV units are more abundant at medium and high mountain ranges; meanwhile climatophylous PNV are mainly in lowlands. Edaphophylous PNV are clearly richer and more diverse ( Table 7), but are geographically mainly restricted (see endemic and minority PNV in Table 5) to the lowlands ( Table 9) and cover small areas (Table 8). This fact indicates that, in the Almería arid lowlands, most PNV diversity consists of plant communities that depend on soil properties representative of a narrow range of climatic conditions. It seems therefore that soils play a relevant role in the habitat heterogeneity of the Almería province (Ibáñez et al. 2016b). Furthermore, the fragmentation index shows that the edaphophylous PNV are highly fragmented in small islands of biodiversity within a more uniform "matrix", where climate or climate and lithology determine the whole of the landscape (    Table 9. Percentages of land covered by potential natural vegetation in the Almería province conditioned by soils, lithology and climate appear in small riverbed spots where the sediment cover allows water upwelling from underlying aquifers, as Shentsis and Rosenthal (2003) demonstrated in other arid lands. This pattern means that the spatial distribution of edaphohygrophylous PNV and associated soil types are stratigraphically/tectonically dependent as it also happened with peasant agriculture in the past (Muñoz-Muñoz 2000).
The polygon size distributions conform to a scaling law, indicating that these types of PNV are distributed in a fractal way, whether individual edaphohygrophylous PNV or the whole set thereof are considered.
Erratic flash floods following torrential rainstorms produced by specific weather conditions called "gota fría" or cold drop, result in intense erosive events, as the large coverage of Regosols and Leptosols testifies (Aguilar-Ruiz et al. 2004;Simón 2005). However, the diverse PNV (highly rich in vascular plants) and the ancient agriculture in the region depend on aquifer recharges and subsequent water upwelling in the rambla riverbeds caused by "gota fría" events. Given the low frequency of "gota fría" rainstorms, sometimes separated by decades, soils develop in the rambla streambeds, causing PNV and soil fauna to become abundant and rich in species (Zaballos 1997;Ortuño et al. 2013). Riverbed sediments and organic materials are often ephemeral as they are exposed to be washed away by repeated flooding. For this reason, they are usually reported in soil maps as miscellaneous areas instead of soil proper. Such pedotaxa have no place in current soil classification schemes, and until very recently this singular habitat type was referred to neither in the ecological nor in the pedological literature (Ortuño et al. 2013). It is therefore interesting and desirable to study them in depth (Freddy Nachtergaele pers. comm.). Rambla PNV are very varied in species composition and water contents, as might be the unstable and discontinuous soil substrate.

Conclusions
The landscapes of the Almería province are typical of arid mountains with certain features belonging to the Mediterranean biome from 800-1 000 m elevation to the mountain tops. Shallow soils such as Regosols and Leptosols are widespread throughout the landscape. These poorly developed soils resulting from longstanding slope erosion hold very little water. In the arid lowlands, scarce and erratic rainfall together with high potential evapotranspiration favours salt accumulation throughout the soils. Thus the lack of water and high temperatures determine the types of soil and vegetation present in the lowlands, including Calcisols, Gypsisols and Solonchaks. Concatenation of bioclimatic belts (BB) conforms a typical arid/semiarid altitudinal catena. Both vegetation and pedotaxa richness are related with surface area following power laws, but with minor differences. The number of pedotaxa increases from highland to coast, whereas Potential Natural Vegetation (PNV) richness is more related with the area covered by each BB. The Shannon diversity index follows similar trends.
While some pedotaxa and soil landscapes are distributed throughout the study area, especially the shallow Leptosols and Regosols, each bioclimatic belt has many idiosyncratic PNV and pedotaxa that are restricted to one or two BB. Thus the results obtained for soils and PNV by BB do not conform to the nested subsets theory in contrast to what happens when the same databases are analyzed by the size of the drainage basins as it was demonstrated in a previous study. This suggests that the results obtained using different fragmentation units (e.g. watersheds, climate belts, administrative units) should not be generalized. Therefore, there are no universal landscape regularities as the conclusions obtained depend on the a priori landscape fragmentation decisions taken by experts.
The areas of the PNV that depend mainly on climate conditions (climatophylous) are the most extensive, followed by the PNV that depend on lithology and climate. In contrast, the areas of PNV related with edaphic factors (edaphophylous) are small and fragmented. Lithological-climatic dependent PNV are more abundant at higher and medium altitudes, whereas climatophylous and mainly edaphophylous PNV occur usually in arid lowlands. Most edaphophylous plant communities within arid units require the presence of soil water (edaphohygrophylous). PNV associated with the driest edaphic media (edaphoxerophylous) appear on saltrich (sodium chloride) and gypsiferous soils and sediments and/or different kinds of rock outcrops (e.g. dolomitic limestones). In view, that richness of edaphohygrophylous PNV is very high, appearing in small widespread plots mainly associated with the most arid climatic belts, where they appear as islands of plant biodiversity within a sea of arid lands.
Richness and diversity of climatophylous, lithological/climatic and edaphophylous PNV follow a pattern opposite to that of landcover. Edaphophylous and specially edaphohygrophylous PNV are the most diverse but due to their small size and to the fact that they are sometimes not recognized as soils in soil surveys, they are not reported in the soil map units surrounding watercourses. Many edaphohygrophylous PNV grow on stony or sandy riverbeds, and these geoforms are not taken into account in current classifications schemes and soil maps, although organisms living in rambla riverbeds are typical soil organisms. It is suggested that the rambla habitat, similar to the wadis of other arid and desertic regions, receive adequate attention in soil surveys. Flash floods are one of the most important driving forces of the landscapes in the Almeria province. Sporadic and intense rains erode the soil cover and generate natural disasters. However, these events also recharge the aquifers that maintain the vegetation diversity in the ramblas, and have allowed traditional agriculture. Hydraulic works to prevent floods can cause the loss of an important part of biodiversity and pedodiversity under edaphohygrophylous PNV in Almería. Pedodiversity and PNV diversity preservation, as part of the natural heritage, will be difficult and controversial due to the particular spatial distribution of the rambla habitat.