Defining the importance of landscape metrics for large branchiopod biodiversity and conservation: the case of the Iberian Peninsula and Balearic Islands

The deficiency in the distributional data of invertebrate taxa is one of the major impediments acting on the bias towards the low awareness of its conservation status. The present study sets a basic framework to understand the large branchiopods distribution in the Iberian Peninsula and Balearic Islands. Since the extensive surveys performed in the late 1980s, no more studies existed updating the information for the whole studied area. The present study fills the gap, gathering together all available information on large branchiopods distribution since 1995, and analysing the effect of human population density and several landscape characteristics on their distribution, taking into consideration different spatial scales (100 m, 1 km and 10 km). In overall, 28 large branchiopod taxa (17 anostracans, 7 notostracans and 4 spinicaudatans) are known to occur in the area. Approximately 30% of the sites hosted multiple species, with a maximum of 6 species. Significant positive co-occurring species pairs were found clustered together, forming 4 different associations of large branchiopod species. In general, species clustered in the same group showed similar responses to analysed landscape characteristics, usually showing a better fit at higher spatial scales.


Introduction
Large branchiopods are a key faunal group of crustaceans characteristic of temporary aquatic systems (with few exceptions) which can be found across all continents . Although nowadays nearly 600 species of large branchiopods are known worldwide (Rogers et al., 2015), a large part of them are known only from a few localities ([50% of the species are known from B10 localities in the case of anostracans; Rogers, 2013), indicating that the faunal group harbours a large amount of biodiversity, and that the knowledge of the group is still far from being complete (Rogers et al., 2015). This is also the case in the western Mediterranean, where descriptions of new species are still being performed (e.g. Korn et al., 2010;Boix et al., 2016) and some efforts to describe the distribution of large branchiopods are made at regional (e.g. Samraoui & Dumont, 2002;Culioli et al., 2006;Marrone & Mura, 2006; and national level (e.g. Thiéry, 1987;Defaye et al., 1998;Mura, 1999;Van den Broeck et al., 2015a;Marrone et al., 2016). The effects of environmental factors on the distribution of large branchiopods are largely known, particularly those related to habitat characteristics (e.g. salinity, turbidity, temperature, surface, depth, altitude, vegetation cover; Alonso, 1998;Boven et al., 2008;Nhiwatiwa et al., 2011;Gascón et al., 2012;Horváth et al., 2013;Sahuquillo & Miracle, 2013;Stoch et al., 2016). Moreover, some spatial variables at regional level are also known to affect species distribution, such as the closeness of sites (which can be related to dispersal or shared environmental characteristics; Nhiwatiwa et al., 2011;Horváth et al., 2013), climatic gradients  or habitat fragmentation (Gascón et al., 2012). The sensitivity of large branchiopods to all this range of environmental and spatial variables makes them especially interesting as biological indicators of the conservation status of temporary wetlands (e.g. Sahuquillo & Miracle, 2015;Van den Broeck et al., 2015b;Lumbreras et al., 2016).
However, the effects of changes in land use and the fragmentation and loss of habitats have been less explored on large branchiopods, although these factors are recognized as some of the main threats to global biodiversity. In the case of the Mediterranean region, larger biodiversity losses are expected in future scenarios, due to the sensitivity of its ecosystems to drivers of biodiversity change, especially those related to land-use change and introduction of exotic species (Sala et al., 2000;Underwood et al., 2009). Although temporary ponds were so far only minimally impacted by traditional agricultural uses, their disappearance and degradation have been highly intensified during the last 50 years, due to infrastructure and urban development, and intensification of agricultural activity (Gallego-Fernández et al., 1999;Rhazi et al., 2012). This increasing degradation and loss of habitats (Belk, 1998;Zacharias & Zamparas, 2010) strongly impact large branchiopod populations, leading to greater distances between remaining populations, eventually leading to the loss of local populations (Eder & Hödl, 2002) and even some species on the regional scale (Martens & De Moor, 1995;Mura, 1999;De Roeck et al., 2007). Therefore, the landscape structure surrounding the wetlands should be taken into account when evaluating the factors affecting large branchiopods assemblages, since some landscape characteristics are known to have great importance for biotic communities, such as landscape heterogeneity or land-use types, acting at several spatial scales (Weibull et al., 2000;Hall et al., 2004;Hartel & von Wehrden, 2013). The evaluation of these landscape characteristics for the conservation of the large branchiopod assemblages has been seldomly carried out, but the type of land use is known to affect assemblage structure of passive dispersers (Hall et al., 2004), as well as some species at different spatial scales , whereas for other species, the land use surrounding the wetland was shown to be irrelevant Horváth et al., 2013). To this end, the integration of landscape characteristics will improve large branchiopod conservation programs at regional and national level, and they will help to better understand the drivers of biodiversity change linked to future scenarios of global change.
The main goal of this article is to describe the overall biodiversity patterns of the large branchiopods present in the Iberian Peninsula and Balearic Islands, establishing a general overview for the conservation of this faunal group in this area. Additionally, we aim to determine if the biodiversity patterns detected are similar across the entire study area, and to investigate the role of land use, landscape structure, population density and protected area networks on the large branchiopod species distribution. These will set a basic framework to understand their distribution and will allow us to assess whether there are particular factors affecting this faunal group that should be taken into consideration in the evaluation of their conservation status.

Study area
The Iberian Peninsula and the Balearic Islands are situated in the western Mediterranean, and cover an area of about 587,000 km 2 . The area is shared by several countries (Spain, Portugal, Andorra, France and Gibraltar), although the scope of this study focusses only on Spain and Portugal. The climate types present in the study area are mainly arid and temperate climates, with cold climates circumscribed to the higher peaks in the mountain ranges (Peel et al., 2007;AEMET & IM, 2011). Arid climate is restricted mostly to southeastern Spain, although there are some areas in northeastern (Ebro valley), central (southern Meseta and Extremadura) and southwestern (Alentejo) Iberian Peninsula, and Balearic Islands (southern Majorca, and the islands of Ibiza and Formentera) with this type of climate. Temperate climates are divided into those with dry summer (which comprises the majority of the area of the Iberian Peninsula, from the western Atlantic coast to the Mediterranean sea, excluding the arid southeastern Spain) or without dry season (restricted to the northern Iberian Peninsula-Cantabrian coast, northern Meseta and the Iberian and Pyrenees mountain ranges).

Species distribution
All known bibliographical and unpublished large branchiopod records from the Iberian Peninsula and Balearic Islands were compiled (dating back to 1916), and geographically referenced to the maximum possible resolution whenever possible (see Online Resource 1 for bibliographical references). This allowed to create the checklist of the large branchiopod fauna, and the global distribution in the Iberian Peninsula and Balearic Islands (see Online Resource 2). However, as it is unreliable to relate old records to present-day landscape metrics, all analyses were limited to records from 1995 onwards. Moreover, some records were not considered for the study of the biodiversity patterns due to particularities in the biology or in the taxonomic resolution of different taxa. Thus, the exotic invasive Artemia franciscana Kellogg, 1906 was not considered as this species competitively displaces the autochthonous species of Artemia. Records of parthenogenetic strains of Artemia with undefined ploidy were also not considered, because parthenogenetic strains with different ploidies (i.e. diploid and tetraploid) were considered as two different taxa, as recent studies suggest that they could be related to different bisexual Artemia species Asem et al., 2016). In the case of genus Triops, we considered all species of the genus as one taxon, due to the difficulties to allocate the bibliographical records of this genus to the 6 existing species present in the Iberian Peninsula (Korn et al., 2006(Korn et al., , 2010; at this moment, syntopic occurrences of different species of Triops have not been recorded in the Iberian Peninsula or the Balearic Islands. For each species, the number of occurrences in the study area allowed to determine their rarity with rankfrequency curves following Siqueira et al. (2012). The inflection point of the curve classified the species as common (those at the left side of the curve) or rare (those at the right side of the curve; in this case those that were present in less than 50 sites; see Online Resource 3). Moreover, those species that contributed with less than 1% to the total occurrences were considered very rare (in this case, those that were present in less than 15 sites).
Biodiversity metrics were calculated using two different approaches. In one hand, for each site species richness and composition was recorded in order to describe biodiversity patterns across the whole study area. We considered co-occurring species as those species that were recorded in the same site, but not necessarily at the same time (i.e. syntopic species). However, data on synchronic species (i.e. those species that were recorded at the same time in the same site) appear in Online Resource 4. On the other hand, for each UTM 100 km grid square we calculated the number of sites with presence of large branchiopods, the mean species richness per site, the total species richness in each grid square and the percentage of sites that showed co-occurring species. Seventy-seven UTM 100 km grid squares had presence of large branchiopods in the Iberian Peninsula and Balearic Islands, but only those with more than 5 sites were included in the analyses (39 out of 77).

Landscape metrics
Landscape metrics were calculated considering human population density (DP, based on local administrative units), protection status (N2000, Natura 2000 network areas), land cover (CORINE land cover) and landscape structure. We assigned to each site the corresponding population density of the local administrative units (LAUs) where it was located. LAUs, as well as their corresponding population, were obtained from 2011 population census (European Statistical System, 2011). Natura 2000 network areas (EEA, 2015) were used as the protection category for the studied region instead of other national protected areas due to its high representativeness of biodiversity (Rosati et al., 2008;Abellán & Sánchez-Fernández, 2015). Sites that were in a Natura 2000 network area were classified as Protected. Land covers were classified in 6 categories based on CORINE Land Cover inventories (EEA, 2006): artificial surfaces (AS), forests and seminatural areas (FS), irrigated agricultural areas (IA), non-irrigated agricultural areas (nIA), wetland areas (W) and water bodies, although this last category was not used to test its effects on the biotic assemblages. Then, for each site, we calculated the relative percentage of surface occupied by each category using 3 different spatial scales (circular areas were constructed using 100 m, 1 km and 10 km radius distances to the central point of each site as centroid). We used ArcMap GIS 10.0 to process all datasets and calculate landscape metrics (ESRI, 1999). Finally, following Weibull et al. (2000), landscape diversity (Ldiv) was calculated as a measure of landscape structure using the Shannon diversity index, Ldiv = -R p i ln(p i ), where p i is the proportion of habitat i in the circular areas used for calculating the land covers. In the case of the UTM grid square approach, we calculated the mean values of the relative percentage of surface occupied by each category of land covers and landscape diversity for all the sites present at each grid square.

Species associations
Species associations were explored using two different approaches: (1) co-occurrence analysis and (2) classification analysis. The co-occurrence analysis allowed us to identify significant co-occurring species, either with positive or negative co-occurrences. Moreover, the effect size (i.e. the difference between expected and observed frequency of co-occurrence) was also obtained in order to quantify the strength of the pattern. In the second approach, a non-hierarchical cluster approach was used to identify species assemblages. Partitioning around medoids (PAM) was the classification technique chosen to cluster species into groups, since it has been pointed out as more robust than more classical approaches like k-means clustering (Borcard et al., 2011). Sørensen similarity was used as the similarity measure. The number of groups tested representing existing species assemblages was from 2 to 16 (number of large branchiopod species minus 1, only considering the species with more than 5 occurrences). Silhouette widths of the resulting groups were obtained to examine group membership's appropriateness, as well as the overall classification result using the Rousseeuw quality index. The silhouette width ranges between -1 and 1, and it is a measure of the degree of membership of an object to its cluster, negative values indicating misclassified objects. Cooccurrence analysis was performed using package ''cooccur'' (Griffith et al., 2016), and PAM analysis was performed using package ''cluster'' (Maechler et al., 2016).

Species and landscape metrics
Hierarchical partitioning was used to identify the relationship's sign (positive or negative) and the relative magnitude of the effects of the landscape metrics tested at different scales (areas of 100 m, 1 km and 10 km radius). This method assesses the independent, joint and total contribution (relative influence) of each predictor variable by averaging a measure of goodness-of-fit over all possible models that include that predictor variable. Therefore, it is less susceptible to multicollinearity problems than are the singlemodel approaches when looking for the ''best'' model (Logan, 2010). In order to evaluate whether the magnitude of the contribution of a variable is significant (Z C 1.65 at the 95% level), a randomization procedure was used in which the independent contributions of each predictor variable were compared to distributions of such contributions generated by repeated randomizations of the data matrix (number of randomizations = 100). While a binomial family error was specified for the species hierarchical partitioning and randomization procedures, a Poisson family error was used for species richness analyses. When over-or under-dispersion was detected, P values were corrected using quasi-likelihood model (Logan, 2010). The sign of the relationship between species presence-absence or richness data and landscape metrics was assessed with generalized linear models (GLM) using all environmental variables as predictors. To assess the model goodness-of-fit, we employed the squared Pearson correlation coefficient between observed and predicted values. Hierarchical partitioning and randomizations were carried out using package ''hier.part'' (Walsh & Mac Nally, 2013).

Biodiversity and landscape metrics
Conditional inference tree (CIT) models were used to assess the relationship between biodiversity (response variables: mean species richness per site, the total species richness in each UTM 100 km grid square and the percentage of sites with co-occurrences) and landscape metrics (explanatory variables: DP, N2000, AS, FS, IA, nIA, W, Ldiv) in the 3 tested spatial scales (100 m, 1 km and 10 km radius). This type of regression displays a binary tree, built by a process known as binary recursive partitioning, which gives a very clear picture of the structure of the data, and provides a highly intuitive insight into the kinds of interactions between variables (Crawley, 2002). CIT are not affected by over-fitting and are unbiased with regard to the types of explanatory variables used (Hothorn et al., 2006;Strobl et al., 2007). CIT models were developed using the ''party'' package for R (Hothorn et al., 2006). Three models, one for each biodiversity metric, were run with the different spatial scales.
All statistical analyses were run using R 3.0.1 (R Core Team, 2013).

Results
In total, 1808 bibliographical and unpublished records of large branchiopods in the Iberian Peninsula and the Balearic Islands were obtained (see Online Resource 1), corresponding to 1167 different sites (i.e. water bodies). One hundred seventy sites had only records before 1995 (Fig. 1 Korn et al., 2010). Fifteen taxa can be considered rare (\50 sites); within this group, 11 taxa are extremely rare, being in less than 15 sites (Table 1). To our knowledge, Cyzicus tetracerus (Krynicki, 1830) is the only taxon that has not been captured again since 1983 , so it was not present in further analyses. The majority of the taxa are considered as Trans-Iberian species with wide distributions across the western Palaearctic, although there are also a large number of Iberian endemic species (10 species, of which 4 can be considered as extremely rare). Spinicaudatans, notostracans and extremely rare anostracans tended to co-occur frequently with other species, with the exception of Tanymastigites lusitanica Machado & Sala, 2013, usually   sites presented co-occurrence of large branchiopod species, with decreasing site frequencies when number of co-occurring species increased (see Online Resource 4). A maximum of 6 species were found co-occurring together in 2 different sites in SW Iberian Peninsula (in Andalusia and Algarve), followed by 5 species found in 7 sites (Andalusia, Extremadura and Alentejo) and 4 species found in 34 sites (Andalusia, Extremadura, Algarve, Alentejo and Valencia). The UTM 100 km grid square approach showed that the general biodiversity patterns were geographically heterogeneous in the study area (Fig. 2).
Number of sites with large branchiopods (mean = 23.3; range 6-125; Fig. 2A) had high values in SW and NE Iberian Peninsula, similarly to the pattern found for the total species richness per UTM grid square (mean = 4.6; range 1-11; Fig. 2C). In the case of the mean species richness per site (mean = 1.4; range 1.0-2.2; Fig. 2B), high values were also more common in SW Iberian Peninsula, but there were also some UTM grid squares with high values in the eastern part of the Peninsula and in the Balearic Islands. The percentage of sites with cooccurrences within a UTM grid square (mean = 26.8;  Abellán et al., 2005) and percentage of sites where each taxon co-occurred with any other taxa. Note that Triops spp. is a heterogeneous group that includes 6 species, 4 of them endemic from the Iberian Peninsula (1 of them extremely rare) C common, EN endemic species, EX exotic species, M Mediterranean species, N northern species (present in Europe and the Iberian Peninsula, but not in north Africa), R rare, RR extremely rare, S southern species (present in north Africa and the Iberian Peninsula, but not in Europe north of the Pyrenees), T Trans-Iberian species (present in Europe, the Iberian Peninsula and north Africa) range 0.0-67.4; Fig. 2D) was especially high in SW, whereas lower values were observed in the NE and some areas in the SE Iberian Peninsula. The probabilistic co-occurrences analysis showed some species pairs with significant positive associations (e.g. Branchipus cortesi Alonso & Jaume, 1991 and Tanymastix stagnalis (L., 1758); Cyzicus grubei (Simon, 1886) and Triops spp.; Fig. 3), whereas other pairs had significant negative co-occurrences (e.g. Branchipus schaefferi Fischer, 1834 and Chirocephalus diaphanus Prévost, 1803; T. stagnalis and B. schaefferi; C. diaphanus and T. stagnalis; Fig. 3). The partitioning around medoids (PAM) analysis obtained similar results, clustering together the species that showed positive co-occurrences (Fig. 4). However, the global average silhouette width (ASW; a measure of how well the species are clustered) was quite low, indicating that in general the affinity among the species in each cluster was not very strong. PAM analysis classified the species in four significant groups: group 1 (including Artemia salina (L., 1758) and the diploid parthenogenetic Artemia), group 3 (including B. cortesi, T. stagnalis and Maghrebestheria maroccana Thiéry, 1988), group 4 (including Streptocephalus torvicornis (Waga, 1842), Triops spp., C. grubei, C. diaphanus and B. schaefferi) and group 6 [including the halophilous Branchinectella media (Schmankewitsch, 1873) and Phallocryptus spinosus (Milne-Edwards, 1840)]. Group 1 had the highest ASW, whereas group 4 had the lowest ASW, pointing that the species in this cluster were not tightly grouped.
Overall, hierarchical variation partitioning showed significant effects of landscape metrics on large branchiopods at several spatial scales. In general, the best fit was observed at the largest spatial scale (10 km). In the case of species richness, although the proportion of the variation explained was low, the presence of non-irrigated agriculture (nIA) revealed positive relationships at lower spatial scales (100 m and 1 km), whereas forests and semi-natural vegetation (FS) showed significant positive relationship at Fig. 3 Pairwise combinations of large branchiopod species in the Iberian Peninsula and Balearic Islands determined by probabilistic cooccurrence analysis. Circles indicate significant cooccurrences, and the size of the circles represents the effect sizes of the association. Only significant co-occurrences are shown. White circles in white background positive cooccurrences; black circles in grey background negative co-occurrences. Abbreviations of species are the same as in Table 1 the lowest spatial scale (100 m) but negative relationships at larger scales (Fig. 5). At the largest spatial scale, presence of wetlands and landscape diversity (Ldiv) showed a positive relationship, and population density showed a negative relationship.
In the case of group 1 of PAM analysis, the proportion of the variation explained was relatively high, especially for the diploid parthenogenetic Artemia. Wetlands (W) and population density (DP) had significant positive relationships on the distribution of both taxa at different scales, whereas the rest of metrics had low relation (Fig. 5). For group 3 of PAM analysis, composed by an assemblage characteristic of the SW Iberian Peninsula, protected areas (N2000) and Ldiv (especially at large spatial scales) showed significant positive relationships on the distribution of B. cortesi and T. stagnalis. FS and nIA showed different effects depending of the spatial scale analysed, being negative at larger scales and positive at smaller ones (Fig. 5). In contrast, the group 4 of PAM analysis (Fig. 5) was characterized by significant positive relationships of nIA, especially at small spatial scales, although the proportion of the variation explained was relatively low. Also, the distribution of some species was defined by being outside of protected areas (C. diaphanus, S. torvicornis and Triops spp.). Interestingly, B. schaefferi was positively related to population density (DP), artificial surfaces (AS), irrigated agriculture (IA), FS and nIA, whereas Chirocephalus diaphanus was negatively related to DP, IA and Ldiv. The group 6 of PAM analysis was composed of Branchinectella media and Phallocryptus spinosus, two characteristic species of saline endorrheic wetlands, and their distribution were explained by being positively related to N2000, Ldiv Average silhouette width of the overall analysis is indicated at the bottom of the plot. j stands for the number of the cluster; n j indicates the number of species within each cluster j; ave ieCj s i corresponds to the average silhouette width for each cluster. Abbreviations of species are the same as in Table 1 c Fig. 5 Plots showing species richness (S) and species groups of PAM analysis (groups 1, 3, 4, 6 and the rest of species not clustered in any group of PAM analysis) in the three spatial scales (100 m, 1 km and 10 km). The first column shows the fit of the binomial models calculated as the Pearson's correlation between observed and expected values (white bars), whereas the rest of columns show the independent contribution of the landscape metrics obtained from the hierarchical partitioning (grey bars indicate positive effects while black bars indicate negative effects). Only significant relationships detected after the randomization process are shown (see ''Materials and methods'' section for more details). Abbreviations of species are the same as in Table 1, and those of landscape metrics appear in Materials and methods (mainly at small spatial scales), nIA and W (both at large spatial scales), while FS had negative significant relationship on them (Fig. 5). The rest of the species formed a heterogeneous group of extremely rare species that did not cluster significantly in PAM analysis. Therefore, there was not a common pattern in the landscape metrics affecting the species distributions, although DP had negative significant relationships on the distribution of all species except Lepidurus apus (L., 1758), whereas nIA had positive significant relationships at larger scales with the exception of the tetraploid parthenogenetic Artemia (Fig. 5).
Wetland and non-irrigated agriculture land covers were the only factors that were related to the biodiversity metrics at the UTM 100 km grid square approach. Significantly higher values of mean species richness per site at a spatial scale of 1 km radius (Fig. 6A) were found in UTM grid squares with high values of wetland land cover (more than 6.9 ha), and also in non-irrigated agriculture land cover (more than 177.4 ha). Similar results were obtained for total Fig. 6 Significant results of the conditional inference trees between biodiversity and landscape metrics for each UTM 100 km grid square at each spatial scale. A Mean species richness per site at a spatial scale of 1 km radius. B Total species richness per UTM grid square at a spatial scale of 1 km radius.
C Total species richness per UTM grid square at a spatial scale of 10 km radius. D Percentage of sites with co-occurrences at a spatial scale of 1 km radius. Abbreviations of landscape metrics appear in ''Materials and methods'' section species richness per UTM grid square at a spatial scale of 10 km radius (Fig. 6C), and for the percentage of sites with co-occurrences at a spatial scale of 1 km radius (Fig. 6D). In the case of total species richness per UTM grid square at a spatial scale of 1 km radius (Fig. 6B), significantly higher values were found only in UTM grid squares with high values of wetland land cover (more than 2.4 ha).

Discussion
This work represents an exhaustive update of the distribution and the biodiversity patterns of the large branchiopod fauna of the Iberian Peninsula and Balearic Islands, similarly to other studies in neighbouring countries (e.g. Defaye et al., 1998;Mura, 1999;Van den Broeck et al., 2015a;Marrone et al., 2016). Since the major works by Alonso ( , 1998 Korn et al., 2010;, 2 more are currently being described (Linderiella sp. and Tanymastix sp.), 1 is new for the Iberian Peninsula (Triops simplex Ghigi, 1921; Korn et al., 2010) and an exotic species was recorded all along the Peninsula (Artemia franciscana; Amat et al., 2005). To date, the number of species found in the study area, 28 species, is similar to neighbouring countries (e.g. 22 species in Morocco (Thiéry, 1986;Amat et al., 2005;Boix et al., 2016;Korn & Hundsdoerfer, 2016), and 25 in Italy (Cottarelli & Mura, 1983;Mura, 1999;Marrone & Mura, 2006;Mura et al., 2006;Scanabissi et al., 2006;Alfonso, 2017)) and slightly higher than others [e.g. 19 species both in Tunisia (Ben Naceur et al., 2013;Marrone et al., 2016) and Algeria (Samraoui & Dumont, 2002;Samraoui et al., 2006;Ghomari et al., 2011), and 18 species in France (Nourisson & Thiéry, 1988;Defaye et al., 1998;Amat et al., 2005;Rabet et al., 2005)]. The taxa with high occurrences were species common in the western Palaearctic (C. diaphanus, B. schaefferi, T, stagnalis, S. torvicornis), together with two Iberian endemics (B. cortesi, C. grubei) widely distributed mainly in the western Peninsula. Also, the genus Triops had high occurrences, but paradoxically it is one of the taxa that gathers more uncertainties about its distribution due to recent taxonomical modifications (until 2006 it was considered as only one species, T. cancriformis (Bosc, 1801); Korn et al., 2006Korn et al., , 2010, that difficult the allocation of all previous bibliographical records to the 6 species present in the Iberian Peninsula. A detailed study of the identities, distributions and conservation status of the genus Triops in the Iberian Peninsula (especially in the northern and eastern regions) and Balearic Islands is therefore strongly recommended. For the extremely rare species, monitoring and conservation programs should be mandatory, due to their low occurrences, in some cases with the major part of their localities outside protected areas. A paradigmatic case is the anostracan Linderiella baetica , of which the only known remaining population is critically endangered by urbanization (García de Lomas et al., 2016), and that urgent measures are needed to avoid extinction. However, surveillance and monitoring programs should be also regarded for not so rare species, because their populations could be declining rapidly, even at supranational scales. This is the case of the autochthonous species of genus Artemia, which are also of great concern. The abandonment of the salterns, together with the expansion of the invasive A. franciscana, are leading to the decline and loss of populations of native species through competition and displacement , and there is an urgent need to update the distributional data and to assess the conservation status of the Mediterranean brine shrimps at national level or higher in order to evaluate the extent of damage to their populations (Mura et al., 2006).
Four large branchiopod associations were detected which corresponded in part to those already described by Alonso (1998). Two associations were characteristic of saline habitats: group 1 is a characteristic association with Artemia present mainly in salterns, whereas group 6 contained the halophile species Phallocryptus spinosus and Branchinectella media, which is also present in other areas with steppic habitats (e.g. Samraoui et al., 2006;Marrone et al., 2016). The other two associations were characteristic of freshwater temporary waterbodies. Group 3 is a characteristic association from the SW Iberian Peninsula with B. cortesi, T. stagnalis and sometimes M. maroccana, partly described in Cancela da Fonseca et al. (2008). This association is partly similar to the one observed in Morocco, where M. maroccana usually also appears with Tanymastix affinis Daday, 1910 and sometimes with B. schaefferi (Thiéry, 1991;Van den Broeck et al., 2015a). Group 4 corresponds to an association that coincides partly with the one described by Alonso (1998), but probably corresponds to two different groups due to the fact that several species of genus Triops are included in this group. In this sense, an association with B. schaefferi, Triops spp. and sometimes S. torvicornis, is present in the arid climates of the north and the east of the Peninsula, and in the Balearic Islands, whereas an association with C. diaphanus, C. grubei and Triops spp. is characteristic of the south of the Iberian Peninsula.
The necessity of increasing our knowledge on large branchiopod ecology should not be restricted only to species or population level, but they should also encompass higher levels of organization. The functioning of biodiversity patterns at different spatial scales must be taken into account for improving conservation actions, because disturbances are scale dependent (Caro, 2010). In general, our study detected that land use had an effect on the species richness and distribution of large branchiopods. However, not all the species were affected by the same landscape metrics. These results are in concordance with those observed by , with T. cancriformis being affected by irrigated agriculture at large spatial scale, whereas Branchinecta orientalis was not. In contrast, other studies found no effects of land use on large branchiopod communities (e.g. Horváth et al., 2013;Mabidi et al., 2016). We found that species clustered together were similarly affected by the same landscape metrics, but that different clusters were affected by different landscape metrics, which could partially explain those contrasting results. Agricultural (irrigated and non-irrigated) and wetland land covers were the metrics that affected to more species (positive or negative, and with different magnitudes). It is known that agricultural land use around temporary ponds can have important effects on their biota (Rhazi et al., 2001;O'Neill et al., 2016), although traditional agricultural practices are a key factor for temporary pond conservation (Bagella et al., 2016). This is in concordance with the positive relationships that showed the non-irrigated agriculture on the majority of large branchiopod species. Similarly, it is known that the extent of wet areas or the pond density also influences invertebrate and plant assemblages (e.g. Boix et al., 2008;Bagella et al., 2010). It is interesting to note that the majority of the relationships had a better fit at larger spatial scales, in accordance with the results found by , suggesting that management at larger scales have also to be taken into account for the conservation of large branchiopods.
Moreover, biodiversity metrics were influenced by the presence of wetlands and/or land use linked to nonirrigated agricultural practices, especially at larger scales. Our results showed that wetlands, floodplains and surrounding areas in the Iberian Peninsula (e.g. Guadiana, Doñana, Cádiz and Empordà wetlands) presented high values of large branchiopod species richness at different spatial scales, probably due to the heterogeneity of the water body types that could host different communities. Moreover, these areas were also important for the presence of rare and/or Iberian endemic species. The importance of this kind of habitats for the biodiversity of large branchiopods has also been observed in other regions (e.g. Eder et al., 1997;Timms & Sanders, 2002;Waterkeyn et al., 2009;Nhiwatiwa et al., 2014;Stoch et al., 2016). However, land-use analyses also revealed that temporary water bodies in non-irrigated agricultural landscapes not linked directly to large wetland areas also hosted a large amount of biodiversity. Traditional, non-intensive agricultural practices are compatible with the existence of highly biodiverse temporary ponds (e.g. Grillas et al., 2004;Robson & Clay, 2005), even in areas where temporary ponds are human-made and the presence of natural temporary water bodies is very rare (e.g. . The large branchiopod fauna associated to these habitats is usually very rich, and they host multiple species (e.g. Boven et al., 2008;Cancela da Fonseca et al., 2008;Marrone et al., 2016;Alfonso, 2017). However, these habitats are rapidly declining due to intensification of agricultural practices, changes in land use and population growth (e.g. Underwood et al., 2009;Rhazi et al., 2012), affecting the biotic populations linked to these habitats (e.g. Euliss Jr. & Mushet, 1999;Beja & Alcazar, 2003;Van den Broeck et al., 2015a).
Since the extensive surveys performed by Alonso ( , 1998 in the Iberian Peninsula and Balearic Islands during the last quarter of the twentieth century, no more studies exist updating the information for the whole study area. In this sense, the present study represents a new step on the knowledge of large branchiopod distribution in this area, not only because it gathers all available information on the distribution of this taxonomic group (historical records and new observations), but also because it analyses the effect of several landscape characteristics on its distribution, detailing the spatial scale at which its influence is noticeable for large branchiopods presence. All this information would be highly valuable for land stakeholders under future conservation scenarios.

Online Resource 2.
Distribution maps of all records (including those records before 1995) of large branchiopods, and for each large branchiopod taxon present in Iberian Peninsula and Balearic Islands. p.