Seed plant features, distribution patterns, diversity hotspots, and conservation gaps in Xinjiang, China

The flora in Xinjiang is unique. Decisions about biodiversity conservation and management based on seed plant diversity hotspots and conservation gaps in Xinjiang are essential to maintain this unique flora. Based on a species distribution dataset of seed plants, we measured seed plant diversity using species richness and phylogenetic diversity indices. Five percent of Xinjiang’s total land area with the highest biodiversity was used to identify hotspots for each index. In total, eight hotspots were identified. Most hotspots were located in mountainous areas, mainly in the Tianshan Mountains and Altai Mountains. Furthermore, we detected conservation gaps for Xinjiang’s seed flora hotspots by overlaying nature reserve maps on to maps of identified hotspots and we designated priority conservation gaps for hotspots by overlaying global biodiversity hotspot maps on to hotspot conservation gaps maps. Most of Xinjiang’s seed plant hotspots are poorly protected; only 10.45% of these hotspots were covered by nature reserves. We suggest that it is essential to promote network function of nature reserves within these hotspots in Xinjiang to conserve this unique flora. Nature Conservation 27: 1–15 (2018) doi: 10.3897/natureconservation.27.23728 http://natureconservation.pensoft.net Copyright Jihong Huang et al. This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited. RESEARCH ARTICLE Launched to accelerate biodiversity conservation A peer-reviewed open-access journal


Introduction
Biodiversity is unevenly distributed; therefore, prioritisation is essential to minimise biodiversity loss (Brooks et al. 2006). Identifying priority areas or hotspots where biodiversity is most threatened is critical for effective conservation (Myers et al. 2000;Olson and Dinerstein 1998). Biodiversity hotspots are characterised by exceptional concentrations of endemic species that experience relatively high rates of habitat loss (Myers et al. 2000;Prendergast et al. 1999). This approach to defining hotspots has been used in many studies (Mittermeier et al. 2005). To adequately define a hotspot, its boundaries must be delineated. A predefined proportion of areas with large amounts of both endemic species richness and habitat loss was previously considered a hotspot (Orme et al. 2005;Prendergast et al. 1993). Recently, the definition of the term 'hotspot' for biodiversity has been generalised. In addition to the amount of endemic species and habitat loss, other aspects, including total number of species, number of threatened species and evolutionary history, have been proposed to identify hotspots (Cadotte et al. 2010;Forest et al. 2007;Rosauer and Laffan 2008). Researchers previously focused taxonomic diversity (Brooks et al. 2006;Reid 1998). However, more recently, many researchers focused on the evolutionary processes (Cadotte et al. 2010;Huang et al. 2016c;Xu et al. 2017). Phylogenies (evolutionary trees) have become increasingly important for identifying biodiversity hotspots. Phylogenetic diversity is often favoured because it shows relationships between extant species and provides more information about evolution over a long time-scale (Omland et al. 2008).
Xinjiang Uygur Autonomous Region (Xinjiang hereafter) is a large (1.64 million km 2 ), topographically varied inland province in northwest China that extends from high mountains (8611 m above sea level) to lowland basins (156 m below sea level) and humid forest to grasslands and dry desert. Consequently, this region has an extreme arid and cold climate and exhibits striking climate, soil and vegetation patterns. Species richness patterns in Xinjiang were previously reported for woody plants (Zhang and Zhang 2014). The vegetation types in Xinjiang are mainly desert and grassland (including steppe and meadow) and forest covers relatively small areas. The constructive species of the grassland are Gramineae, Compositae and Leguminosae, whereas the constructive species of the forest include the coniferous Larix, Picea, Abies, and Pinus; broadleaved Betula and Populus; and the native fruit species of Malus, Armeniaca and Juglans, which are endemic to Ily Valley (Xinjiang Investigation Group of Chinese Academy of Sciences 1978). Although many studies have been conducted on Xinjiang's flora, data are still not sufficient to provide systematic guidance for regional biodiversity conservation.
In this paper, based on native seed plant species in Xinjiang, we identified seed plant hotspots using both species richness and phylogenetic diversity indices. Moreover, Xinjiang's seed plants conservation gaps, which are hotspot areas that are not covered by China's nature reserves, were detected by overlaying the seed plant hotspot map on a map of Xinjiang's nature reserves; and Xinjiang seed plant priority conservation gaps, which are areas within both Xinjiang seed plant hotspots and global biodiversity hotspots that are not covered by Xinjiang's nature reserves, were detected by overlaying the Xinjiang seed plant conservation gap map on a global biodiversity hotspot map. Finally, we discussed the potential causes of the patterns revealed in this study and proposed suggestions for Xinjiang flora protection.

Data sets
We compiled a species distribution database that included all seed plant species distributed in Xinjiang based on a large number of herbarium specimens (http://www. cvh.org.cn/), a list of Chinese seed plant species and distribution information at provincial level (Wu et al. 1994(Wu et al. -2012 and other literature sources that mainly included Flora Xinjiangensis (Commissione Redactorum Florae Xinjiangensis 1992-2011), Sylva Xinjiangensis (Yang 2012), Desert Plants in China (Lu et al. 2012) and Rare Endangered Endemic Higher Plants in Xinjiang of China (Yin 2006). The county was used as the basic spatial unit for species distributions and geographical coordinates were included for each specimen. To improve the quality of the distribution data at the county level, we overlaid county distributions over altitudinal range distributions for each species. We divided the entire region into 50×50 km grids and then overlaid these grids on to the distribution polygons or points. Finally, we acquired presence/absence data for each species in each grid. In our study, data on 3715 seed plant species were obtained based on a total of 145781 occurrence records.
A georeferenced or spatially referenced dataset of 33 nature reserves in Xinjiang was constructed. The main body of nature reserves was downloaded from the World Database on Protected Areas (http://www.protectedplanet.net/). The global biodiversity hotspots dataset was downloaded from Conservation International (http://www. conservation.org/How/Pages/Hotspots.aspx).

Diversity measures and phylogenetic tree construction
To detect the distribution patterns of Xinjiang seed flora, we calculated seed species richness and phylogenetic diversity of seed plants. Phylogenetic diversity is a simple measure of evolutionary diversity that is determined based on the sum of the lengths of all branches that are members of the corresponding minimum spanning path (Faith 1992). The phylogenetic diversity equation is: where C indicates the set of branches in the minimum spanning path joining the species to the root of the tree, c indicates a branch in the spanning path, and c L indicates the length of branch c. Phylogenetic diversities were calculated by using the phylogenetic diversity algorithm in PHYLOCOM (Webb et al. 2008).
PHYLOMATIC (Webb and Donoghue 2005) was used to construct a phylogenetic supertree for Xinjiang's seed plants. The backbone of the supertree in PHYLO-MATIC was based on the latest plant phylogeny group classification. The BLADJ algorithm (Webb et al. 2008) was used to adjust the branch lengths of our phylogenetic tree by using known molecular and fossil dates (Wikstrom et al. 2001).

Hotspot identification and ranking and conservation gap analysis
Both diversity indices were calculated using two packages, Ape and Picante, in R 3.3.2 (R Core Team 2016) and the resulting maps were generated with ARCGIS 9.0 (ESRI, Redlands, CA, USA). The relationship between the patterns of diversity characterised with the different indices was calculated using Spearman's correlation coefficient for each pair of indices. Hotspots were identified with each diversity index as indicated hereafter. We accumulated all of the cells with the highest values of the two diversity indices until the total area reached 5% of Xinjiang's total land area. It was well known that complementarity (Vane-Wright et al. 1991) plays an important role in biodiversity hotspot identification (Zhang and Ma 2008). Here, we did not use complementarity (Vane-Wright et al. 1991) to identify hotspots, because our method ensures that all local areas (i.e. cells) within the hotspots have more species than any other areas outside the hotspots. Additionally, our grid cell size was 50km × 50km, which is large enough for local authorities to manage.
We then identified 'conservation gaps ' (Huang et al. 2016a) and 'priority conservation gaps ' (Huang et al. 2016a) for Xinjiang seed plants. We synthesised species richness, hotspot size, the proportion of protected areas that overlap with hotspots and global role in biodiversity conservation to rank all hotspots identified in our study for conservation. Species richness was calculated by species number in each hotspot; hotspot size was determined by the area of each hotspot; the proportion of protected areas that overlap with hotspots was calculated based on the each hotspot divided by the area covered by nature reverses; and global role in biodiversity conservation was calculated based on the area of each hotspot divided by the area covered by global biodiversity hotspots. These four indices were ranked and divided by total number of hotspots then summed to produce a total value. All hotspots were ranked according to their total value.

Xinjiang seed plant features
Based on our data, Xinjiang has 3715 seed plant species (including 66 subspecies and 224 varieties), which account for 12.71% of all seed plant species in China. These seed plants belong to 767 genera and 113 families. Of these species, gymnosperms, dicots and monocots accounted for 1.10%, 80.57% and 18.33%, respectively, of all Xinjiang seed plant species (Table 1).

Xinjiang seed plants distribution and similarity at the regional level
Species richness at the regional scale generally increased with decreasing latitude in Xinjiang. In the 10 geographical regions (Fig. 1), region 9 (South-Central Tarim Basin) had Table 1. Numbers of seed plant orders, families, genera and species in each geographical region of Xinjiang. Regional division of Xinjiang is shown in Fig. 1  the lowest species richness (only 264 species) and region 2 (South of Junggar Basin) had the highest species richness (2445 species) (Table 1). In terms of regional occurrence, only 61 species (1.64% of total species) were found in each of the 10 geographical regions. By comparison, 1046 species (28.15% of total species) occurred in only one of the 10 geographical regions. The number of species that were present in a given number of regions decreased as the number of regions increased, which demonstrates that seed plants species with narrow spatial ranges were more common than those with wide or transcontinental spatial regions (Fig. 2). With the hyperbola model (G=61.07+1060.20/RT), where G was the number of species and RT was the total number of regions ( Fig. 2), the spatial range of the genetic distributions of seed plants explained 89% of the variation in the number of species with varied distribution ranges. The similarity indices for pairwise comparisons of regional species composition ranged from 0.07 to 0.58 (Table 2). In general, region 9 was the least similar to all other regions. This trend was, in part, due to the low number of species in this region. The lowest similarity indices were found between the regions 9 and 1 and 9 and 3, whereas the highest was found between the regions 1 and 2 (Table 2).

Diversity hotspots of Xinjiang seed plants
Distributions of Xinjiang seed plant species richness and phylogenetic diversity were uneven and concentrated in the northern regions of Xinjiang (i.e. north of the Tianshan Figu 2. Species richness distribution for 3715 species in Xinjiang across the 10 spatial range classes. Spatial range was measured based on the number of geographical regions, without specific reference to particular regions. The line represents the least-squares regression using a hyperbola model with 61.07 and 1060.20 as the a and b parameters, respectively (R 2 =0.89, P<0.001). Table 2. Jaccard similarity indices for pairwise comparisons of geographical regions using all Xinjiang seed plant species. Regional division of Xinjiang is shown in Fig. 1 Mountains; Fig. 3a-c). Both species richness and phylogenetic diversity indices calculated for each spatial unit across Xinjiang were highly significantly correlated with each other (R=0.99, P<0.0001). The diversity hotspots identified based on Xinjiang seed plant species richness were distributed in seven distinct centres: H1-H7 (Table 3, Fig.  3d). The diversity hotspots identified based on phylogenetic diversity included eight centres, which included all of the above seven hotspots and H8 (Table 3, Fig. 3e). All hotspots were located north of the Tianshan Mountains (Fig. 3f ).  Figure 3. Geographic distribution of Xinjiang seed plant diversity patterns and hotspots and major mountains ranges in northwest China. a Geographic distribution of species richness b geographic distribution of phylogenetic diversity c major mountains ranges d diversity hotspots identified by species richness e diversity hotspots identified by phylogenetic diversity and f diversity hotspots identified by both species richness and phylogenetic diversity. Albers projection. Hotspot centre codes are consistent with those listed in Table 3.

Priority conservation gap areas for Xinjiang seed plant hotspots
The spatial distributions of 33 nature reserves and eight hotspots of Xinjiang seed plant species significantly differed (Fig. 4). These reserves and hotspots overlapped across 11238.63 km 2 , which represented 4.18% of the total area of Xinjiang's nature reserves; 89.55% of the hotspot areas were not covered by nature reserves. Therefore, most of the seed plant hotspots are not well protected and thus represent conservation gaps (Fig. 4).  Table 3.
Four of the eight hotspots overlapped with global biodiversity hotspots (H2, H3, H5 and H7; 60000.00 km 2 ); these overlapping areas accounted for 55.81% of the total area of all hotspots and contained 59.11% (2192 species) of all Xinjiang seed plant species. The remaining four hotspots harboured 1836 seed plant species, which accounted for 49.42% of all Xinjiang seed plant species. A priority conservation area that could alleviate Xinjiang seed plant conservation gaps is centre H5 (Fig. 4). Xinjiang's nature reserves, seed plant hotspots and global biodiversity hotspots overlap across 3319.79 km 2 and only account for 1.23% of the total area of Xinjiang's nature reserves.

Xinjiang flora uniqueness and complexity
Xinjiang is located in the contact area of several natural geographical units, including the Altai Mountains, Tianshan Mountains, Pamirs, Kunlun Mountains, Altun Mountains and Northern Tibetan Plateau. Moreover, phytogeographically, Xinjiang also has intersections of Eurasian Forest, Eurasian steppe, Central Asian desert and Himalayan Chinese plant sub-regions. These characteristics have led to floral uniqueness and vegetation complexity in Xinjiang (Xinjiang Investigation Group of Chinese Academy of Sciences 1978). It was reported that there are 4081 species of vascular plants in Xinjiang (including subspecies), which account for 13.4% of the total number of plants in China (Yin 2006). Although, the number of species in Xinjiang is relatively low, the proportion of species that are locally endemic to Xinjiang is high (Huang et al. 2016b;Huang et al. 2014).

Xinjiang seed plant hotspots
Species richness is higher in Xinjiang mountains than basins (Li et al. 2013). The eight Xinjiang seed plant hotspots identified in this study cover key mountain ranges in Xinjiang, including both the Tianshan Mountains and Altai Mountains. Wang et al. (1993) proposed that these two terrestrial areas are key for conservation of China's biodiversity and these areas represent two of the 32 terrestrial priority areas for biodiversity conservation in China developed by the Ministry of Environmental Protection (Editorial Committee of China National Biodiversity Conservation Strategy and Action Plan 2011) and one of 19 centres of Chinese endemic seed plant hotspots (Huang et al. 2016a). Thus, Xinjiang seed plant hotspots encompass key areas for biodiversity conservation in China. This finding indicates that Xinjiang's seed plants represent a valuable group for identifying hotspots or priority areas that can be used to aid Xinjiang biodiversity conservation.

Hotspots in major mountain ranges of Xinjiang
Seed plant species are unevenly distributed across Xinjiang; all their hotspots are located in mountain areas mainly within the Tianshan Mountains, Altai Mountains and western mountainous areas in Junggar. These areas approximately correspond to the montane forest and cold temperate coniferous forest in temperate steppe and temperate desert regions (Chen 2014). Mountainous areas were potentially important for shaping the distribution of Xinjiang seed plant species; the high level of heterogeneity in the physiognomy could have offered a diversity of environments and thus promoted formation of the substantial species diversity (Qian and Ricklefs 1999) and caused aggregation of seed plants in the mountain regions of Xinjiang. The Tianshan Mountains and Altai Mountains are two natural physical barriers that alter atmospheric movement and consequently cause changes in precipitation and temperature. These conditions maintain a greater variety of habitats and thus probably enhance speciation, differentiation and preservation of the species that live in these areas, because species diversity increases with habitat variation (Latham and Ricklefs 1993;Rosenzweig 1995). These natural physical barriers also restrict north-south and east-west plant migration and interchange by forming a geographically isolated area (Wu and Wang 1983), which thereby facilitates the speciation and differentiation of endemics. Alternatively, other factors may play prominent but not decisive roles in some local regions, such as transitional areas (Axelrod et al. 1996;Wu 1980;Wu and Wang 1983) and islands isolated from the continents (Lomolino et al. 2006). Consequently, the abundance of seed plant species in mountainous forests results from a combination of effects and is mainly attributable to the heterogeneous terrain, geological history and contemporary environments in this region.

Suggestions for Xinjiang flora protection
Xinjiang is the largest province in China and accounts for one-sixth of China's total areas (Xinjiang Investigation Group of Chinese Academy of Sciences 1978). Xinjiang has extensively varied geographical and topographical features and a complex and ancient geographical history. Although plants are not abundant (Yin 2006), Xinjiang has unique flora and complex vegetation. Therefore, it is necessary and urgent to protect plants in Xinjiang. Based on the results of this study, some suggestions for conservation of the Xinjiang flora are proposed.
Firstly, we need to create a unified conservation plan for the Xinjiang flora. Conservation biologists have proposed many conservation suggestions and schemes across China (Li et al. 1993;Liu et al. 2003;Sang et al. 2011;Wu et al. 2011). However, these suggestions have unfortunately not been perfected and implemented at the provincial level, mainly because solid data support for plant diversity conservation in Xinjiang is lacking. Although research on species diversity in Xinjiang is very limited, researchers have conducted a large number of taxonomic investigations on plants in Xinjiang (CRFX 1992. Secondly, we need to establish a provincal georeferenced (spatially referenced) database for biodiversity that includes plant, animal and nature reserve distributions. In particular, efficient conservation of Xinjiang's plant resources requires accurate data on the current distribution and threat status of plant species in the country (Sang et al. 2011). A large amount of plant diversity information has been digitised in China, including through the National Specimen Information Infrastructure (http://www.nsii. org.cn/), Chinese Virtual Herbarium (http://www.cvh.org.cn/), Chinese Field Herbarium (http://www.cfh.ac.cn/), the Flora of China project (http://hua.huh.harvard.edu/ china), Scientific Database of China Plant Species (http://db.kib.ac.cn/eflora/Default. aspx). Some datasets have also been published with Chinese administrative counties as the spatial units, such as in studies on wood plants in China (Fang et al. 2009) and Chinese endemic seed plants (Huang et al. 2014). Additionally, by the end of 2014, China established 2729 nature reserves (Ministry of Environmental Protection the People's Republic of China 2015). Although we produced a spatially referenced dataset of 2139 nature reserves (the World Database on Protected Areas: http://www.protectedplanet. net/), the boundaries of many nature reserves were not accurate, because they did not have clear boundaries. Moreover, although a vegetation map of China was published in 2008 (Editorial Committee of Vegetation Map China 2008), many vegetation field survey projects were completed 30-50 years ago. All of the described studies provide an important data for research on biodiversity conservation across China. However, the data provided by those projects could not meet the precision requirements for biodiversity conservation research at the provincial level; therefore, it is necessary to develop additional solid sources of information for Xinjiang plants.

Conclusions
Different diversity indices revealed similar distribution patterns of plant diversity in Xinjiang Province. A total of eight seed plant hotspots in Xinjiang were identified by species richness and phylogenetic diversity. All of these hotspots are located in the mountainous areas within the Tianshan Mountains and Altai Mountains. Most centres of Xinjiang seed plant hotspots are not covered by the current nature reserve system. How large spatially-explicit optimal reserve design models can we solve now? An exploration of current models' computational efficiency

Introduction
Anthropogenic activities have caused tremendous impacts on ecosystems all over the world. These impacts include but are not limited to environmental pollution, habitat loss and fragmentation, invasion of exotic species and climate change (e.g. Aplet and McKinley 2017). These impacts led to an increased rate of extinction of species and impaired the services that ecosystems have been providing to humans. Establishing nature conservation reserves have been adopted across the globe as a direct way to restore ecosystems and protect species. However, resources devoted to this purpose such as land and budget have always been scarce and they often face conflicting demands from other sectors of society such as housing and industry. With the aim of designing ecologically effective and economically efficient nature reserves comes the science of reserve design (Kingsland 2002).
The reserve design problem is defined as selecting a subset from a larger set of candidate sites to assemble nature conservation reserves to provide adequate protection to targeted species or habitats. This problem has been studied using various approaches, including site scoring (e.g. Pressey and Nicholls 1989), gap analysis (e.g. Scott et al. 1993), heuristics (e.g. Pressey et al. 1997) and formal optimisation (e.g. Possingham et al. 2000, Williams et al. 2004). Gap analysis and heuristics have been widely used both in academic research and conservation practices, primarily because they are intuitive, easy to use and computationally convenient. Some popular nature reserve design software such as Marxan and C-plan use heuristics as their major solution algorithm (Ball et al. 2009, Sarkar 2012. However, it has long been demonstrated that these techniques find optimal solutions only occasionally and often generate solutions which significantly deviate from the optimal ones (Önal 2003, Williams et al. 2005). This means that they may not be allocating scarce conservation resources in the best possible manner. The problem can also be formulated as an optimisation model, specifically as a linear mixed integer programme (MIP), to determine the best selection of sites while meeting the conservation goals. The resulting model can be solved using MIP solvers such as CPLEX and GUROBI integrated in off-theshelf mathematical modelling software such as GAMS and AMPL (AMPL 2016, GAMS 2016. Finding the best possible solution helps guarantee efficient allocation of conservation resources. Due to this advantage, optimal reserve design models have received substantial attention since the 1990s. The early optimisation formulations of the reserve design problem adapted the set covering and maximal covering frameworks well known in the operations research literature (Cocks and Baird 1989, Underhill 1994, Camm et al. 1996, Church et al. 1996, Ando et al. 1998). These initial formulations neglected spatial attributes of the reserve, usually ending up with a solution where selected sites are highly scattered over the landscape and far apart from each other. A reserve with such a spatial configuration is not very meaningful in practice because both species survival and reserve management benefit from a spatially coherent reserve. The second generation optimum reserve design models incorporated spatial attributes to improve the coherence and functionality of the reserves (Williams et al. 2005). Amongst the spatial attributes discussed in literature, the most important ones are spatial contiguity (or connectivity) and compactness of the selected sites. Spatial contiguity means that any two sites in the reserve can be linked by a path of selected sites. Compactness is a more difficult concept to define, but, in general, it means that the selected sites must be close to each other (or grouped together) to the greatest possible extent. Some studies used the total distance between pairs of selected sites and the boundary length of the reserve as proxy measures for compactness. These two spatial attributes have been successfully modelled using graph theory, network flow theory and mathematical techniques (see Williams et al. 2005, Moilanen et al. 2009, Billionnet 2013 for reviews of those studies). Incorporating these spatial attributes into the optimisation frameworks requires more complex formulations, however, including a large number of additional variables and constraints. Thus, the reserve selection problem gets computationally more difficult when spatial configuration is of concern. For instance, a set covering formulation involving thousands of candidate sites could be solved in seconds (Rodrigues and Gaston 2002), but a compactness model involving about 700 sites could not be solved to optimality within 10 hours (Fischer and Church 2003).
Several modelling approaches have been introduced in the reserve design literature to incorporate spatial aspects in site selection. Due to their differences in structure and size, those models differ in their computational performance. Moreover, individual models were run on different platforms with different CPU speeds and solved using different MIP solvers and, therefore, it is not clear how those models compare to each other in terms of their computational efficiency. A fair and accurate comparison of the convenience of these models would be important for on-theground conservation practitioners. This is the main motivation for this article. For this purpose, we selected eight optimisation models incorporating (or could incorporate easily by proper modification) contiguity and compactness and conducted computational experiments with these models by running them on the same computer and using the same data sets that were generated randomly. We aimed to answer two questions: i) how do these models perform in solving contiguity and compactness problems with different sizes and ii) what are the likely factors affecting their computational efficiencies?
We carried out the tests on a Lenovo desktop computer with Intel Pentium (R) CPU of 2.8 gigahertz and RAM of 8.0 gigabyte and the operation system was Windows 64-bit. As the MIP solver, GUROBI 5.0 integrated in GAMS 23.9 was used. Table 1 summarises the eight models we tested. The 'Algebraic formula' column of the table presents the basic algebraic formula used to ensure spatial contiguity.

Contiguity models
Some of the models were not originally proposed for reserve design purposes or did not consider the cost of site selection. We made necessary modifications to those models so that they are applicable to the reserve design problem we address here. The objective of all the models is stated as the minimisation of the total cost of site selection, namely: where U i is a binary variable which equals 1 if site i is selected and 0 otherwise and c i is the cost of selecting site i.
Those models, not originally proposed for reserve design problems, are extended by adding the species coverage constraint described below: where e i is the population size of target species in site i; and p is the targeted (minimum viable) population of target species. This constraint is supposed to ensure a longterm survival chance for target species.
The Miller et al. (1960) model is modified as below to control the formation of arcs between selected sites.
where X ij is a binary variable which equals 1 if an arc directs from i to neighbour j and 0 otherwise and N i is the set of i's neighbour sites. For details, see Williams (2002).
We first generated ten hypothetical data sets. Each data set considers a 10*10 grid partition including one hundred sites. The distribution of species over these sites was generated by using a uniform distribution function with the range of [0, 5], which means that the population size (e i ) in each site takes an integer value between 0 and 5. We generated the selection cost (c i ) similarly using the range of real numbers [1, 10]. We assigned a different value to the target population size (p) in each data set, where p is calculated by p=α*n where n=100 is the number of candidate sites and α is a popula- Table 1. Mechanisms used to ensure contiguity and algebraic formulae of the eight models used in the computational efficiency tests.

Model number
Mechanism to ensure contiguity Algebraic formula References 1 An arc can direct from node i to neighbour node j only if the value assigned to i is less than the value assigned to j, so an arc directing away from a node will not return to it and therefore loops are prevented.
(1)-(6); Two trees are formed, one in primal graph, the other in dual graph, the interwoven structure of these two trees prevents the formation of loops.
(1)-(6); (2)-(4) in Williams (2002). Williams (2002) 3 One selected site serves as a sink, every other selected site provides one unit of supply to the flow and all flows finally reach the sink.
(1) and (2) The tail length of a selected site equals the number of arcs linking to it and that length is strictly larger than that of the site's preceding sites.
(1)-(2); (3)-(7) in Önal and Briers (2006). Önal and Briers (2006) 5 Solve a relaxed problem first with no contiguity constraints. If a cycle occurs, apply an inequality (a cut) to that cycle and solve again until there is no cycle in the solution.
Önal and Wang (2008) 6 An amount of flow is injected into the partition from an outside node, each selected site consumes one unit of flow, residual flow is absorbed by a predefined variable.
(1) and (2); (5)-(8) in Conrad et al. (2012). Conrad et al. (2012) 7 Capital flows from the supply node to trans-shipment nodes and for each selected site, capital can come into it from only one of its adjacent sites.
Jafari and Hearne (2013); Jafari et al. (2017) 8 If a site is selected then at least one of its adjacent sites, which is nearer (structurally or functionally) to the centre, must also be selected.
(2),(4), (6) and (8) Wang and Önal (2016) tion coefficient ranging from 0.25 to 2.5 with an increment of 0.25 (thus, a population between 25 and 250 individuals must be protected in the reserve and a larger α implies a larger reserve area). By designing reserves with different areas, we intend to test the impact of the number of selected sites on computational efficiency.
To eliminate the impact of input data on the test result, for each value of α we generated 50 data sets, each with a different distribution of species and selection cost across the 100 sites and used in a different run. To save time in the test and avoid tying the computer to a long computational time, in each run we imposed a one-hour maxi-mum processing time limit. The total time limit for all runs was limited to 2 hours (i.e. the programme was terminated if it took more than 2 hours to complete the 50 runs). Average CPU times in completed runs were reported.

Contiguity and compactness models
The Önal et al. (2016) model and the Williams (2002) model were modified and used in the contiguity and compactness test because these two models were found to be the most efficient models in the above contiguity test. We modified the objective function of the Önal et al. (2016) model as below so it incorporated selection cost which is an important factor in reserve design practices: where X ki is a binary variable which equals 1 if site i is selected into the reserve which is centred at site k and d ki is the distance between centre site k and selected site i. The first part of (7) is the total distance between the cluster centres and selected sites and the second part is selection cost. These two components can be given different weights to reflect their importance in site selection; for computational purposes, we used the same weights (=1.0).
In the Williams (2002) model, we made the following modifications to incorporate compactness assuming that compactness is measured by total pair-wise distance between selected sites (Önal and Briers 2002). s.t.: where W ij is a binary variable which equals 1 when sites i and j are both selected and 0 otherwise. The first part of (8) represents reserve compactness while the second part is the cost of site selection. Constraint (9) relates W ij to the site selection variables.

Contiguity: the 100-site problem
The computational results of the 100-site contiguity problem are summarised in Table 2 (for clarity reasons, we only report results with populations (p) of 50 to 250 with an increment of 50). Except Williams (2002) Table 2. Computational efficiency of alternative contiguity models in solving a 100-site problem.

Models
Population size to be covered in the reserve 50 100 150 200 250 Miller et al. (1960) 136.9 >3600.0 a >3600.0 a >3600.0 a >2400.0 b Williams (2002) 9.0 50.8 20.2 3.0 0.1 c Shirabe (2005) 64.2 >969.1 a >477.1 a 47.0 0.1 c Önal and Briers (2006) >2024.3 a >3600.0 a >3600.0 a >3600.0 a >2400.6 b Önal and Wang (2008) >505.6 a >944.4 a >824.2 a >2405.3 a >1224.9 b Conrad et al. (2012) >1854.8 a >2551.1 a >2026.5 a >123.5 a 0.2 c Jafari and Hearne (2013) >1905.5 a >3600.0 a >3600.0 a >311.8 a 0.9 c Önal et al. (2016) 12.0 5.3 3.3 4.2 6.6 c a: 2 to 50 runs were completed within the total time limit (2 hours). These values are averages of the CPU times used in completed runs. b: 3 to 11 runs were completed within the total time limit (2 hours), amongst which 1 to 5 runs were infeasible due to high coverage requirement and/or inadequate species occurrence. The reported values are averages of the CPU times used in feasible runs. The ">" indicates that at least one run could not be completed within 1 hour. c: 22 runs were infeasible. These values are averages of the other feasible 28 runs.
reported by Williams (2002) also. The Önal et al. (2016) model was computationally steady regardless of the population size specification (the run times were in the range of 3.3-12.0 seconds).

Contiguity: larger size problems
Having observed that, in the 100-site test problems, the Williams (2002) and Önal et al. (2016) models are substantially more efficient than the other models, in the remaining test runs we compared these two models against each other in solving larger problems that involved 200 to 2000 sites. We relaxed the optimality requirement by setting the relative gap (percentage difference between the incumbent solution and the best possible solution) at 1 percent level (instead of 0 percent, no gap) to avoid long processing times that might possibly be needed to solve large problems to exact optimality (this is done in GAMS by writing a simple code option optcr=0.01). Table 3 reports the results for different model sizes.
The Williams (2002) model solved the problem for almost all sizes when the population coefficient α was high (2.5 and 2.0), although the processing times increased as the number of candidate sites (n) was increased. For instance, when n=200 and α =2.5, the model completed all feasible runs (in some cases the model was infeasible due to high population requirement and/or inadequate population distribution in the candidate sites) with an average processing time of 0.2 second. For n=2000, the average processing time increased to 4.2 seconds. A similar pattern is observed when the population coefficient was decreased to α =2.0. For n=200, an average of 4.0 seconds was used, while for n= 2000, only 14 runs were completed within the 2 hour total time limit. For other population settings, the model encountered computational difficulty more often and could not find a solution after running 1 hour in at least one run (indicated in Table 3 with ">"). The computational difficulty has become more severe as the problem size increased. For the setting α =1.5 and n=200, for example, the model completed 14 runs (average processing time =552.3 seconds). For α =1.5 and n=500, the model completed only 4 runs (average processing time = 2124.9 seconds). The model completed only 2 runs when n=1000, 1500 and 2000. The worst performance was observed when α =1.0 and 0.5, where the model ran out of the 1 hour time limit almost in each run even with the smallest test problem (n=200).
The Önal et al. (2016) model displayed a relatively steady computational performance as in the 100-site problem. In general, the model used less time when a smaller number of sites was selected (Table 3, n=200). However, it encountered more difficulties than the Williams (2002) model in solving large problems even when the population settings were high. For instance, when n=500, the model completed only two runs for each population setting; and when n=1000, 1500 and 2000, it stopped without finding any solution within the 1 hour time limit. A more serious problem was encountered in some cases where the computer ran out of memory due to the large size of the model (millions of variables and equations, see Table 3).

Contiguity and compactness
In this section, we compare the computational efficiency of the Williams (2002) and Önal et al. (2016) models in solving contiguity and compactness problems. We called the model based on Williams (2002) but with (8) as the objective function and (9) as an additional constraint "W-CM" (CM for contiguity and compactness) and the model based on Önal et al. (2016) but with (7) as the objective function "Ö-CM". The number of sites in the test runs was set as n=100, 200 and 500 and we defined the population coefficient α as in the contiguity tests. When n=500, both models encountered severe computational difficulty for most population coefficient settings. Table 4 summarises the results. a: Population coefficient. The size of population to be covered in the reserve equals to α times the number of candidate sites (n). b: The first number is the number of variables, the second number is the number of equations in the model. c: 17 to 29 runs are infeasible due to high coverage requirement and/or inadequate species occurrence.
Values are average of times used in feasible runs. d: 2 to 19 runs were completed in 2 hours. The ">" indicates at least one run was out of time limit (1 hour). e: 2 or more runs were completed, at least one of which terminated with no solution after running 1 hour or with integer solution but gap not available. f: Only 1 run was completed which took more than 7 hours and found integer solution with gap over 80 percent or no solution returned. g: Out of memory after running about 2000 seconds or integer solution found after running nearly 3 hours but with gap greater than 70 percent.
The W-CM model encountered difficulty for some population settings ( α ) when solving the 100-site problem. For α =2.5, the solution time was only 1.4 second, but for the settings α =2.0 and 1.5, it jumped to 178.2 and 1692.3 seconds. For α =1.0, only 3 runs were completed within the 2 hour total time limit. However, for α =0.5, the solution time was decreased dramatically to 90.6 seconds, showing the fast-slowfast pattern observed when solving the contiguity problems. For n=200 and 500, the model completed only 2 runs within the 2 hour total time limit for most population settings. In general, larger models require more solution times. However, raising the population requirement level improved the computational efficiency. For example, for α =2.5, the average solution times were only 4.8 and 78.3 seconds for n=200 and n=500, respectively.
The Ö-CM model used only about 2 seconds in the 100-site problem for all population coefficient settings. When n=200, the model used about 20 minutes or less for all population settings, again better than the W-CM model. When n was increased to 500, this model encountered difficulty too, where only 2 runs could be completed within the 2 hour total time limit.
It is interesting that the Ö-CM model is computationally more efficient than the W-CM model although it has a larger size both in terms of the number of variables and  Table 3 for the explanations of α , model size and symbol ">". b: Some runs were infeasible due to high population requirement and/or inadequate species occurrence.
Values are averages of times used in feasible runs. c: 2 to 41 runs were completed in 2 hours. Values are averages of times used in completed runs.
the number of equations (see Table 4). This indicates that model size may not be the only factor determining the computational efficiency of MIP models. Finally, we note that the Ö-CM model (which includes an additional distance component in the objective function, see (7) Table 3), while the Ö-CM model used only 1135.4 seconds (see Table 4). This, again, indicates that the model size is not the only factor that may affect computational efficiency of MIP models. We discuss these factors in the next section.

Discussion
Based on the computational experiments we conducted, we identified four important factors that affect the computational efficiency of MIP models in nature reserve design problems with contiguity and compactness considerations: 1) Model size, 2) The share of selected sites, 3) Model structure and 4) Input data. We discuss these factors below.
Model size is the most important factor affecting the computational efficiency of these models. As the model size (especially the number of binary variables) increases, solution time may increase substantially. Figure 1 depicts an example using the Williams (2002) model. As the number of candidate sites n increased from 100 to 1500, which gradually increased the number of variables from 6,091 to 18,831 and the number of constraints from 5,772 to 13,142, the solution time increased exponentially from 3 seconds to nearly 400 seconds. This is largely due to the number of nodes solved during the branch and bound procedure used by most MIP solvers. In general more nodes have to be solved when the model has a larger size.
Our test runs showed that the solution times also depend significantly on the share of selected sites in the optimal solution, which in turn depends on the specification of the population coefficient ( α ). This can be seen clearly, for instance, in the results obtained with the Shirabe (2005), Conrad et al. (2012) and Jafari and Hearne (2013) models for n=100 (see Table 2) and in the Williams (2002) model results for n=200 (see Table 3). Generally, these models follow a fast-slow-fast pattern, that is, as the share of selected sites increases, the solution time increases first, reaches a peak (where the share of selected sites is about 20-30 percent of all the candidate sites) and then decreases. This may be explained by the number of binary variables that have to be fixed at the values of 0 or 1 during the solution process. When the share of selected sites is small (or large), a large portion of binary variables will be fixed at 0 (or 1), leaving fewer combinations (a smaller branch and bound tree) that need to be searched for improving the objective function value (A brief description of the branch and bound algorithm can be found in Önal 2003). Similar results have been reported earlier by Williams (2002).
The Önal et al. (2016) model seems to be an exception. As noted above, it showed a relatively steady pattern in solution times for all population coefficient settings. The solution time tends to decrease only when a very small number of sites are to be selected.
Model structure is the third important factor that affects computational performance of the models we tested. Although the number of variables and constraints in a model may be substantial, due to the algebraic model structure, a large number of those variables and constraints may actually be redundant. Those redundant elements are eliminated during the preliminary solution procedure (Presolve in some GAMS solvers such as CPLEX and GUROBI) and the model sizes are much smaller after the Presolve. This is one of the main reasons for the superior computational advantage of the Önal et al. (2016)  In addition, some algebraic forms may be integer friendly and some may not. For example, the graph-theoretic contiguity model introduced by Önal and Briers (2006) has the following constraint: (10) where X ij and U j are binary variables as defined before and N j is the set of j's neighbour sites. This constraint states that if site j is selected, then at most four arcs can direct to it from its neighbour sites. The computational efficiency of the model can be improved substantially by modifying this constraint as in constraints (3) and (4). A comparison of the original Önal and Briers (2006) model and the modified model is depicted in Figure 2.
It is clear that replacing constraint (10) with constraints (3) and (4) slashes the solution time. We performed a paired samples t-test for each value of α to test the significance of the difference between these two models' processing times. The result reported a p value less than 0.05 for each α value, indicating that the difference between the processing times is statistically significant. This result should be attributed to the fact that a constraint like (10) is integer unfriendly while the constraints (3) and (4) improve the model's computational efficiency by helping construct a unimodular constraint matrix (ReVelle 1993, Williams 2002. Indeed, having a unimodular constraint matrix is one of the main reasons for the Williams (2002) model's good efficiency.
Input data such as land price and species distribution across candidate sites are also crucial factors that can affect the model's computational efficiency. Figure 3 depicts the processing times used by the models of Önal et al. (2016) and Williams (2002) when solving a 100-site problem with population coefficient α =0.5. The problem size remained unchanged across the 50 runs while species population and site costs in each run were varied. indicate that alterations in the input data may increase or decrease the solution time substantially even if the model size remains the same. This can also be explained by the workings of the branch and bound procedure. One data set may yield a 'good' integer solution in an early iteration and eliminate a large portion of the branch and bound tree; therefore, a fewer number of nodes needs to be solved (or just the opposite may occur). This seems to have happened in runs 4, 17 We should note that none of these factors alone can determine the performance of a MIP model, instead they affect the model's performance together. Solution techniques have been designed to speed up the computational process, such as supplying a cut value (Önal 2003, Cerdeira et al. 2005, Carvajal et al. 2013 or using a hybrid two-phase solution method (Conrad et al. 2012). We did not explore the effects of these techniques on the reserve design models tested here.

Conclusions
We have compared the computational efficiencies of eight representative optimal reserve design models that explicitly incorporate contiguity and compactness considerations. We identified the models presented by Williams (2002) and Önal et al. (2016) as the most efficient models. They both solved a 100-site contiguity problem easily, but encountered difficulty in solving the 200-site and larger contiguity problems. When compactness is also included, the Williams (2002) (2002) model when both contiguity and compactness are considered, but still the model could not be solved within 1 hour when 500 sites are involved. We identified model size, share of selected sites, model structure and input data as the four important factors affecting the computational efficiency of spatiallyexplicit optimal reserve design models. When developing such a model, one should try to define as few variables and equations as possible and to construct the model in such a way that its constraint matrix is unimodular (ReVelle 1993). These findings may be useful to conservation planners who select and use models in reserve design or management practices. This work may also provide guidance to academic researchers who try to develop computationally efficient optimal reserve design models.
The current nature reserves, whether terrestrial or marine, may not be adequate in providing sufficient protection to biodiversity. For instance, only over 1.89% of the marine protected areas (MPAs) worldwide is covered by exclusively no-take MPAs, far from the Convention on Biological Diversity's (CBD) Aichi Target 11 of 10% MPA coverage by 2020 and even further from the 30% target recommended by the IUCN World Parks Congress 2014 (IUCN 2018). It is almost certain that the current systems of nature reserves should be expanded. Acquiring more protected areas, however, would have to be done in a world where the human economy struggles to grow and the ecosystems and nonhuman species demands protection. Indeed, there is a conflict between economic growth and wildlife conservation (Czech 2014). Various disciplines such as ecological economics and environmental philosophy discussed this conflict and offered principles and guidelines to mediate the conflict and achieve sustainability. The science of nature reserve design or, more specifically, the optimal reserve design models, can contribute too. These models aim to solve practical conservation problems by incorporating ecological and socioeconomic considerations and help best prioritise conservation effort by guaranteeing optimal solutions. With the remarkable improvements in optimisation methods and computational power, it is reasonable to anticipate that spatial optimisation models will likely play a more significant role in the future than it does today in analysing and acquiring new protected areas.

Introduction
Freshwater habitats cover less than 1% of the earth's surface, but they support almost 10% of the known species on the planet (Strayer and Dudgeon 2010). However, freshwater species are in a worldwide decline (Howard et al. 2015), resulting thus in the decrease of freshwater biodiversity, which is mainly threatened by over-exploitation, water pollution, hydrological alterations, destruction and degradation of habitats and invasion by exotic species (Dudgeon et al. 2006). The Balkan region is rich in freshwater fauna due to high endemism (Griffith et al. 2004, Glöer et al. 2007). The rich fauna of freshwater molluscs is considerably threatened (Albrecht and Wilke 2008) and 29 species of the Balkan Peninsula have already become extinct (Régnier et al. 2009), representing 21% of the 140 known extinctions of freshwater mollusc species worldwide. In this region, many endemic gastropod species exist with a restricted range only extending to small hydrographic systems (rivers, lakes and springs). The karst landscape, which occurs in a large proportion of the Balkan Peninsula, also contributes to this high degree of endemism (Régnier et al. 2009). These species are vulnerable not only due to habitat degradation caused by human activities, but also by their very limited range, which increases the risk of accidental extinction (Régnier et al. 2009). For Greece, a similar crucial decline of population densities and losses of some endemic mollusc species have been reported (Albrecht et al. 2006), suggesting the necessity for urgent actions for their conservation.
The distribution of gastropod species is often associated with biotic and abiotic habitat parameters (e.g. Covich 2010), which are related to specific geographical characteristics, physicochemical conditions, biological interactions and historical and random factors (Lodge et al. 1987). Thus, the distributional patterns of endemics can be often better explained based on records of key environmental parameters, which drive the structure of biological communities (Pérez-Quintero 2012). Such knowledge may also help to identify and protect specific regions for the conservation of biodiversity of little-known freshwater bodies (Abell et al. 2007), especially in the Mediterranean region which has been neglected.
Out of those, two pyrgulinid species, belonging to the genus Dianella Gude, 1913, are known to be present in lentic systems of Greece: Dianella schlickumi Schütt, 1962 andDianella thiesseana (Kobelt, 1878). The first species is known from Lake Amvrakia and the latter from lakes Trichonis and Lysimachia (Aitoloakarnania, western-central Greece) (Schütt 1962, Radoman 1973, Szarowska et al. 2005, Szarowska 2006). Both species are narrow endemic species, with major gaps in historical records, while no information is available about their abundance, population structure, life cycle or genetic diversity. In the IUCN Red List of Threatened Species, Dianella thiesseana appears under the status "Critically Endangered", while D. schlickumi has been declared as "Possibly Extinct" (IUCN 2017). Apart from Dianella, the aquatic gastropod communities of these lakes include several other endemic hydrobiids, such as Islamia graeca Radoman, 1973, Islamia trichoniana Radoman, 1979, Pseudoislamia balcanica Radoman, 1979, Trichonia trichonica Radoman, 1973(Radoman 1983. D. schlickumi has not been observed for more than 30 years, despite several sampling efforts performed in Lake Amvrakia (Reischütz and Reischütz 2002, Szarowska et al. 2005, Albrecht et al. 2006, 2011a, which is considered as the type locality of this species. Only fresh empty shells of this species were recorded in an irrigation channel, near Acheloos River (close to Lake Amvrakia) by Reischütz et al. (2008). D. thiesseana has only been found at the northeast bank of Lake Trichonis during the last decade (Albrecht et al. 2009(Albrecht et al. , 2011b. The external morphology, the anatomy and the phylogenetic position of the majority of Pyrgulinae have been extensively studied and discussed (Schütt 1962, Radoman 1983, Riedel et al. 2001, Szarowska et al. 2005, Szarowska 2006, Anistratenko 2008, Wilke et al. 2013. The combination of anatomical and phylogenetic data suggests that the Pyrgulinae subfamily has a very close (possibly even sistergroup) systematic relationship with Hydrobiinae (Szarowska et al. 2005). Concerning the two species of Dianella, it has been reported that they have similar external morphology and soft body anatomy, but different body size (Radoman 1983). However, the inter-specific differences have been questioned by Szarowska et al. (2005), which has made the taxonomic distinction between D. thiesseana and D. schlickumi unclear and questionable. Although the soft body anatomy of D. thiesseana is described and/or depicted in detail (Radoman 1983, Szarowska et al. 2005, Szarowska 2006), various anatomical characters of D. schlickumi remain unknown, un-described or unpublished.
Hence, the aims of the present study were to: (a) report the rediscovery of the seemingly extinct pyrgulinid D. schlickumi, (b) describe its unknown anatomical characters, (c) examine its morphometric discrimination from D. thiesseana and (d) reveal the environmental parameters that probably drive the distribution of both species.

Study area
We studied three natural lakes located in the western-central part of Greece (Figure 1). These lakes are considered as the last remains of a former big lake, which by several processes, has been split into four smaller lakes (Albanakis et al. 1995). The Acheloos River separates Amvrakia and Ozeros lakes from Trichonis and Lysimachia lakes (Figure 1). Lake Amvrakia (Table 1) is a semipolje karst lake of tectonic origin with high sulphate concentrations, situated in Mesozoic limestone (Verginis and Leontaris 1978, Overbeck et al. 1982). Its water level shows significant seasonal fluctuations, which may be attributed to various factors such as the hydraulic communication with karst aquifers, high evaporation rates, especially during summer and intensive use of water for agricultural purposes (Danielidis et al. 1996). Lake Trichonis (Table 1) is the deepest and the largest natural lake in Greece, which is also a karst water body (Zacharias et al. 2002), situated in a tectonically active area (Poulimenos and Doutsos 1997). It is supplied by precipitation, approximately 30 intermittent streams and sub-aquatic karst springs . At its western end, a narrow canal connects Lake Trichonis with Lake Lysimachia (Table 1), which is shallower (Figure 1) and it was receiving the untreated urban wastewaters from the city of Agrinio until the year 2000. At its western end, Lake Lysimachia drains to the Acheloos River via an artificial canal (Figure 1).

Field sampling
Macroinvertebrate samplings were conducted in spring and autumn 2014 by boat, at the sublittoral and profundal zones of each lake, using an Ekman-Birge grab (three replicates per station, 225 cm 2 sampling area). A total of 66 samples were collected in the studied lakes [Amvrakia: 14 stations (7 stations at each zone), Trichonis: 13 (5 and 7 stations at the sublittoral and profundal zones, respectively) and Lysimachia: 6 (3 stations at each zone); Figure 1]. Samples were sieved with a 500 μm mesh and preserved  Danielidis et al. (1996) in 70% ethanol. After sorting, zoobenthos was identified to the lowest possible taxon and abundance was converted to density (ind./m 2 ). Additionally, at each station, water samples from 1 m above the bottom of each lake were collected, using a Niskin-type sampler. Environmental parameters such as water temperature (WT, 0 C), pH, conduc- tivity (Cond, μS/cm) and dissolved oxygen concentration (DO, mg/l) were measured in situ just above the bottom sediment of each station using the Aqua Read AP-2000 probe (Kent, GB). Secchi depth (Secchi, cm) was also recorded.

Anatomy
Genitalia of male and female specimens and other soft body characteristics of D. schlickumi from Lake Amvrakia were studied and compared with those of D. thiesseana collected from Lake Trichonis during the present study and with some additional individuals collected in 2016 (Figure 2). Before dissection, the shells of the specimens studied were removed by soaking in Perenyi solution, following Hershler and Ponder's (1998) morphological terminology. The characterisation of the nervous system was based on Davis et al. (1986) and its concentration was measured as the RPG ratio (Davis et al. 1976).

Morphometry
Shell morphometry was recorded in both D. thiesseana (Trichonis: 33 individuals and Lysimachia: 5 individuals) and D. schlickumi (Amvrakia: 44 individuals) in well preserved shells of adult specimens, using the software ImageFocus v3.0.0.1. Seven linear measurements were taken, specifically: shell height (H), shell width (W), shell aperture height (Ha), shell aperture width (Wa), spire height (SH), body whorl height (BWH) and penultimate whorl height (PWH) (Figure 3). Measurements were then expressed as ratios: shell height/shell width (H/W); shell aperture height/shell aperture width (Ha/Wa); spire height/body whorl height (SH/BWH), penultimate whorl height/body whorl height (PWH/BWH) and body whorl height/shell width (BWH/W). Prior to further analyses, all the linear measurements and the ratios were logarithmically transformed for normalisation. Subsequently, the t-test for independent samples was performed for detecting significant differences (p < 0.05) between the two species using the statistical software SPSS v22.
The studied populations were also tested for morphological diversification using a landmark-based method of acquiring geometrical data of shape and size. Shell variation has been traditionally quantified through linear measurements and ratios to distinguish between individuals and populations, amongst and within snail species. Recently, geometric morphometrics (GMs) have been employed for examining shells, both to provide direct size-free analyses of shell shape (Conde-Padin et al. 2007b, Hayes et al. 2007) and to answer broader evolutionary questions (Pfenninger and Magnin 2001, Conde-Padin et al. 2007a, Haase and Misof 2009, Páll-Gergely et al. 2012, Giokas et al. 2014. GMs have the advantage over traditional morphometrics in being, theoretically, free of effects due to size, position, rotation and scale (Rohlf and Marcus 1993).
For that purpose, we used sub-samples of the undamaged and well preserved shells of adult specimens (i.e. 15 D. schlickumi and 9 D. thiesseana specimens), using the aperture features. Specimens were set down on the same plane, with the aperture facing up and we took digital photographs of them. Geometric morphometric variables of the shells were obtained with 18 landmarks (LM) representing the outline of the shell and of the aperture as shown in Figure 4. Landmarks were digitised using tpsDig (Rohlf 2006) and occupied the same positions over the sum of the specimens. Two of these (LM1 and LM18) were of type I (developmentally homologous), 12 landmarks (LM2-LM13) represented suture lines and were considered type II (geometrically homologous). The remaining landmarks were either extreme points on surfaces or were placed halfway between other landmarks, making them type III (Bookstein 1991, Zelditch et al. 2004).
All geometric morphometric analyses were performed with MorphoJ (Klingenberg 2011). The coordinates of shape to be used for further statistical analysis were obtained with Procrustes generalised least squares superimposition. This method excludes the impact of size on the shape of the shell and, theoretically, variations independent of shape are removed using this analysis (Zelditch et al. 2004). Size variables   (centroid size) for each specimen were also generated. We examined, using pooled within-group regression analysis, the effect of centroid size on shape, regressing to the resulting Procrustes coordinates on centroid size. That regression score was not significant (Permutation test, 10,000 rounds against the null hypothesis of independence, % predicted = 6.170, p = 0.2183). Therefore, we eventually used the Procrustes coordinates, which were analysed with Principal Component Analysis (PCA) and Canonical Variate Analysis (CVA) in order to explore the morphological variation between the two Dianella species. The reliability of the discrimination was assessed by leave-oneout cross-validation.

Analyses of abundance and environmental data
Detrended Correspondence Analysis (DCA) (indirect gradient analysis) (programme CANOCO version 4.5.1; Ter Braak and Šmilauer 1998) was performed to abundance data in order to verify if they corresponded linearly to the environmental gradient or if abundances peaked around an environmental optimum (unimodal response) (Ter Braak and Šmilauer 1998). In this study, the species abundance data corresponded roughly linearly to the gradient, because the length of the gradient of the first axis was three times less than the range of the within-sample standard deviation (Ter Braak and Šmilauer 1998). Thus, Redundancy Analysis (RDA) was further applied to link the environmental parameters to species abundance data. The Monte Carlo permutation test was performed to select the statistically significant (p <0.05) environmental parameters that explain most of the variance of species abundance, while the variation inflation factor (>20) was used to assess the interrelation of the parameters.

Species distribution
Dianella schlickumi was found at the south-east part of Lake Amvrakia, in 5 (S7, S11, S12, S13, S14; Figure 1) out of the 14 stations (36%) surveyed, at depths ranging from 5 m to 13 m. No individuals were recorded at its western part (S9 and S10, Figure 1), not even empty shells. The highest abundance was recorded at station S15 (607 ind./ m 2 ) at depth 5 m (Table 2). D. thiesseana was found in 2 (S6 and S13) out of the 13 stations (15%) in Lake Trichonis, located on the northeast side of the lake, at depths of 13 m to 18 m ( Figure 1, Table 2). During the current study, only empty shells from D. thiesseana were collected in Lake Lysimachia. However, the presence of a live D. thiesseana was confirmed in the east side of Lake Lysimachia by surveys executed during the National Water Monitoring Programme in spring 2015 (Mavromati, personal communication). Both species were present during the spring surveys while none was found during autumn.

Soft body anatomy of D. schlickumi
Ctenidium-Osphradium: Ctenidial filaments broader than high; osphradium elongate, approximately opposite to the middle of ctenidium. Nervous system ( Figure 5A): Cerebral ganglia similarly sized, white; supraoesophageal and suboesophageal ganglia white; supraoesophageal connective much longer than the suboesophageal. Nervous system extremely elongate, RPG ratio 0.75. Gastric caecum: Large and elongate gastric caecum on the posterior stomach chamber. Rectum: The U-shaped intestine loop in the pallial cavity roof was wide and the faecal pellets were packed sideways ( Figure 5B). Female reproductive system ( Figure 6A): Pallial oviduct glands with straight border; bursa copulatrix medium-sized, ovoid-globular, positioned posteriorly to the albumen gland; bursal duct long; renal oviduct unpigmented, coiled tightly with more than two bends; typical seminal receptacles absent; a beak-shaped dilatation of oviduct in the position of seminal receptacle rs2, anteriorly to the bursa. Penis ( Figure 6B): Penis large, gradually tapered, proximal part folded, penial base of medium width, its attachment area central and well behind the head, penial duct almost straight located at the outer edge of the penis. Egg capsule ( Figure 6C-D): Several lenticular and transparent egg capsules (each one containing a single egg) in various growth stages were found on the external shell surface of two specimens.

Morphometry
Seven out of the 12 morphometric measurements taken, differed statistically (p < 0.001) between the two species. Generally, individuals of D. thiesseana exhibited higher mean values compared to D. schlickumi, except for SH/BWH ratio, which was higher in D. schlickumi (Table 3).  Geometric morphometric analysis revealed a clear size and shape distinction between the two Dianella species. Shells of D. thiesseana were significantly larger (Centroid Size) than those of D. schlickumi (F 1, 22 = 38.01, p < 0.0001, Figure 8). The landmark-based analyses of individual shells also confirmed a shape distinction between the two species. ANOVA of Procrustes coordinates also showed a clear shape difference between the two species (F 32, 704 = 20.78, p < 0.0001). PCA of Procrustes coordinates, which aimed at finding linear combinations maximising the total variance, revealed three principal components (PC1-PC3) explaining 88.1% of the total variance. The distinction between the two species was quite evident along the PC1 axis (Figure 9). In the CVA, which maximised differentiation amongst species (or predefined groups), the first and only extracted canonical variate  (CV1) explained 100% of total variation and showed a much clearer separation than PCA between the two morphotypes along the CV axis ( Figure 10). Clear shape discrimination was also found between D. thiesseana and D. schlickumi shells after a permutation Discriminant Function Analysis (T-squared = 236.3713, p (parametric) = 0.8335, p (permutation) < 0.0001). However, after a permutation test with 1,000 replicates, only 54.2% of "unknown" individuals (leave-one-out cross-validation) were correctly classified. Shape changes along the CV1 are shown in Figure 11 and portray the tendency for D. thiesseana to be slimmer compared to D. schlickumi and also to have a larger, i.e. elongated, aperture.

Environmental parameters
The Monte Carlo test in RDA analysis revealed two out of the six environmental parameters examined as significantly different between the two species (p < 0.05): conductivity and Secchi depth. Depth and water temperature were excluded from the analysis   due to their high inflation factor (> 20). The first and all canonical axes were statistically significant (p = 0.044 and p = 0.002 respectively). The first two ordination axes of RDA explained 99.5% of the total species variance (Monte Carlo test, p < 0.05). Axis I (eigenvalue 0.929, p < 0.05) was related to conductivity and DO (intra-set correlation values 0.952 and -0.788 respectively) and Axis II (eigenvalue 0.066) to pH (intra-set correlation value 0.428). D. schlickumi specimens from Lake Amvrakia were ordered at the positive side of Axis I, having the highest conductivity values while D. thiesseana specimens from the other lakes were ordered at the negative side of Axis I and correlated with lower values of conductivity and higher values of DO (Figure 7, Table 2).

Discussion
Dianella schlickumi was recorded in the sublittoral zone of Lake Amvrakia in 1962 (Schütt 1962), but no exact details about this location were provided. The latest record confirming the presence of the species in the lake was before 1983 (Radoman 1983). After that, D. schlickumi was assumed to be extinct, as reported by Reischütz and Reischütz (2002), Szarowska et al. (2005) and Albrecht et al. (2006), based on sampling efforts in 1983, 2002 and 2003 respectively. Albrecht et al. (2011a) mentioned that the originally recorded area of D. schlickumi is the northern extent of Lake Amvrakia. At this part of the lake water, level fluctuations were high during the last years, due to insufficient connectivity with the main body of the lake. In fact, during extremely dry summers, this area is completely dry or water remains only in small pools. In such cases, the predation pressure (mainly from fish) is intense, resulting in local extinction of molluscs (Nilsson et al. 2008, Downing et al. 2010. Such rough conditions have probably negatively affected the population size and the distribution of D. schlickumi in this lake. Nevertheless, in our study, we were able to collect specimens of various ages at several sampling sites, indicating that this snail is still extant and thriving in its type locality. Under the IUCN Red List guidelines (IUCN 2017), D. schlickumi is considered as "Critically Endangered (Possibly Extinct)" [following the criteria B1ab (i, iii)] (Albrecht et al. 2011a). However, since a viable population of this species was demonstrated in the present study, the tag "Possibly Extinct" should be removed from its assessment. The rediscovery of the species does not however improve its conservation status, as its extent of occurrence is estimated to be less than 100 km 2 , restricted to the sublittoral zone of Lake Amvrakia (<10 km 2 ). Moreover, D. schlickumi is known to exist only in a single location, where habitat degradation was also noticed due to agricultural practices and pesticide contamination (Thomatou et al. 2013 a, b). Consequently, the classification of Critically Endangered should be maintained, as the IUCN criteria continue to be met. As for D. thiesseana, it should also continue to be considered as Critically Endangered [B1ab (ii,iii) and 2ab (ii, iii)] since it has an extent of occurrence of 94 km 2 and an area of occupancy of 2 km 2 . Its habitat is also degraded because of intensive agricultural practices, urban sewage, stock grazing land and small industries (Bertahas et al. 2006).
Both species thrive on soft substrate (Albrecht et al. 2011a, b). D. schlickumi seemed to prefer slightly shallower habitats, from 3 m to 15 m, while D. thiesseana had a range of distribution from 2.5 m to 27 m. The species absence from autumn surveys may be due to depth migrations often observed in gastropods (Brown 2001). Even though this study has confirmed the rediscovery of D. schlickumi, it does not clarify in detail the accurate distribution of the species, which might be the subject of future research.
The anatomical characters of D. schlickumi known from literature (Schütt 1962) (i.e. the radula) and those described in the present paper (i.e. the ctenidium and osphradium, the central nervous system, the gastric caecum, the rectum and the male and female reproductive organs) are similar to those reported by Radoman (1983) in the diagnosis of the family Pyrgulidae Brusina 1881 (nominal subfamily Pyrgulinae) and of the genus Dianella Gude, 1913[(type species Dianella thiesseana (Kobelt, 1878]. Similarities have also been found between D. schlickumi and D. thiesseana in the nervous system, the gastric caecum and the female reproductive organs as they are described for the latter species by Szarowska et al. (2005, pp. 22-31 and related figures). However, a notable difference exists in the shape of the penis between the two species: the penis of D. thiesseana is wide, with a triangular shape due to its expanded base and has a clear outgrowth on its left side (Szarowska et al. 2005, p. 29, figs 32-37). The penis of D. schlickumi does not have an expanded base and is gradually tapering without any distinct outgrowth.
D. schlickumi and D. thiesseana are discriminated due to their shell dimensions; the first one is smaller than the second (Radoman 1983). In our study, traditional and geometric morphometric analyses support the assumption that D. thiesseana and D. schlickumi are quite distinguishable in terms of external morphology. However, neither in the original description of the two taxa (Kobelt 1878 andSchütt 1962 respectively) nor in the subsequent studies (Radoman 1983, Szarowska et al. 2005, shell shape appears to be significantly different between the two taxa. Moreover, Szarowska et al. (2005) reported that the shell characters and dimensions of D. schlickumi, described in the literature and those of D. thiesseana from Lake Lysimachia, lie within the range observed for D. thiesseana from Lake Trichonis. Even though, the intra-subspecific morphometric variation is quite high, the two species are actually distinct in size and shape (Figures 8-10). Therefore, in the "Dianella" case, geometric morphometric analysis has been proved to be a useful tool for exploring patterns of species diversity and seems promising for examining morphological adaptations of these species. Specifically, it is clear that the Centroid Size, derived from Geometric Morphometric Analysis, can be used safely to distinguish those two species.
In some cases, the environmental parameters set limits (Pyron and Brown 2015) and may influence (Økland 1983, Pip 1986, Falniowski 1987) the distribution of freshwater snails. The analysis performed revealed conductivity, DO and pH as the most significant parameters explaining the current species variance. D. thiesseana, according to RDA, seems to prefer more oxygenated waters with lower conductivity values and pH than D. schlickumi. Conductivity, dissolved oxygen (e.g. Çabuk et al. 2004, Kazibwe et al. 2006, Cloherty and Rachlin 2011 and pH (Økland 1983(Økland , Pip 1986(Økland , Çabuk et al. 2004) are also referenced as key environmental proxies explaining gastropods distribution in water bodies. Biological interactions and/or historical factors, which need further surveys and analyses, may also be responsible.
Both Dianella species require urgent conservation measures to be enforced and suitable management plans to be implemented in the whole studied area, focusing on the protection of these species and the improvement of their habitats. The studied lakes form an ecologically important complex as Special Conservation areas (Council Directive 92/43/EEC), listed in the NATURA 2000 network. However, the integrity of local mollusc populations in these ecosystems is threatened due to the presence of invasives; Ferrissia fragilis (Tryon, 1863) and Physa acuta Draparnaud, 1805 in lakes Trichonis and Lysimachia (Albrecht et al. 2014) and of Potamopyrgus antipodarum (J.E. Gray, 1843) in Lake Trichonis (Radea et al. 2008). The invasive impacts upon native molluscs may be direct, through competition for food and space (Carlsson et al. 2004) or indirect, through changes in ecosystem function and/or parasitism (Brown et al. 2008). The antagonistic behaviour of P. antipodarum and P. acuta is already wellknown (e.g. Schreiber et al. 2003, Zukowski andWalker 2009), although no such outcomes have been documented for F. fragilis until now (Albrecht et al. 2014). Available knowledge on the effects of invasive species on native mollusc species and especially on D. thiesseana, is insufficient to estimate the possible threats.
One of the main pressures identified in lakes Amvrakia, Trichonis and Lysimachia is eutrophication from agriculture run-off and wastewaters (Bertahas et al. 2006, Thomatou et al. 2013. Dianella inhabits the sublittoral zone of these lakes, which is mainly affected by eutrophication (Pilotto et al. 2012). The development of water quality indices based on zoobenthic communities and related to eutrophication will help to assess the ecological status of these water bodies and prevent their further deterioration. The impacts could be minimised by preserving natural vegetation or establishing buffer zones. Only low-impact activities should be allowed near the lakes and harmful practices should be prohibited or regulated to a more distant zone. Furthermore, it is crucial to restore the type locality of D. schlickumi by maintaining the connectivity of the northern extension of Lake Amvrakia with the main lake water body, thus preventing the species' historical distribution area to become dry. Additional measures to reduce water abstraction could help the maintenance of the ecological water level in the lakes, securing slight seasonal fluctuations and preventing aquatic biodiversity in general (Zohary and Ostrovsky 2011). Further temporal and spatial surveys which will identify the ecological features, the ecosystem functions and the habitat preferences of the studied species could significantly support conservation efforts. Environmental education and citizen awareness could additionally support conservation efforts for both endemic gastropods.

Conclusion
Our results confirm the rediscovery of the narrow endemic species Dianella schlickumi in Lake Amvrakia since past sampling efforts did not detect it for more than 30 years. Moreover, they provide basic knowledge on its anatomical, morphological and distributional patterns and its discrimination from Dianella thiesseana. The two species are discriminated by their shell size and shape. Moreover, the presence of Dianella thiesseana is confirmed in lakes Trichonis and Lysimachia. We conclude that further conservation measures have to be implemented for effective protection of both species.

Breeding site fidelity, and breeding pair infidelity Introduction
Carnaby's Cockatoo Calyptorhynchus latirostris is endemic to south-western Australia. Formally common throughout its range, over the second half of the last century, extensive clearing of native vegetation for the development of broadscale agriculture and urban development saw the species contract significantly in range and abundance. As a result of the loss of breeding and foraging habitat (Fig. 1), the status of the species changed from being classified legally up to the 1960s as vermin with a bounty on its bill, to being listed as endangered (Saunders 1990)  Beginning in 1968 and continuing, the ecology and behaviour of the species has been the focus of one of the longest running vertebrate studies in Australia (Saunders 1982, Saunders and Ingram 1998, Saunders and Dawson 2017. One breeding population of the species at Coomallo Creek ( Fig. 1) was studied intensively from 1969 to 1977, then monitored intermittently using the same protocols until 1996, and then every year from 2009 (Saunders and Dawson 2017) until the present. As a result of this study, the breeding biology of the species is well known. Once the birds form pairs and reach breeding age (usually at four years), the pairs remain together until one of the partners dies. They return to the same breeding area each year in winter, with rainfall in the first half of autumn influencing the commencement of egg-laying; the wetter the period, the earlier breeding commences. The female selects a nesting hollow in a large eucalypt, usually lays two eggs, on average up to eight days apart. Both eggs usually hatch with the younger nestling dying within 48 hours of hatching. In around six per cent of cases both nestlings are successfully fledged. The female nests in the same hollow used the previous breeding season, provided she successfully fledged a young the previous year, and the hollow is not already occupied when she returns to breed. If she has an unsuccessful breeding attempt she usually moves to another hollow nearby, and makes another breeding attempt, either during the same breeding season or next season (information from Saunders 1982, Saunders et al. 2013, 2014a. The fact that Carnaby's Cockatoo is rare, and difficult to breed in captivity (Saunders 1976) made it a lucrative target for poachers who took nestlings illegally from hollows used by wild breeding pairs. In accessing the hollow floor, poachers often cut a large hole in the side of the nesting tree rendering the hollow unsuitable for Carnaby's Cockatoo. In dealing with the problem posed by poachers, the Western Australian Government authority responsible for conservation, in conjunction with staff at Curtin University, investigated the efficacy of using analysis of DNA to establish relationships between putative parents and the nestlings poachers claimed they had bred in captivity. White et al. (2009) developed 20 microsatellite markers for population and forensic applications suitable for Carnaby's Cockatoo. White et al. (2012) demonstrated the successful application of these markers in three case studies, one of which demonstrated clearly that a Carnaby's Cockatoo nestling claimed to have been bred in captivity bore no relationship to the pair the breeder claimed had produced the nestling, and further that the nestling was a full sibling of a nestling known to have been fledged successfully from a nest hollow monitored in the wild. This unequivocal result was used in a successful prosecution for poaching. White et al. (2012) demonstrated that analyses of DNA provided reliable discriminatory power for identification of individuals, testing kinship and identification of breeding populations, making possible the determination of illegal trade or harvest of live animals from the wild, and providing a method to monitor both wild and aviary populations. They pointed out that use of DNA is a minimally-invasive technique, even less so when moulted feathers are used as described by Hou et al. (2018). White et al. (2012) also noted that the black cockatoo database they described constituted one of the most extensive Australian wildlife forensic databases.
In this paper we examine the use of DNA taken from nestlings of one breeding population over seven breeding seasons (2003,2005,(2009)(2010)(2011)(2012)(2013) to investigate aspects of the breeding biology of Carnaby's Cockatoo. We compare the results obtained from DNA analyses with earlier research  on the population using observations of individually marked birds to test the efficacy of analyses of nestling DNA for studies of breeding behaviour.

Methods
The Coomallo Creek (Fig. 1) study area is described in detail in Saunders (1982), Ingram (1998), andSaunders et al. (2015). The cockatoos nest in hollows in large eucalypt trees that are distributed along a 9km narrow belt of woodland surrounded by agricultural land and uncleared Kwongan (native shrub land sensu Beard 1984). There was no further clearing of the eucalypt woodland during the study period, but the number of available nest hollows fluctuated as trees were lost through natural attrition (Saunders et al. 2014b), and new hollows entered service, and were located by us.
From 1969 to 1996 the area was visited for 21 of those years for at least one week in early September, and one week in early November. During each visit every known nest hollow was inspected, and the identity of breeding females established if possible by trapping or by using a telescope to read the number on the bird's leg band/ring. Both methods were time consuming, and it was not possible to identify all banded females. Searches were also undertaken to find any nesting hollows in use not registered in the study area. Any nestlings older than three weeks were measured (length of the folded left wing and body mass) and banded with a stainless steel band with a unique number. Nestlings were aged using the methods of Saunders (1986), and laying dates extrapolated from the age of the nestlings.
The area was visited once in November 2003 and 2005, and each year from 2009 to 2013 the area was visited for one week in early September and early November, with a three-day visit in early January to establish nesting success, and band any nestlings laid late in the season. The identity of breeding females was established by photographing their legs, and checking for bands using the methods of Saunders et al. (2011). This was a much more successful method for establishing the identity of banded females nesting in the area than was possible in the earlier study. Searches were also undertaken to find any nesting hollows in use not registered in the study area. Nestlings were measured (length of the folded left wing and body mass) and banded, and laying dates extrapolated from nestling age. Nestlings were sexed on shape and colour of their cheek patch, and at least three pin feathers taken from their chests for subsequent DNA analysis. The body of small dead nestlings were taken for DNA analysis, as was at least one toe from larger dead nestlings or dead adults.
Material taken for DNA analysis was plucked with tweezers or cut with scissors, and placed in a vial of preservative (20% dimethylsulphoxide saturated with sodium chloride). The vial was labelled with species identification, collection date, name of collector, name of the study area, nest hollow number, and the description of material taken (e.g., 4 feathers, small dead nestling, toe of dead nestling). All equipment used for handling DNA was washed in disinfectant, 90% ethanol, and rinsed with water between samples to ensure no cross-contamination of DNA.
A total of 309 DNA samples were collected from Coomallo Creek between November 2003 and October 2013 for DNA extraction with Qiagen Blood and tissue kit (Qiagen). Samples were subjected to a real-time quantitative PCR (qPCR) assay to assess the DNA extracts for quality and quantity of nuclear DNA prior to microsatellite genotyping. Genetic profiles were generated from 16 microsatellite markers as previously published by White et al. (2009White et al. ( , 2012White et al. ( , 2014. The software COANCESTRY 1.0.1.7 (Wang 2011) was used to estimate the pairwise relatedness (r xy ) between individual birds. A total of 873 birds were included in the analysis, 587 birds were DNA profiled and published in White et al. (2014), and an additional 286 DNA samples were collected across the breeding range of the species for pairwise relatedness examination. COANCESTRY implements seven methods to estimate the pairwise relatedness between individuals. These methods are: the triadic likelihood estimator (denoted as TrioML; Wang 2007); five moment estimators denoted as Wang (Wang 2002); LynchLi (Lynch 1988, Li et al. 1993; LynchRd (Lynch and Ritland 1999); Ritland (Ritland 1996); QuellerGt (Queller and Goodnight 1989); and a dyadic likelihood estimator (DyadML, Milligan 2003). The relatedness classifications for the genealogical relations were set at r xy = 0.5 for parent-offspring and full siblings, r xy = 0.25 for half siblings, and r xy = 0 for unrelated individuals, with 1000 bootstraps to calculate 95% confidence intervals. For the purposes of this study, we report the genealogical relations of parent-offspring and full siblings from the 309 DNA samples collected at Coomallo Creek.
Distances between nest hollows used successively by the same breeding pair were calculated using the "show distance between waypoints" function of OziExplorer® software. The density of nest hollows in the study area was only calculated for the period 1972-1990 and 2009-2013 due to the fact that nest hollows were being located for the first time at a high rate during the first two years of the study (1970)(1971).

Results
Between 1970 and 1990, observations of individually marked females provided data on nest hollow use by 92 breeding pairs, each making two or more breeding attempts. Pairs recorded making two or three breeding attempts provided 68.5% of the data. During that period one female was recorded making 12 breeding attempts (Table 1). Twenty of these breeding pairs had both birds individually marked. In all but one case the cockatoos remained paired with same partner, one pair being together for five years  1970-1990 and 2009-2013. 2  3  4  5  6  7  8  9  10 11 12 Total   1970-90  45 18 15  7  2  1  1  0  2  0  1  92  2009-13 24 11 4 2 1 0 0 0 0 0 0 42  (1972)(1973)(1974)(1975)(1976)(1977), during which time they made seven breeding attempts. There was only one case of pair separation, where both partners were subsequently seen breeding with other partners. Between 2009 and 2013, DNA taken from nestlings provided data from 42 breeding pairs making two or more breeding attempts (Table 1). Most data (83.3%) were from females making two or three breeding attempts over the five year period, but one female made six attempts, one each year 2009-2013, and a second breeding attempt in 2013 after the first attempt failed. All attempts were with the same partner. One hundred and forty-two relationships between nestlings were established, ranging from two nestlings from the same breeding attempt to nestlings from six breeding attempts. Data were obtained from 172 successive breeding attempts from the 92 pairs between 1970 and 1990, and 63 successive breeding attempts from 42 pairs between 2009 and 2013 (Table 2). One hundred and thirty-three (93.7%) successful breeding attempts between 1970 and 1990 were followed next season by the pair breeding in the same hollow or a different one if the previous one was not available. Between 2009 and 2013, 51 (92.7%) successful breeding attempts were followed next season by breeding in the same hollow or a different one if the previous one was not available. Between 1970 and 1990, 28 (93.3%) unsuccessful breeding attempts were followed by breeding attempts in a different hollow. Between 2009 and 2013, seven (87.5%) unsuccessful breeding attempts were followed by breeding attempts in a different hollow (Table 2). There was no difference (χ 2 =1.21, d.f.= 4, p=0.88) between the relative rates of success or failure of breeding between the two phases of the study.

Number of breeding attempts made by individual females
During the period 1970-1990, data on movements to different hollows between successive breeding attempts were obtained from 60 females that moved 94 times, and 2009-2013 from 30 females that moved 40 times (Fig. 2). The mean distance pairs moved between successive breeding attempts was 0.66±0.74km (n=94; range 0.02-4.05km) for the period between 1970 and 1990 and 1.25±1.31km (n=40, range 0.07 -5.50km)  (Table 3). There was a significant difference in the mean distances between successive nest hollows for the periods 1970-1990 and 2009-2013 (Mann-Whitney U-test, one-tailed: Z=-2.81, p=0.002). The distance between nesting hollows used in successive years increased despite the average density of nest hollows remaining C similar between the two phases of the study (Table 4; Mann-Whitney U-test, one-tailed: Z=-1.29, p=0.10).
For the period 2009-2013, DNA was obtained from both nestlings from 22 breeding attempts. Of these, six (27.3%) of the sibling pairs had different males fertilising the two eggs. Both nestlings fledged successfully in 13 cases where both eggs were fertilised by the same male (13/16 = 81.3%), and in four cases where the eggs were fertilised by different males (4/6 = 66.7%). Of the five cases where one or both of the nestlings died, both eggs were fertilised by the same male in three (60%) cases, and two (40%) by different males. The ages of three of the adult females for which there were DNA data on both nestlings from a single breeding attempt were known as the females were banded. Band number 210-03089 was at least 22 (she was banded as an adult), and both eggs were fertilised by the same male; 210-01876 was at least 26, 27 and 28, and the same male fertilised both eggs in all three clutches; and 210-01694 was 20, and the two eggs were fertilised by different males.
The average time between the laying of the two eggs in each case of siblings for which laying dates were known was 10.6±3.9 days (n = 15; range 5-17 days) where the same male fertilised both eggs, compared with 14.8±5.9 days (n = 5; range 9-24 days) where the second egg was fertilised by a different male. There was no significant difference in the mean interval (in days) between the laying of successive eggs in nests fathered by the same male and different males (t-test: t=-1.71, d.f.=17, p=0.053).

Discussion
Nest site fidelity, breeding site tenacity or breeding philopatry is a well-known phenomenon, demonstrated by many vertebrate species (Greenwood 1980). Among birds it is practiced by both migratory and sedentary species, from a wide range of taxa, both passerines and non-passerines, and in most species it is females that display breeding site fidelity (Greenwood 1980). Advantages accruing as a result of breeding site fidelity are Table 4. Density of nest hollows (km -2 ) used by Carnaby 's Cockatoos 1972-2013 Mean±s.d. density of nest hollows (km-2) Range 1972Range -1990 14.0±4.9 4.5-23.8 2009-2013 16.9±4.8 11.2-23.8 reduction in susceptibility to predation (Burger 1982) and other adverse conditions, such as micro-climatic impacts, and knowledge of rich food resources (Blancher andRobertson 1985, Ehrlich et al. 1994). In such species, breeding success in one season is often followed by use of the same site the next season (Blancher and Robertson 1985). There is tendency for breeding site fidelity to increase with age of individuals, based on the fact that older birds tend to be more experienced breeders, and therefore successful, so moving sites less frequently (McNicholl 1975). McNicholl (1975) pointed out that breeding site fidelity is strongly developed in highly stable habitats (e.g., cliff faces), but greatly reduced in highly unstable habitats (e.g., grasslands, mud beaches or bars in rivers). Carnaby's Cockatoo breeds in highly stable habitats (sensu McNicholl), in that tree hollows are usually available for many years (Saunders et al. 2014b). So it is unsurprising that analysis of successive breeding attempts based on observations of individually marked females at Coomallo Creek between 1970 and 1990 confirmed breeding site fidelity, with 93.7% of successful breeding attempts being followed the next season by the birds breeding in the same hollow, or a hollow nearby if the hollow was in use when they returned to breed. Analysis of nestling DNA found the same result for the period between 2009 and 2013, with 92.7% of successful breeding attempts being followed the next season by the birds breeding in the same hollow, or a hollow nearby if the hollow was in use when they returned to breed. Such breeding site fidelity is shown by other species of Psittaciformes (Rowley 1990, Murphy et al. 2003, Berkunsky and Reboreda 2009. Analyses of data on hollow use following breeding failure from observations of individually marked females and nestling DNA confirmed Saunders' (1982) finding that after a failed breeding attempt, females usually moved to another hollow nearby for their next breeding attempt, either in the same breeding season or the next season. Saunders (1982) established that when Carnaby's Cockatoo formed breeding pairs, in most cases they were monogamous, pairing for the life of one of the partners. This was based on observations of birds individually marked with patagial tags. Not only did they remain paired during the breeding season, but they remained together throughout the year (Saunders 1980). Analysis of nestling DNA confirmed this fidelity of breeding pairs from 2009 to 2013, but it also revealed extra-pair fertilizations. In 27.3% of the 22 cases where DNA samples were obtained from both nestlings from the same breeding attempt, the second egg was fertilised by a male not paired with the female. While these extra-pair fertilizations were surprising in light of observations of individually marked breeding Carnaby's Cockatoo, this extra-pair fertilization is not uncommon in birds (Petrie andKempenaers 1998, Forstmeier et al. 2011), and mammals (Crawford et al. 2008). Extra-pair fertilization has been recorded at varying rates in some larger parrot and cockatoo species, for example: 9% extra-pair fertilization in nests of Blue and Yellow Macaw Ara ararauna (Caparroz et al. 2011); 20% in nests of Red-tailed Amazona Amazona brasiliensis (Fernandes 2015); and 40% in nests of Monk Parakeet Myiopsitta monachus (Martinez et al. 2013). In parrot species with multi-male breeding systems, extra-pair fertilizations have been suggested in the Echo Parakeet Psittacula eques (Taylor and Parkin 2009), confirmed in Eclectus Parrot Eclectus roratus (Heinsohn et al. 2007), and confirmed in 100% of nests of Vasa Parrot Caracopsis vasa (Ek-strom et al. 2007). However, in other species such as the Burrowing Parrot Cyanoliseus patagonus (Masello et al. 2002) and Crimson Rosella (Platycercus elegans; Eastwood et al. 2018) there is no evidence of extra-pair fertilization.
The chances of observing extra-pair copulations of individually marked birds were low. While many breeding females were individually marked with patagial tags and/or leg bands, far fewer males were individually marked. While copulations were observed, those few cases where both sexes were identified were all with known mates, reinforcing the conclusion of monogamy. In 2010 a documentary about the breeding of Carnaby's Cockatoo was filmed (http://www.abc.net.au/tv/guide/abc1/201203/programs/ NH1001W001D2012-03-13T203229.htm accessed 5 April 2018). The documentary was based on one breeding pair, and after the first egg was laid, film was taken of the female leaving the hollow when her unbanded mate was not around. There were two adult males nearby, and both immediately started courtship displays towards the female. One of the males had a leg band, and that could have only been placed on the bird seven years previously or longer. Earlier observations and this film clearly indicate that males will readily court breeding females other than their mates (Saunders 1979), and analysis of nestling DNA indicates that females will accept copulation from males other than their mate. It is highly unlikely that the first egg would result from copulation with a male other than the mate. This is because the pair remain together from the time they return to the breeding area, select and prepare their nesting hollow, and the first egg is laid. At that time the female incubates the egg, and depends on the male for food (Saunders 1982), and males are aggressive to other males approaching their mate (Saunders 1979). During incubation, and until the older nestling is about three weeks old, the female remains either in the hollow or around the nest tree. The male is around the nest tree at dawn, mid-morning when he feeds the female, and again when he feeds her at dusk (Saunders 1982). On occasions the female may be out of the hollow when the mate is not around. At this time extra-pair copulations are possible. The second egg is laid on average eight days after the first, but may be laid up to 24 days later (Saunders 1982, Saunders et al. 2014a. At Coomallo Creek, both nestlings have fledged from hollows when there was 24 days difference in age, with nestlings growing at the same rate; i.e., the younger nestling does not grow faster than the older in order to fledge at the same time. These differences in time between the eggs provide ample opportunity for extra-pair copulations. The limited data available from analysis of nestling DNA indicate that older, more experienced females may indulge in extra-pair copulations, and that second nestlings, the result of such copulations may fledge along with the older nestling, the result of copulation with the male of the pair. These data came from banded females, however nothing is known about such rates in younger breeding females. Who were the males responsible for these extra-pair copulations? Of the six cases known to us, in three cases, there were no data indicating the second nestling was a half-sibling to any other nestlings from the area that breeding season for whom we had DNA. In each of the other three cases (one in 2009, 2010 and 2012), there was a nestling from three other breeding pairs that indicated a half-sibling relationship. The average distance between any of these hollows and the hollow with the extra-pair copulations was 2.8±1.2km (range 1.7-5.2km). If it was one of the males producing the half-siblings, they weren't neighbours. Without DNA from the males, unfortunately we cannot identify any of the males responsible.
The conservation implications for such extra-pair fertilizations in Carnaby's Cockatoo are difficult to discern. Extra-pair fertilizations within local breeding populations would result in broader genetic diversity within breeding populations than would result from purely monogamous pairings.

Conclusions
Comparisons of results obtained from long-term studies based observations of individually marked Carnaby's Cockatoo breeding adults, and those from analyses of DNA from Carnaby's Cockatoo nestlings showed that analyses of DNA are useful for establishing stability of breeding pairs and breeding site fidelity. Analyses of nestling DNA demonstrated that aspects of mating behaviour that would be extremely difficult to show using observations of individually marked adults could be revealed by studies of DNA. Analyses of DNA from nestlings over subsequent breeding seasons would be necessary to investigate the identity and ages of the adults involved in extra-pair copulations, and provide evidence of the advantages to the species of breeding pair infidelity. White  Stephanopachys linearis (Kugelann, 1792) (Coleoptera, Bostrichidae) in Poland

Introduction
Stephanopachys linearis (Kugelann, 1792) is a representative of the forest fauna characteristic of boreal and alpine regions of the Palaearctic. It belongs to the family of horned powderpost beetles (Bostrichidae), of which 29 native species occur in Europe (Borowski 2007), 9 of them -including S. linearis -having been reported from Poland (Dominik 1958, Burakowski et al. 1986). S. linearis was discovered and described in 1792 by the German entomologist J.G. Kugelann, working as an apothecary in Ostróda in the (then) Olsztyn regency of East Prussia (now NE Poland). Besides short Latin and German descriptions of Apate linearis, Kugelann gave a diagnosis allowing its distinction it from the now common bark beetle, Anisandrus dispar (Fabricius, 1792) (Curculionidae, Scolytinae), with information that the new species had been found only once, on an old hoarding. Kugelann's publication has become the basis of the inclusion -after World War II -of S. linearis in the Polish fauna (Dominik 1958, Burakowski et al. 1986, Borowski 2007, but during the 220 years since its description, it has never been reported from this country again and, for example, the renowned German faunist Adolf Horion (1961) considered its occurrence in Poland or Baltic countries as doubtful.
The current area of the distribution of S. linearis extends over the boreal zone of the Palaearctic, from Norway and Denmark, through Sweden, Finland and Siberia to the Far East (it was recently discovered in northern Manchuria - Zhang et al. 1995), reachingthrough Estonia, Latvia and Lithuania -as far south as NE Poland, Belarus and Ukraine; relict localities are also spread over the mountainous regions of central and southern Europe (France, Switzerland, Germany, Austria, Italy, Slovenia and Bohemia -Borowski 2007), as well as the Caucasus (Horion 1961, Geis 2002) and on Corsica (Sainte-Claire Deville 1902). Its recent discovery in northern Iran (Liu et al. 2016, Nardi andAudisio 2016) seems highly doubtful: the locality, albeit close to the border of Azerbaijan, is practically devoid of vegetation, with the climate unsuitable for the development of this species, so the specimen seems either misidentified or mislabelled; one of us (JB) contacted the Iranian author of the respective publication for possible verification of its taxonomic identity, but the specimen proved unavailable not only for loan but even for taking a photograph.
In the years 2009-2017, in the course of the study of saproxylic insects on a burned site in Białowieża Primaeval Forest, the existence of this species has been discovered. Information concerning this discovery and the data supplementing the current knowledge on the ecology of the development of S. linearis are provided below.

Locality and collection data
In 2009, the research programme was launched in Białowieża Primaeval Forest (Faliński 1986) to study the dynamics of the changes in species composition and abundance of (especially saproxylic) beetles after the disturbance caused by the forest brushwood fire on 28.04.2009. The studied area of the burned site (ca. 7 ha), in the compartment 105B of the Białowieża National Park, included natural stands with high proportion of old (120 years or more) trees in humid forest (60.5% of the disturbed area), boggy forest (24.5%) and humid mixed forest (13.2%). The site has been left without economic activity, with actions having been restricted to scientific research. The observations started a few days after the fire. Beetles were caught in 2009-2011 and 2015-2017 using 9 non-baited Moericke's and 18 Netocia-type traps (Piętka and Borowski 2015) as well as 3 black funnel traps each baited with α-pinen, ethanol, fuscumol, ipsenol and ketol (https://www2.gov.bc.ca/assets/gov/farming-natural-resources-and-industry/forestry/ forest-health/forest-health-docs/spruce-beetle-docs/spruce_beetle_funnel__traps.pdf ). All traps were suspended 1-2 metres above ground.
In June 2017, we measured the diameter at breast height and extent of trunk scorching (from the ground level) of 140 living pine-trees in the burned area and evaluated (in a three-level scale: slight, medium, strong) the degree of insolation of each trunk. All trees had been carefully examined for the presence of S. linearis larval galleries.

Statistical methods
To assess the preferences of Stephanopachys linearis for particular types of trees, we used a Generalised Linear Model (GLM), with binomial distribution for response variable (0 = tree not occupied by beetle, 1 = tree occupied by beetle). In the model, we used three explanatory variables: diameter of a tree degree of burning of the trunk (continuous, log-transformed variables) and exposure to the sun (categorical, three-level variable). All analyses were done using R (version 3.2.2) software.

Results
The first specimen of Stephanopachys linearis (Fig. 1) was caught between 26.04.2015 and 11.05.2015 in the Netocia-type trap on living but scorched pine; it was identified in autumn 2016 and we decided to look for the species at the same place in 2017.
On 9.05.2017, we searched for the species on living insolated pines with visible scorching of the outer layer of bark (Fig. 2). The trees did not show any signs of weakness and, had well developed crowns with green needles. On the black, scorched surface of bark, we saw small (1.4-2.0 mm in diameter) rounded exit holes of S. linearis from the previous year (Fig. 3). Whittling the outer bark layer away exposed galleries in the inner layer or under the bark (Fig. 4). Tortuous, variously directed, frequently inter-crossing or overlapping galleries 1.0-2.2 mm in width, filled with fine brown sawdust, did not enter the xylem. Humidity in the feeding places was moderate with a tendency to low. In the galleries, one dead individual, remnants of five others and a living larva of S. linearis were found. Galleries occurred only in the scorched part of the trunk, at a height of 0.3 to 2 m, all around the tree.
Thus, it seems that Stephanopachys linearis prefers thinner trees (it has not been found on the thickest ones), most frequently those of 15-34 cm diameter at breast height (14 trees) and only once it was found feeding on a pine of 38 cm diameter. On the basis of our observations, this species feeds only on living trees, on small subcortical surfaces, where fire had locally damaged the cambium and phloem; it heavily infests scorched pines (up to the height of 135-350 cm), avoiding those blackened only below 130 cm.
S. linearis was accompanied by the deadwatch beetle Ernobius mollis mollis (Linnaeus, 1758) (Ptinidae, Ernobiinae); their galleries were almost identical, but the anobiid prefers those places where the fire had killed the cambium, its larvae feed deeper, frequently damaging the cambium and outer layer of xylem. Its exit holes are somewhat larger and less regular.

Bionomy of Stephanopachys linearis
Stephanopachys linearis is bound biologically to conifers: in the northern zone mainly with Scotch pine (Pinus sylvestris L.), in the alpine areas on various species of pines, larches, much less frequently spruces and firs (Brustel et al. 2013, Ranius et al. 2014, Nardi and Audisio 2016. As a pyrophilous species, it primarily infests living trees damaged by fire, being only rarely found on those showing mechanical wounds like bark flayings or local necrosis. Before intensive forest management, the species probably proliferated on areas of natural fires and, when fire did not occur, on lightningdamaged trees. Like the majority of bostrichids, swarming in Poland occurs in spring, from late April to June. As a representative of the subfamily Dinoderinae, females probably gnaw into the bark to lay eggs. Larval galleries, filled with brown cortical sawdust mixed with faeces, terminate with ovate pupal chambers, from which imagines gnaw their round exit holes through bark. The beetles may develop on the infested tree during at least several years (Wikars 2006). The galleries of S. linearis run at the contact surface between bark and phloem, but larvae feed only in bark, not in living phloem, selecting the layer of dead dry phloem previously killed by fire and then overgrown with new, living tissue. The species prefers the unusual environment produced by, for example, fire. As thick bark protects phloem from being killed by fire, the beetles infest mainly thin-barked trees where the subcortical tissues could have been more easily "boiled". Infestation occurs only after partial regeneration of phloem has occurred, usually 2-5 years after the fire (Wikars 2006).

Status of Stephanopachys linearis in Europe
Stephanopachys linearis is the 20 -th -first of the family Bostrichidae (Gutowski and Przewoźny 2013) -of the beetle species listed in Annexes I and II of the Habitats Directive (Council Directive 92/43/EEC), which have been found in Poland. As a rare species, important for the European Union, it has been placed in Annex II of the Habitats Directive (code: 1926) (Gutowski and Przewoźny 2013; Council Directive 97/62/ EC), as well as in the IUCN Red List and European red list of saproxylic beetles (Nieto and Alexander 2010) -in both cases with the category LC -and Red list of near-extinct and threatened animals in Poland (Pawłowski et al. 2002) as probably extinct (EX?). S. linearis could probably be considered at the Polish regional scale as "Critically Endangered (CR)" following the IUCN criteria B1ab (i,ii,iv) (IUCN 2014).
Apart from S. linearis, considered as an obligatory pyrophile, many other beetle species are associated with burned areas, which -as important causes of natural ecological disturbances -make a significant element of the strategy for preservation of biodiversity (Wikars 1992). In the past, ground fires were frequent in the Białowieża Primaeval Forest (Niklasson et al. 2010, Zin et al. 2015, as in Scandinavia (e.g. on the Swedish island Gotska Sandön - Niklasson 2015). The occurrence of many rare pyrophilous beetles (Gutowski and Jaroszewicz 2001, Gutowski et al. 2012, Niklasson 2015 in these areas is probably a consequence of that "igneous" past.

Concluding remarks
Unfortunately, despite its being listed in the EU Habitats Directive, S. linearis remains poorly known with regards to its biology of development, and many data from various reports placed from time to time in the internet are contradictory, making them difficult to interpret reliably.
The confirmation of the occurrence of S. linearis in Poland will probably stimulate more focused attention to the trees -especially pines -damaged by fire. The authors are convinced that the species can also be found elsewhere in the northeastern parts of Poland, perhaps on larger surfaces with more numerous infested trees (in Białowieża the relevant area was relatively small: 7 ha, 16 pines with galleries of S. linearis). This would allow more detailed studies hopefully resulting in the clarification of its bionomy, at the level of other, well recognised species listed on the EU Habitats Directive.

Different forest cover and its impact on eco-hydrological traits, invertebrate fauna and biodiversity of spring habitats
Martin Reiss 1 , Peter Chifflard 1

Introduction
Mountainous headwater springs are mostly small water bodies where dominant groundwater occurs at the surface to form an intermittent or permanent discharge influenced by subsurface interflow and overland flow (Frisbee et al. 2013). Spring habitats are ecotones in two spatial dimensions. Firstly, springs are vertical aquatic below-above-groundinterfaces between groundwater and surface water. Secondly, springs are above-groundhorizontal terrestrial-aquatic interfaces between different adjacent areas within a continuum of aquatic, semi-aquatic or amphibian to hygrophilous terrestrial spring areas similar to the riparian zone ecotone (Gibert et al. 1990). Therefore, springs are normally nearly constantly cold and wet island habitats without clarifying benchmarks on the delimitation of hot springs (e.g. Pentecost et al. 2003). Glazier (2009) introduces the term "ambient springs" to characterise the water temperature regime with low minimum and maximum amplitudes and an annual mean related to the local yearly mean air temperature (Thienemann 1924). The relative constantly cold water temperature conditions or isothermy (Odum 1971, Teal 1957 lead to specific fauna adaptation of cold-stenothermic specialists and glacial relicts (Nielsen 1950, Ward and Stanford 1982, Fischer 1996. However, fauna composition is much more heterogeneous in springheads and is often explained as a consequence of a high degree of substratum heterogeneity, which can result in high species richness (Bonettini and Cantonati 1996, Lindegaard et al. 1998, Williams and Williams 1999, Staudacher and Füreder 2007, Koperski et al. 2011). The springhead can be built of different substrate types that build heterogeneous mosaic-like structures called patches, hence, springs which can also be termed patchy ecosystems on the micro scale (Pickett and White 1985). The first quantitative approach for analysing the microhabitat-fauna-relationship in springheads and for calculating ratios of invertebrate substrate preferences regarding the ecotone characteristics was completed by Reiss (2011) and considered mountainous headwater springs in Central Germany. The objective of this study is to determine effects from different forest cover by comparing spring habitats in deciduous and coniferous forests on eco-hydrological structures and biodiversity. This research focuses on impacts from forest types as a determinant of the occurrence of corresponding microhabitat types, its substrate type composition and diversity, as well as its specific colonisation by invertebrates. Here, acidity is also considered, because conifer forests could be contributing to acidification of surface waters (Nisbet 1990, Cannell 1999.

Study area
Study sites are located in 6 different forested parts of the Low Mountain Ranges in Central Germany (Fig. 1, Tables 1, 2). They were originally chosen to guarantee a wide range of hydro-morphological structures within diverse substrate types as microhabitats for inver-  tebrates. In total, 86 springheads were analysed, split amongst 61 springs in deciduous forest land cover dominated by Fagus sylvatica (Beech) and 25 springs in coniferous forest land cover dominated by Picea abies (Spruce).

Data analysis and field methods
Data analysis was based on a data set from a study on the microhabitat-fauna-relationship and the importance of invertebrate substrate preferences (Reiss 2011). Here, the evaluation of possible effects from different forest cover on abiotic parameters and biodiversity was not an objective of the previous research. However, data regarding adjacent biotope types were mapped. The determination of new taxa was also added. Adjacent biotope field mapping was done by observing a length of 100 metres from the springhead and by considering four separated quarters, orientated by compass directions (Fig. 2) (Zollhöfer 1997).
Hydro-morphological structure mapping and invertebrate sampling was conducted using a novel integrated technique for multi-habitat sampling (Reiss 2011, Reiss andChifflard 2015). The technique was based on a modified approach and substrate type terminology was used similar to the AQEM/STAR procedure (Hering et al. 2003). Here, the substantial difference is the subdivision of mineral and organic substrate types within the determination of the coverage ratios of microhabitats. Mineral and organic layers were considered individually in a 2-layer approach by taking the area of the whole springhead habitat as a reference surface (5-10 m 2 ). Distribution of sub-sample replicates in each layer should be done according to the estimated share of microhabitats (100 percent coverage means 20 samples in total per layer: 1 sample taken for 5 percent coverage). For example, if the sampling layer of mineral microhabitats consists of 50% psammopelal (mud), 25% microlithal (coarse gravel), 25% macrolithal (coarse blocks and head-sized cobbles) then 10 replicates of benthic macroinvertebrate samples should be taken in the psammopelal, 5 replicates in the microlithal, 5 replicates in the macrolithal. Fauna sampling corresponds with each 5 percent coverage fraction by sampling invertebrates for 2 minutes over a 10 cm x 10 cm reference area for each fraction. Invertebrates were preserved in ethanol alcohol (90%) for transport and subsequent identification in the laboratory.
Data preparation was carried out by filtering and sorting all data related to 100% of the deciduous and coniferous forest cover within the four-segment approach of the adjacent biotope type field mapping procedure. Multiple methods were used to characterise diversity indices: Margalef richness (d) after Margalef (1958), Shannon index (H') after Shannon and Weaver (1949) and Pielou's evenness (J') after Pielou (1966), using the statistical analysis tool PRIMER-E (Clarke and Gorley 2006). Box-Whisker-Plots were generated with BoxPlotR using the online version (shiny.chemgrid.org/boxplotr) (Spitzer et al. 2014).

Hydro-physical-chemical parameters
The results from hydro-physical-chemical in situ measurements in springs (Table 3) indicated that there is not much difference between deciduous forest and coniferous forest land cover for water temperature (Wtemp) and electrical conductivity (EC). In contrast, differences are obvious for pH, oxygen concentration (O 2 Conc.) and oxygen saturation (O 2 Satur.).
The acidification of spring water seems to be one major effect in conifer forest springheads, showing that the pH measurements giving a much lower median value in coniferous forest with pH 5.4, than in deciduous forest with pH 7.7 (Fig. 3).

Microhabitats
Microhabitats, which cover ratios for deciduous and coniferous forest springs (Fig. 4), highlight differences within the two dissimilar forest types. Generally, 13 substrate types are equally present in deciduous forest and coniferous forest springs. However, cover ratios of substrate types in coniferous forest springs are obviously rarer, so that most of the microhabitats, especially organic substrate types are under-represented. Otherwise, certain microhabitats, such as coniferous litter and moss cushions, show slightly higher cover ratios in coniferous forest than in deciduous forest springs. Emergent macrophytes in coniferous forest springs are considerably under-represented in comparison to springheads in deciduous forests.

Fauna and Biodiversity
Overall, 52 taxa and 2617 individuals in springs of deciduous forest land cover and 33 taxa and 326 individuals in springs of coniferous forest land cover were found. Tables 4 and 5 show taxa lists of found invertebrates including individuals per taxon,   occurrence of a taxon (no. of springs), the ratio of individuals and springs (occurrence of a taxon) and the classification of the substrate preference. Typical spring dwelling taxa are in bold letters. In deciduous forest, the ratio of the total number of individuals in a taxon and their level of occurrence (no. of springs) is 43 individuals per spring and, in coniferous forest, 13 individuals per spring. The average number of springs, colonised by individuals of a taxon, is much higher for deciduous forest springs (Ø 10 no. of springs) than it is for coniferous springs (Ø 3 no. of springs). The identification of substrate preferences is also higher for springs in deciduous forest (49 identified with 12 very good substrate type preference identification), than in springs in coniferous forest (40 identified with 7 very good substrate type preference identification). Taxa, with an organic substrate type preference, are much more dominant in both forest types than taxa with a mineral substrate type preference. A clear decline of cold stenothermal species, such as the two groundwater amphipods Niphargus schellenbergi, Niphargus aquilex and the spring-dwelling flatworm Crenobia alpina, is particularly apparent in coniferous forest springs regarding the individuals per spring ratios. The number of spring dwelling taxa -as a group of well-known crenobionts and crenophilous taxa -is almost equal in both forest types. In deciduous forest springs, 11 spring dwelling taxa and, in conifer forest springs, 10 spring dwelling taxa occur. Certainly, the median of the individuals per spring ratio is different; while in deciduous forest springs, it is 4.7 individuals per spring and in conifer forest springs, it is 2.8 individuals per spring. Overall, the amount of spring dwelling taxa is higher in proportion to the total number of taxa in conifer forest springs than in deciduous forest springs. Characterising biodiversity, using biodiversity indices like Margalef richness (d) and Shannon index (H') with deciduous and coniferous forest cover, shows only marginal differences (Table 6). Overall, biodiversity is slightly  higher in deciduous forest springs. The distribution of individuals within occurring invertebrate taxa is more homogenous in coniferous forest springs regarding Pielou's evenness (J').

Discussion
Results show three main forest land cover impacts on eco-hydrological structures and biodiversity in mountainous headwater springs concerning deciduous and coniferous forest: 1) Acidification in conifer forest springs; 2) Higher cover ratios of organic substrate type based microhabitats in deciduous forest springs; 3) Higher biodiversity in species richness and total number of individuals as well as abundances in deciduous forest springs. Acidification of springs in coniferous forests of the German Low Mountain Ranges is a well-known problem, especially in siliceous springs with poor acid buffer capacities and can causes a decrease in the total number of species and individuals (Klös 1984, Lükewille et al. 1984, Meijering 1984, Mauden 1994, Orendt and Reinhart 1997, Audorff and Beierkuhnlein 1999, Hahn 2000. Therefore, it can be expected that the lower values of the total number of species and individuals in conifer forest springs is mainly a consequence of acidification. This interpretation is congruent with findings from a multivariate statistical analysis to characterise the main components of environmental factors influencing the absent and present data on spring fauna using a principal component analysis in the entire data set (Reiss 2011). It is important to note that only in-situ measurements (on-site pH values) were taken into account and no further investigations could be executed to analyse the dynamics or additional reasons for low pH values in order to deduce more accurate and detailed interpretations about acidification.
Higher cover ratios of organic substrate type based microhabitats in deciduous forest springs compared with springs in coniferous forest can be seen as a higher hydromorphological representation of specific microhabitats (e.g. emergent macrophytes or CPOM). However, conifer forest springs have the same quantitative share of occurring substrate types. CPOM (Coarse particular organic material), mainly composed by leaf litter from deciduous trees, is more important as allochthonous organic material as the most important food source for spring fauna communities with a high proportion of primary consumers in particular like herbivorous shredders (e.g. Gammarus fossarum, Gammarus pulex) or grazer and scrapers (e.g. Bythinella dunkeri, Bythinella compressa) (Teal 1957, Mninshall 1967, Whiles and Wallace 1997, Kemp and Boynton 2004. There is no microhabitat preference in coniferous litter, because it is commonly known that the chemical composition of needles tend to retard biological activity in decomposition (Olson 1963). Here, detritus processing takes a much longer time (Sedell et al. 1974, Whiles andWallace 1997) and can be characterised as a low quality detritivore food resource (Taylor et al. 1989, Klemmedson 1992, Friberg and Jacobsen 1994. Furthermore, comparing deciduous and coniferous spring microhabitats, one main conclusion is that substrate types are more important than substrate heterogeneity in contrast to a more general finding by Kubíková et al. (2012). However, in this study, it is not possible to distinguish between acidification and substrate types for describing a prevailing key factor causing taxa and individuals' richness in forest spring habitats.
Higher biodiversity in species richness, total number of individuals as well as abundances in deciduous forest springs in contrast to conifer forest springs, are the most obvious impacts caused by different forest stands. The result of a higher amount of spring dwelling taxa in conifer forest springs than in deciduous forest springs, in proportion to the total number of all found taxa, is consistent with findings by Hahn (2000). Hahn (2000) concluded that spring dwelling species are more tolerant to low pH-values than rhithrobiontic ones due to their adaptation to special crenal living conditions. However, conifer forest springs in our investigation, with their tendency to acidification, are characterised by considerably less individuals per spring ratios and a total number of individuals considering only spring dwelling taxa. Higher cover ratios of organic substrate type based microhabitats in deciduous forest springs compared with springs in coniferous forest seem to be less crucial for biodiversity. Nevertheless, although organic substrate type based microhabitats are of importance, they are more decisive on taxa composition. Some taxa which are under-represented in deciduous forest springs become more quantitatively relevant in conifer forest springs. The abundance of the shredding plecopteran Leuctra spp. is much higher in conifer forest springs, which is in line with the findings in a study about leaf litter decomposition in macroinvertebrate communities in pine and hardwood headwater catchments of the southern Appalachian Mountains in North Carolina, USA (Whiles and Wallace 1997). Leuctra spp. is known as an acid-tolerant aquatic invertebrate (Dangles and Guérold 2000), which is an indicator for acidification as an important factor on taxa composition. Otherwise, acid-sensitive amphipod Gammarus fossarum is similarly abundant, which means no acidification shift in taxa composition is noticeable. It can be inferred that both factors -acidification and microhabitat -govern fauna assemblages or biodiversity and no strong key variable can be characterised. However, this study indicates impacts from different forest cover when comparing spring habitats in deciduous and coniferous forest on eco-hydrological structures and biodiversity.

Conclusion
Different forest land cover causes considerable contrasts in microhabitat structures; obvious organic substrate type composition and cover ratios; as well as differences in species richness and invertebrate abundance of spring habitats in deciduous and coniferous forest. This means, land cover as an ecological mesoscale, properly determined by different forest types, has an impact on eco-hydrological structures and biodiversity on the micro scale. It implies an essential consideration for adjacent biotope type mapping and is an important integrative parameter for spring habitat assessment approaches. Furthermore, the recognition of substrate preferences of invertebrates, within an ecotone based assessment approach, characterises microhabitats explicitly for all parts of a springhead, regarding aquatic and terrestrial spring habitat zones. Here, the importance of forest land cover and the substrate type diversity relationship is taken into account within an ecological spring habitat assessment methodology and characterises its consequences on invertebrate biodiversity. Therefore, negative effects from forest management practices (e.g. forest conversion) within a nature conservation perspective can be included in decision-making and action plans to realise national or regional strategies on biodiversity.