Microhabitat selection and communal nesting in the insular Psychedelic Rock Gecko, Cnemaspis psychedelica, in Southern Vietnam with updated information on trade

The Psychedelic Rock Gecko, Cnemaspis psychedelica, was described in 2010 and certainly belongs to the most spectacular gecko discoveries worldwide. The species is endemic to two small offshore islands in Rach Gia Bay. Its striking colour pattern makes the species highly attractive for the international pet market. The existent Cnemaspis population is negatively affected by habitat degradation and predation by introduced macaques. We herein provide the first characterisation of microhabitat selection of this species, including seasonal variation on Hon Khoai and Hon Tuong islands, Ca Mau Province, Vietnam. We found that characteristics of the selected microhabitat, such as substrate type, temperature and canopy cover slightly differed between the wet and dry seasons. We also demonstrated age-related differences in the selection of perch heights. Communal nesting was, for the first time, reported for C. psychedelica, as well as natural predation by a snake species (Lycodon capucinus). In addition, we documented ongoing habitat destrucNature Conservation 31: 1–16 (2018) doi: 10.3897/natureconservation.31.28145 http://natureconservation.pensoft.net Copyright Hai Ngoc Ngo 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
The old world gecko genus Cnemaspis Strauch is considered as one of the most speciesrich paleotropical gekkonid genera, encompassing more than 100 species with a widespread distribution throughout tropical Africa, South Asia and Southeast Asia (Smith 1935, Wickramasinghe and Munindradasa 2007, Ganesh et al. 2011, Uetz and Hošek 2017, Vidanapathirana et al. 2014, Srinivasulu et al. 2015. To date, six species of Cnemaspis have been recorded from Vietnam, namely C. aurantiacopes, C. boulengeri, C. caudanivea, C. nuicamensis, C. tucdupensis and C. psychedelica. Two of these from Vietnam are endemic to islands Ngo 2007, Grismer et al. 2010). One of them is the insular Psychedelic Rock Gecko, Cnemaspis psychedelica, which was only described from isolated Hon Khoai Island, Ca Mau Province, Southern Vietnam in 2010 (Grismer et al. 2010). More recently C. psychedelica was also recorded on a further offshore island, namely Hon Tuong isle, which only covers a small area of ca. 300 m 2 (Ngo et al. 2016). Still, the species is assumed to be endemic to islands of Rach Gia Bay, southern Vietnam.
The known wild population of the species was estimated to comprise approximately 500 mature individuals (Ngo et al. 2016). Although the wild population of C. psychedelica appears to be relatively stable and actively reproducing, it is currently suffering from increasing habitat degradation on the small Hon Khoai Island (Ngo et al. 2016). Particularly, those granite formations, which represent essential microhabitat sites for C. psychedelica are frequently blasted by dynamite in order to flatten several areas on the island for construction of roads or artificial ponds (Ngo et al. 2016). Additionally, poaching for the international pet trade has already been recorded shortly after description of the species. C. psychedelica has been observed in European reptile fairs and on the online markets for prices of up to 3,000 EUR/pair (Auliya et al. 2016, Ngo et al. 2016). In addition, the introduction of invasive Long-tailed Macaques (Macaca fascicularis) on Hon Khoai Island poses another potential threat to C. psychedelica, as macaques were observed to feed on geckos and theirs eggs (Grismer et al. 2010, Ngo et al. 2016). The extremely restricted habitats, together with a low reproduction rate, make the species especially vulnerable to external stressors and limit its capacity to recover from threats such as harvesting. As a result, the protection needs of C. psychedelica have received attention from all around the world. Based on a first population and risk assessment, the species has recently been classified as "Endangered" on the IUCN Red List of Threatened species  and was subsequently listed in Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES 2016). Furthermore, immediate conservation measures have been initiated including the establishment of an in-country ex situ assurance population for C. psychedelica ). However, detailed information about microhabitat requirements of the species, which are essential for both adequate habitat protection measures and to successfully implement conservation breeding programmes, are still lacking.
The study provides detailed data on microhabitat selection of C. psychedelica within its natural habitat. We examined differences in microhabitat selections of lizards between the wet and dry seasons, as well as variations in habitat use amongst age classes. Furthermore, we provide new information about egg deposition behaviour, natural predation and ongoing habitat destruction on Hon Khoai Island, as well as new evidence for local trade in Cnemaspis psychedelica, highlighting the need for immediate measures to protect the remaining populations.

Field surveys
Study sites were selected based on known occurrences of C. psychedelica on Hon Khoai and on Hon Tuong islands, Ca Mau Province, southern Vietnam, according to previous surveys (Grismer et al. 2010). The region is characterised by a monsoon sub-equatorial climate with constant annual temperatures, but with distinct climatic difference between a dry and a wet period with heavy rains (Fig. 1) (General Statistics Office of Vietnam 2016).
Field surveys were carried out during the wet season in November 2015 as well as during the beginning of the dry season in January 2016. Another short survey on Hon Khoai Island took place in January 2017. We mainly conducted night excursions after sunset between 19:00 h and 24:00 h. In order to determine the sex of the geckos and for taking measurements, individuals were captured by hand and subsequently released at the same spot in the morning or afternoon of the next day between 10:00 h and 17:00 h. While releasing geckos, further sighted individuals were also recorded, measured and released. Each gecko was measured with a digital slide-caliper to the nearest 0.1 mm. Abbreviations are as follows: snout-vent length (SVL), measured from tip of snout to anterior margin of cloaca; tail length (TL), measured from posterior of cloaca to tip of tail in geckos (Grismer et al. 2010).
In order to investigate the biodiversity of Hon Khoai Island, we carried out respective field surveys in January 2016. During these surveys, a few amphibians and reptiles were also collected and subsequently deposited in the collection of the Institute of Ecology and Biological Resource (IEBR), Hanoi, Vietnam. Stomach contents of collected snakes were analysed after dissection in order to identify whether they might represent natural predators of C. psychedelica.

Microhabitat assessment
Microhabitat parameters were recorded for each sighted individual, including substrate type (classified as cliff, rock, branch, leaves, forest floor), perch height [in m] (vertical distance between captured animal and ground), percentage of vegetation or rock coverage above animal, position (resting outside, under a rock or within a crevice), substrate surface condition (dry or wet), activity (resting, feeding, foraging), air temperature [°C], substrate temperature [°C], animal's body surface temperature [°C] and relative air humidity [%]. Air temperature and relative air humidity were measured with a digital thermo-hygrometer at a vertical height from the ground up to 2m height at each microsite (TFA Dostmann/Wertheim Kat.Nr.30.5015), while substrate and body surface temperatures of animals were determined with an infrared thermometer (Measupro IRT20).
To identify intraspecific differences in microhabitat selection of C. psychedelica, individuals were classified into different age classes according to Ngo et al. (2016) based on snout-vent lengths (SVL): SVL < 58 mm = juvenile and SVL ≥ 58 mm = adult.
Chi-square tests and t-tests with P ≤ 0.05 were performed to determine significant differences in selected microhabitat parameters with categorical and continuous variables, respectively, between wet and dry season and amongst age classes. Statistical analyses were applied with the programme SPSS Version 16.0 (SPSS Inc., Chicago).

Threat assessment
To obtain an overview on the availability of and evidence for trade in Cnemaspis psychedelica in Vietnam, we visited several local pet markets (one in Ca Mau Province where the species occurs, one in Dong Nai Province, one in Ho Chi Minh City and one in Ha Noi City) and further investigated different internet platforms, reptile Facebook pages and Forums. We also interviewed three local dealers offering the lizard online and two local keepers in March 2018, in order to obtain information on origins, commercial prices and networks of illegal trade in this species. Additionally, we conducted field surveys at previous survey sites located on Hon Khoai Island, Ca Mau Province in January 2017 to evaluate anthropogenic threats. Nearly 10 road workers on Hon Khoai Island were interviewed to determine the local use of the species. Names of interviewees were kept anonymous to ensure data privacy rights and internet links were not disclosed to prevent misuse.

Microhabitat use
The Psychedelic Rock Gecko was found active during both the wet and dry seasons. Regarding daily activity pattern, C. psychedelica was observed active during any sighting between 10:00 h to 24:00 h. The vast majority of lizards were found in the shade, even if a patch of direct sunlight was in close proximity. During the wet season, C. psychedelica was mainly found on loose granite rocks, followed by cliff and leaves on trees (61%, 38%, 1%, n = 413, respectively, see Fig. 2A), but never on branches or on the forest floor. During the dry season, the main selected substrate types were cliffs, followed by loose granite rocks, leaves on trees and branches (65%, 19%, 14%, 2%, n = 156, respectively, Chi 2 = 36.4, df = 3; P < 0.001, see Fig. 3A-D). We mainly encountered lizards resting under rock formations during the wet season, while animals frequently resided outside of caves or rock shelters during the dry season (outside of caves: 38.5% during the wet season vs. 56.7% during the dry season, respectively, see Fig. 2B, Chi 2 = 19.1, df = 2, P < 0.001). We further observed C. psychedelica at mean heights of 0.71 ± 0.6 m (0.01-3 m, n = 569) (see Table 1). There was no difference in the height selections of lizards between wet and dry season (t-Test, t = 0.86, df = 567, P > 0.05) (see Fig. 2E). While animals tended to reside under a mean canopy coverage of 79.6 ± 30.5% (n = 156) during the wet season, a higher percentage of canopy cover of 91.2 ± 20.1% (n = 413) above the animals was recorded during the dry season (t-Test, t = 4.4, df = 208, P < 0.001) (see Fig. 2D, Table 1). We observed that the majority (about 75.9%, n = 432) of lizards were resting during surveys, whereas only about one quarter of individuals (n = 137) were found actively foraging.

Oviposition
Egg depositions of Cnemaspis psychedelica happened on flat surfaces of cliff walls, without direct sunlight. Egg deposition sites were located under high vegetation coverage during both seasons (see Fig. 4A). We did not find any eggs attached to branches or leaves. Newly laid eggs were bright and clear white in colouration, then changed to pinkish (early development of embryos) and then to a slightly grey colour (with developed embryos). Eggs were almost round in shape, except for the flattened side attached to the substrate. Clutches of C. psychedelica consistently comprised two eggs. Clutches of different individuals of C. psychedelica were commonly deposited in close proximity to each other, forming clusters. Furthermore, new clutches were frequently found to be deposited on top of remains of eggshells from previous clutches. We measured the min- imum distance amongst different clusters as about 25 cm. Furthermore, the maximum observed number of unhatched eggs within a cluster was 10 eggs/5 pairs (see Fig. 4A, B). The clusters were generally deposited at heights ranging from 0.3 to 3.5 m above the ground. We frequently found several adult individuals in close proximity to the aggregation of clusters. Solitary clusters comprising only two eggs were usually observed higher than 2.5 m above the ground and were scarcely attended by adult individuals.

Threats
During the survey in January 2017, the optimal habitats of C. psychedelica on Hon Khoai Island were found to have experienced further degradation due to expanding road construction and the building of further artificial ponds (see Fig. 6A, B). In particular, destruction with dynamite has entirely removed the natural habitat of C. psychedelica at one survey site and further extended to other survey areas.
Trade in living C. psychedelica has been recorded by our team in local pet shops from northern Vietnam (Ha Noi City) and southern Vietnam (namely, Ca Mau Province where the species is distributed, further in Dong Nai Province and Ho Chi Minh City), as well as online in Vietnamese Reptile Forums, Zalo Online and Facebook. During our visit to two pet shops, one in Ca Mau Province and one in Dong Nai Province in March 2018, we observed 20 and only one remaining C. psychedelica, respectively, which were kept in small boxes (Fig. 5B). A large number of animals had already been previously exported. Local dealers confirmed that C. psychedelica have been caught from the wild for sale since 2015. According to interviews with local dealers, 15 to 70 specimens from Vietnam have been frequently offered for sale and export to Thailand and Indonesia (for USD 100-300 per individual) and afterwards smuggled to Russia and Europe for USD 700 per individual, while prices achieved in the local trade were comparably lower (USD 20-40 per animal). According to one private keeper in Ho Chi Minh City, C. psychedelica have been offered for USD 450 per pair from another local pet shop in Ho Chi Minh City since 2016. In March 2018, he had 15 breeding specimens for sale for USD 200 per individual and several eggshells were observed in his terrarium. Regarding international trade, at least 23 specimens of C. psychedelica were found being offered for sale online (such as on Facebook) in Europe between September 2017 and March 2018 (Altherr et al. pers. comm.).
During the survey in January 2017, we found some reptile species such as Cnemaspis psychedelica, Cyrtodactylus leegrismeri, Draco maculatus, Hemidactylus frenatus, Gehyra mutilata, Gekko gecko and Ahaetulla prasina soaked together in rice wine and used for traditional medicine by road workers on Hon Khoai Island (Figs. 5C and D). According to interviews with road workers, C. psychedelica individuals have been kept as pets in their terrariums at home since 2015.

Microhabitat
We found that selected microhabitat characteristics such as substrate type, temperature and canopy cover slightly differed between the wet and dry seasons. Our data showed that the species is more frequently found under granite rock formations during the wet season, probably to avoid heavy rains (10 of 11 survey days). During the dry season, geckos were observed in remarkable numbers outside of caves or rock formations. The species was found strongly associated with granitic rock formations and cliffs during all life stages. While C. psychedelica has never been found on the forest floor and only scarcely within the vegetation, clutches were also unexceptionally deposited on rock formations at heights of at least 0.3 m above the ground. We assume that vertical cliffs might provide shelter from ground-dwelling predators.
Furthermore, animals were found at spots which were more densely covered with canopy during the dry season than during the wet season, which could be their behaviour to avoid direct sun exposure. Accordingly, fewer animals were found along transects with only slight vegetation coverage during the wet than during the dry season. These findings suggest the dependence of C. psychedelica on the availability of shaded habitats. Egg depositions were also exclusively found on cliffs, which are facing away from the sun.
Our study revealed that the body surface temperature of C. psychedelica was positively correlated with the substrate temperature (r s = 0.51, P < 0.001, n = 567). Thus, as in other ectotherms, basic physiological functions of C. psychedelica, such as locomotion, growth and reproduction are determined by the environmental temperature. Since tropical lizards are considered to have narrow temperature optima and only few options for behavioural and physiological compensation, they are assumed to be, in particular, vulnerable to extinction by climate warming (Deutsch et al. 2008;Doody and Moore 2010;Huey et al. 2009;Vié et al. 2009).

Communal nesting and spatial distribution
We frequently observed aggregated oviposition sites and the placement of fresh eggs on top of or close to previous oviposition sites. Communal nesting is defined as "nonincidental deposition of eggs at a shared nest cavity by two or more co-specifics" (Espinoza and Lobo 1996). Such communal oviposition has been described in numerous lizard species including other Cnemaspis species (e.g. Kalaimani 2015;Lima et al. 2011;Magnusson and Lima 1984;Vitt 1986;Soares de Oliveira et al. 2015;Somaweera 2009;Srinivasulu and Srinivasulu 2013a, b;Werner 2002) and is assumed to offer potential benefits such as protection, predator-satiation and metabolic heating (Gurgel de Sousa and Freire 2010, Sönmez 2018). Mateo and Cuadrado (2012) experimentally proved that hatching success was significantly higher in communal than in solitary clutches of the gecko Ptyodactylus oudrii and significantly decreased when adult lizards were excluded from the oviposition site. Accordingly, we also observed high densities of C. psychedelica specimens of different ages and sexes, accumulated in close vicinity to communal oviposition sites following the principle of clumped distribution. The maximum observed density was 23 individuals, including 18 adults (9 males), per 4 m 2 . The clumped presence of individuals, independent of age and sex in close proximity to clutches, has also been observed by Mateo and Cuadrado (2012), assuming a high population density and restricted availability of suitable oviposition sites to trigger communal nesting.
With a closer look on the spatial distribution of C. psychedelica, we found age related differences in perch heights, namely juveniles occurring at significantly lower heights than adults. Similar habitat divergences between juveniles and adult individuals have also been reported for Cat Ba Tiger Geckos and Crocodile Lizards in Vietnam (Ngo et al. in press, van Schingen et al. 2015) and gekkonids from New Caledonia (Snyder et al. 2010). Van Schingen et al. (2015) suggested that climbing might be a trade-off in lizards between saving energy costs and the reduction of interactions with competitors and predators. However, this hypothesis needs to be tested in future studies. Furthermore, adult specimens were usually found at similar heights to the closest communal clusters, which might indicate some kind of parental care behaviour.

Predation
Reptiles, including lizards are commonly prey of other vertebrates, such as mammals and reptiles (Whiles and Grubaugh 1993). The first case of a snake as a predator of Cnemaspis psychedelica, as herein reported, is thus nothing extraordinary (e.g. Herman 2017, Nogueira et al. 2013, Vrcibradic and Eisfeld 2016, Whitaker and Captain 2004. However, Long-tailed Macaques (Macaca fascicularis) as predators, as reported by Grismer et al. (2010), is a worrying fact, as this species was introduced and thus does not belong to the normal ecosystem of Hon Khoai Island. The artificially enhanced predation probably impacts the population dynamics of C. psychedelica.

Anthropogenic threats and implications for conservation
Lizards with limited geographic distributions, such as Cnemaspis psychedelica, may be extremely vulnerable to local habitat loss or alteration (Sodhi and Ehrlich 2010). The limited availability of habitats for the species, namely only two offshore islands in the Rach Gia Bay, Vietnam, makes the species particularly vulnerable to habitat destruction. The optimal microhabitats of C. psychedelica on Hon Khoai Island were found to have experienced ongoing degradation due to expanding infrastructure. If the present trend continues, not only C. psychedelica, but also other species and yet undescribed or just recently discovered diversity on the island (Nguyen et al. 2018), will become potentially threatened by extinction. Furthermore, these constructions will probably further facilitate hunters to illegally approach Hon Khoai Island. C. psychedelica has already been observed on international markets, but detailed information about illegal traffic networks has been lacking (Altherr 2014;Auliya et al. 2016, Ngo et al. 2016. New evidence from interviews with local reptile dealers indicates that they paid local fishermen and visited Hon Khoai Island as tourists to collect live specimens of C. psychedelica without any permits. Live specimens of C. psychedelica have been frequently offered for sale in southern and northern Vietnam and are mainly smuggled to Thailand as a middle country and then exported to be sold in Russia and Europe for relatively high prices. However, the number of smuggled animals has allegedly declined to fewer than a total of 20 individuals per deal which might be an effect of the improved control of trade in the species due to the recent inclusion in Appendix I of the Convention on International Trade in Endangered Species of Wild Fauna and Flora (CITES) as well as promulgating strict decrees that have restricted harvest and reduced the stress from over-exploiting of a number of lizard species (Sodhi and Ehrlich 2010). However, traders still continue to harvest wild individuals and the species is still found to be offered in local and European markets probably due to high demands and commercial value. Our findings indicate that there are currently no legal commercial harvests and trade in the species in Vietnam, implying that wild specimens entering the international trade have been illegally harvested. We highly recommend that the Vietnamese authorities more strictly control the illegal collection of C. psychedelica in the natural habitat. Furthermore, wildlife law enforcement along key border regions should be highly improved to alleviate the illegal exploitation and trade, not only in wild C. psychedelica, but also in other taxa. Nguyen TQ, Ngo HN, Ziegler T, van Schingen M (2016)  Tropical forest succession may be related to the structural and floristic changes, such as changes in plant species richness, stem density and plant species composition (Gibson et al. 2014), which should lead to changes in food resource availability for primates (Bryson-Morrison et al. 2016;Pinotti et al. 2015). The intermediate successional stage is markedly more diverse and abundant than the early stage, but it is not significantly different from the later stage (Kalacska et al. 2004). Tree species composition tends to increase (Kalacska et al. 2004) or changes slightly (Dent et al. 2013) in similarity to old growth forests with forest age. Compared to primary forest, secondary forest is also diverse in terms of primate food species, but large fruit tree are rare (Bryson-Morrison et al. 2016). Understanding which forest types related to different successional stages providing different food resources is vital for the long-term survival of primates.
The Hainan gibbon, Nomascus hainanus (Thomas), is the world's rarest ape and one of world's most endangered mammal species (Bryant et al. 2015;Geissmann and Bleisch 2008;Stone 2011;Zhou et al. 2005). The last surviving population (approximately 26 individuals) is confined to the areas of remaining habitat within Bawangling National Nature Reserve (BNNR), on western Hainan Island, China (Deng et al. 2017;Stone 2011). The species is regarded as "Critically Endangered" in the IUCN Red List of Endangered Species due to an estimated population decline of over 80% within the past several decades, primarily resulting from hunting and habitat loss (Geissmann and Bleisch 2008). Since 1989, the Hainan gibbon has been a First Class Nationally Protected Species under the Chinese Wildlife Protection Law (Zhang et al. 2010). The population was estimated at over 2000 individuals in the 1950s and lived throughout the tropical forest of Hainan Island (Zhou et al. 2005). Deforestation and other anthropogenic activity have led to the species today being restricted to BNNR. Zhou et al. (2005) indicate that the species is unlikely to survive unless efficient conservation policies and strategies are put in place immediately. The apparent threat of extinction that has been known since 2003 led to a series of conservation projects such as reforestation of degraded habitat including planting gibbon food trees that were initiated by local forestry departments, Kadoorie Farm & Botanic Garden (KFBG) and other nonprofit organisations (Stone 2011). Although an initial success of the present restoration effort can be seen, basic knowledge regarding the distribution of food trees that are important for Hainan gibbon long-term survival remain unclear (Fellowes et al. 2008).
Accordingly, the overall goal of this study was to assess diversity patterns of food trees and variation in community composition across different successional stages of natural forests, as well as plantations, to improve our understanding of food distribution and succession in BNNR. Our study provided basic information on distribution patterns of food trees species for the Hainan gibbon to ensure the long term conservation and survival of this unique species.

Study area
BNNR (18°57.15'-19°11.21'N, 109°03.32'-109°17.51'E) is located in west Hainan Island, south China ( Fig. 1), at the northern edge of the Asian tropical forest zone (Ding et al. 2012). It covers an area of 482 km 2 , with its highest peak being 1654 m above sea level ( Fig. 1). BNNR was established as a National Nature Reserve in 1988 with the main objective of conserving the Hainan gibbons and their habitat. The region is characterised by a tropical monsoon climate with a distinct wet season from May to October and a dry season from November to April. The average annual precipitation is 1750 mm and average annual temperature is 23.6°C. Soils at higher elevations are mixtures of red loam and yellow loam, while latosol, developed from granite, dominates soils at lower elevations (Zang et al. 2010).
The landscape in BNNR is a mosaic of natural vegetation and exotic species plantations in which natural tropical forests dominate. Due to severe and repeated anthropogenic disturbances over the past 40 years, forest landscapes in BNNR have become increasingly fragmented and old growth forests have been progressively replaced by secondary forests in various development stages, shrub/grass land and rubber or pine plantations (Zang et al. 2010). Logging operations may modify the diet of Hainan gibbons directly by cutting down important food trees and indirectly by transforming original forests into secondary vegetation or plantations (Morgan and Sanz 2007). In summary, it can be said that the past anthropogenic activity in the region has gradually destroyed habitats and food resources for Hainan gibbons (Zhou et al. 2005).

Data collection
Gibbon food trees were defined as tree species whose edible parts (leaves, flowers, fruits or seeds) are known to be part of the gibbon's diet in the area. We used food plant lists from long-term studies of the Hainan gibbons in the tropical forest region of BNNR (Chen et al. 2009;Lin et al. 2006;Wu 2007;Zhou 2008) to identify food plants amongst the surveyed species. We investigated a total of 149 sample plots (20 × 20 m) based on 1 × 1 km or 1 × 2 km grid across the BNNR forest region (Fig. 1). At each grid point, an inventory plot was established and positioned by a hand-held global positioning device (GPS) for the vegetation survey. Sample plots cover a total area of 5.96 ha in four forest types. The forest types were classified based on time since the harvest or the successional stage of the forest was within the period of the inventory. For the classification of successional stages, information available from timber harvest archives of the BNNR was used. In addition, we employed experienced loggers to identify the successional stages in each sample plot. Of the total number of plots, 33 fell within the young natural secondary forests (YSF; < 25 years since disturbance), 26 in middleaged natural secondary forests (MSF, 25-60 years since disturbance), 76 in old growth forests (OGF, > 60 years since disturbance) and 14 in Pinus merkusii plantation forests (PF, 20-35 years).
In each sample plot, all free standing woody stems with diameter at breast height (DBH) ≥ 1 cm, were counted, measured and identified to the species level. Species nomenclature follows Flora of China (http://www.efloras.org). In total, we surveyed 60,930 woody stems belonging to 587 species, 275 genera and 82 families.
Adult trees generally produce a greater quantity of food resources than younger trees (Chapman et al. 1992) and only those trees with a DBH ≥ 10 cm are large enough for Hainan gibbon to forage in, so we used tree DBH as an indicator to determine suitable food trees. Tree DBH is commonly used to describe primate resources (Worman and Chapman 2006). In order to account for differences in species richness and individual density of food tree size classes amongst the four forest types, all woody stems were classified into three size classes: saplings (1 cm ≤ DBH < 5 cm), young trees (5 cm ≤ DBH < 10 cm) and adult trees (DBH ≥ 10 cm).
We treated each plot (20 × 20 m) as a unit for all subsequent analyses, except for the accumulation curves. For each plot, we quantified species richness as well as stem density of all food tree species.

Data analysis
To compare food species diversity amongst the four forest types, we computed samplebased rarefaction curves constructed using the analytical formula implemented in Esti-mateS 9.0 (Colwell 2013). Rarefaction curves can be used to standardise samples that differ in terms of individual size or plot size. This procedure was repeated for each of the 200 bootstrap iterations. Patterns of species dominance were compared between different forest types using species rank-abundance plots. The significance of differences amongst forest types with respect to the food species richness for three stem size classes was assessed by one-way ANOVA, followed by a Tukey HSD multiple comparisons test when significant differences (p < 0.05) were found. Prior to the statistical analysis, all data were log-transformed.
We assessed the completeness of each forest type by calculating the number of observed species as a percentage of the total richness and coverage of each forest type by calculating the number of species recorded as a percentage of the average estimated richness, which was estimated based on the average of two abundance-based nonparametric estimators: Chao 2 and jack 2 (Colwell 2013). Both estimators were used to reduce the bias in underestimating species richness and have been shown to be accurate to true species richness relative to other non-parametric estimators, especially when based on small sample sizes (Hortal et al. 2006).
We used ANOSIM (analysis of similarities) to test whether species composition differed amongst four forest types for each DBH size class (saplings, young trees and adult trees) using the ANOSIM function in the Vegan package. ANOSIM was implemented with a maximum of 999 permutations. Pairwise similarities showing community overlay between forest types for each DBH size class were also calculated based on Sørenson's quantitative index (Chao et al. 2005) using the SIM function in the Vegetarian package. All the statistical analyses were performed with R 3.3.2 (R Core Team 2016) unless otherwise indicated.

Results
A total of 85 species belonging to 55 genera in 35 families were identified amongst 16,882 collected food plants (Suppl. material 1: Table S1). The families with the highest number of species were Lauraceae (13%), Moraceae (12%) and Myrtaceae (8%). The most abundant families were Rubiaceae (38%), Lauraceae (12%) and Myrtaceae (6%). Distribution of food tree species amongst forest types was different, ranging from 21 species present in all forest types to 11 species exclusive to one single forest type each (Table 1).

Changes in species diversity patterns
Food trees were at least 70.5% of the estimated total number of species in each forest type (Table 1). This result displays the validity of our comparisons of species richness amongst four forest types. At least 83.5% of all trees were recorded in natural forests with various successional stages, while only 25.9% of all trees were found in plantations. Habitat-exclusive species had the highest richness in OGF, followed by MSF, while we did not find habitat-exclusive species in PF. At plot scale, stem density and species diversity were higher in the natural forests (YSF, MSF and OGF) and the Table 1. Diversity of total sampled food trees suitable for Hainan gibbon in young secondary forests (YSF), intermediate secondary forests (MSF), old forests (OGF) and plantation forests (PF) of BNNR, Hainan Island, China. N = total number of sampled plots, n = number of individuals sampled, S tot = total number of species recorded, Coverage = number of species recorded as a percentage of the average estimated richness, S plot = mean richness per plot ± standard error (ANOVA F 3,145 = 15.3, p < 0.001), D plot = mean density per plot ± standard error (ANOVA F 3,145 = 14.5, p < 0.001), Ŝ MaoTao = species richness at a standardised sampled size, Completeness = number of species recorded as a percentage of the landscape total, Ex tot = total number of habitat-exclusive species. Values designated by the different letters within each variable are significant at p < 0.01.  highest stem density and species richness were found in MSF (p < 0.01) ( Table 1). The analysis between forest types using a standardised sample size (Ŝ MaoTao for N = 14) showed species richness across four forest types were similar to the plot-scale diversity patterns (Table 1).

Forest type
Rarefaction curves for all four forest types were nearly asymptotic, with species richness being highest in MSF and lowest in PF (Fig. 2a). Plantations were dominated by a few very abundant species (Fig. 2b). The four dominant species Psychotria rubra, Syzygium chunianum, Adinandra hainanensis and Machilus chinensis accounted for 75.8% of the sampled individual trees. Natural forests exhibited lower species dominance than plantations and the differences of species dominance amongst various recovery stages were not significant (Fig. 2b).
There were significant differences amongst forest types regarding species richness and abundance of various DBH size classes (p < 0.001 in all cases). All size classes of food trees in plantations displayed significantly lower species richness than each of the three natural secondary forest types (Fig. 3a). MSF showed considerably more food trees richness (Adult: 9.62±0.83; Saplings: 17±1.33; Young: 7.76±0.67) across all size classes per plot compared to OGF (Adult: 6.66±0.57; Saplings: 12.11±0.75; Young: 5.14±0.40). The abundance of adult food trees was higher in MSF (19.08±1.71) and OGF (14.00±1.18) than in YSF (7.57±2.25) and PF (4.25±2.93) (Fig. 3b). No significant differences were found regarding the abundance of saplings and young trees amongst the different forest successional stages (p > 0.05). However, saplings and young trees of food species were more prevalent in natural forests (YSF, MSF and OGF) than in PF (Fig. 3b).

Variation in species composition
Differences in species composition amongst forest types were significant for adult trees, young trees and saplings (ANOSIM test, R > 0.3, p < 0.05 in all cases). Furthermore, pairwise comparisons revealed that PF had distinct species assemblages for all size classes and OGF exhibited strikingly different adult trees (ANOSIM test, R = 0.27, p < 0.05) or saplings species (ANOSIM test, R = 0.21, p < 0.05) composition compared to YSF. There was no significant difference in young trees species composition between different successional stages of natural forests. Floristic similarity for young trees species was, however, slightly higher between MSF and YSF than between OGF and MSF or between OGF and YSF (Fig. 4).

Discussion
BNNR currently remains the last refuge for Hainan gibbons in the world (Stone 2011). The food tree diversity and composition assessment in the BNNR, presented in this study, highlights distribution patterns amongst different forest types and displays important information for the successful protection and conservation of Hainan gibbons in the future. Food abundance and forest type likely have a stronger influence on primate habitat selection than any other factors (Furuichi and Hashimoto 2004).
Compared to natural forests in different successional stages (YSF, MSF and OGF), plantations (PF) have unfavourable attributes, such as a lack of soil humidity, higher solar radiation, lower litter decomposition rates and a greater density of ruderal species (Lemenih et al. 2004;Lohbeck et al. 2013;Nagaike 2002) that result in lower richness and abundance of food trees, most of which are shade-tolerant species (Table 1, Figs 2a,3). PF also had a distinct species composition amongst all size classes (adult trees, young trees and saplings) due to the differences in forest structure and the aforementioned characteristics. Different plantation management practices can also affect food tree diversity. The thinning of objective trees in an appropriate period for timber production can efficiently minimise the decrease in the abundance and richness of the food resource, especially shade-intolerant trees (Sakamaki et al. 2011). However, clear-cutting, most common in BNNR, is likely to reduce the number of food trees and proportion of larger food trees during timber harvesting (Table 1, Figs 2a, 3). In BNNR, planted forests are generally found in the areas below 760 m where anthropogenic disturbances occur more frequently than at higher elevations (Zhang et al. 2010;Zhou et al. 2005). The lack of food resources, lower numbers of large food trees and simple forest structures in plantations have severely limited Hainan gibbon's ability for locomotion, feeding and resting. To date, the gibbons have never been seen in plantations (Zhang et al. 2010). Therefore, converting plantations back into forest, useable by gibbons, is of highest importance.
The development of species diversity along successional gradients may follow this hypothesis: that species diversity increase with vegetation recovery time, but peaks at mid-succession (Auclair and Goff 1971). As expected, a similar pattern of total food species diversity distribution across successional stages was identified in our study as well, regardless of whether or not stem density per plot and the number of species (S plot and Ŝ MaoTao for N = 14) were considered (Table 1, Fig. 2a). Similar results for each DBH size class distribution were also found (Fig. 3). Food species diversity patterns could be explained by trade-offs in dispersal and competitive ability (Howard and Lee 2003). The trade-offs limit species presence at the two successional extremes (YSF and OGF). However, in the mid-successional stage, most tree species have already dispersed and established, but competitive exclusion has not yet eliminated the shade-intolerant food species, resulting in the highest food species diversity. Similarly, the decline in food species diversity in OGF may have resulted from reduced light availability (Lohbeck et al. 2013). As trees become larger, canopy cover becomes denser and low light levels may inhibit food species recruitment from seed or seedling banks, thus limiting food species richness and stems with small DBH class (Fig.  3). Our results could also be explained by the intermediate disturbance hypothesis (IDH), which indicates that species diversity is maximised where intermediate scales of disturbance are experienced (Connell 1978 ) and have therefore been specifically targeted during past harvest operations. As a result, the large food trees preferred by gibbons were lost in young secondary forests and are restricted to primary or old-growth forest adjacent to ravine areas (Liu et al. 1989;Wu 2007). After 25-60 years of succession, the adult food tree diversity increased with recovery age and were equivalent to, or slightly higher than, OGF levels ( Fig. 3), suggesting that prohibition of mechanical logging since 1994 apparently stimulated the growth of large trees that provide food for gibbons.
Our study revealed that differences between OGF and YSF were larger than those between MSF and YSF or between OGF and MSF regarding food species composition in each DBH size class (Fig. 4), indicating that food species composition in secondary forests became more similar to that in old growth with increasing age. Our results showed many similarities to a study by Howard and Lee (2003) on Barro Colorado Island in Panama, where species composition and availability of food resource in secondary forests converged towards those characteristics found in old-growth forests with stand age. This hints at the fact that it is possible to successfully recreate lost habitats. In the tropical forest of north-eastern Costa Rica, Letcher and Chazdon (2009) also found that secondary forests (> 30 years) have similar species composition to the old growth forest. The changes of food species composition across successional gradients may be attributed to the land-use history (Letcher and Chazdon 2009) and conversion of plant strategies from acquisition in the early successional time to conservation in the late successional stage (Lohbeck et al. 2013). Other factors, such as the extent and condition of surrounding forests, can also impact the distribution of species composition (Chazdon 2003). The driving forces influencing food species composition variation, however, need to be further tested based on abiotic factors and measure-based functional traits data in future studies, to ensure a fast and successful conversion of depleted forests back into suitable habitat.
In this study, eleven food tree species were detected exclusively in one single forest type (Table 1). They are very infrequent and an important source of fruits fed on by Hainan gibbons (Chen et al. 2009). For example, Ficus championii was identified as one of the preferred food resources for gibbons (Chen et al. 2009;Wu 2007). Pouteria annanmensis was consumed in the driest months (February-April). It was reported that there were only seven food species consumed by Hainan gibbons during the dry season (Chen et al. 2009;Liu et al. 1989;Wu 2007), suggesting that the presence of Pouteria annanmensis may have helped to improve seasonal resource availability in resourcepoor months. OGF had higher habitat-exclusive food species richness (7) than MSF (3) and YSF (1), although the differences of species dominance between various successional stages were not significant. It is crucial to consider these food species, especially in OGF, as a priority for conservation to develop effective management and restoration plans for Hainan gibbons.

Implications for conservation
This study allows informed decision-making regarding what forest types this highly endangered species needs to survive (Tweheyo et al. 2004). At present, suitable forest habitat for Hainan gibbons is fragmented and mainly restricted to a very small region of ca.14-16 km 2 (Zhou et al. 2005). Fellowes et al. (2008) revealed that gibbons frequently forage in narrow strips of habitats adjacent to degraded shrub lands and monoculture pine plantations. So one promising activity would be to expand the area suitable for the use by gibbons through planting food tree species in shrubgrass lands, old logging tracks and skid rows as well as plantations where food species diversity is currently the lowest (Table 1). Suitable species include those consumed most frequently by the apes such as Endospermum chinense, Bischoffia javanica, Polyalthia laui, Ficus championii and Nephelium topengii and species which fruit in late winter months when the natural food source is scarce, especially Schefflera octophylla, Ficus altissima, Elaeocarpus subglobosus, Podocarpus neriifolius and Pouteria annanmensis (Chen et al. 2009;Wu 2007). Other studies have shown that primates prefer primary forest (Fan et al. 2011;Furuichi et al. 2001;Zhang et al. 2010). However, with increasing recovery time, secondary forests may help the survival of Hainan gibbons and improve habitat suitability because of the high density and richness of important adult foods trees (Fig. 3). We believe this approach to be very promising because the preference for certain vegetation types is influenced mainly by the abundance of food resource rather than forest type per se (Furuichi et al. 2001). Pinus merkusii plantation forests could also serve as corridors between suitable forest patches or as refuge vegetation for Hainan gibbons if high-quality food species were not removed during low-impact logging or supplemental planting programmes were adopted. Accordingly, future priority actions in BNNR have to be aimed at (1) further protecting large food trees and habitat-exclusive food species in each habitat; (2) promoting recovery and regeneration of natural forest fragments and plantations where most of the important food tree species are associated; (3) implementing special measures (such as directional felling) to avoid damaging food trees and reducing the degree of forest canopy damage so as to preserve important food resources and nesting sites for Hainan gibbons and (4) extending our research to predict the potential distributions of important food tree species and analysing the effect of food availability on Hainan gibbon. Of courses, besides food availability, other factors, such as genetic and anthropogenic factors, are also likely to constrain Hainan gibbon population (Bryant et al. 2016;Bryant et al. 2017) and conservation actions should also focus on managing these factors. If these measures are adopted in the near future, there is a realistic chance that this rare species can survive and remain amongst those species that are characteristic to Hainan Island.

Conclusions
Intermediate secondary forests, which appeared to have a more similar food species composition to old forest and higher adult food species diversity, could serve as an alternative habitat for Hainan gibbons in the short term, given that further anthropogenic disturbances are successfully ruled out. However, Hainan gibbons will find limited resources in secondary forests less than 25 years old, while non-suitable structural characteristics and limited food resources make Pinus merkusii plantations unsuitable for Hainan gibbons.

Behavioral responses of Australian fur seals to boat Abstract
In Australia, a multi-million-dollar industry is based on viewing the Australian fur seal (Arctocephalus pusillus doriferus), predominantly through boat visits to breeding colonies. Regulation of boat approaches varies by site and no systematic investigations have been performed to inform management guidelines. To investigate possible effects of disturbance, experimental boat approaches were made to a colony at Kanowna Island in northern Bass Strait and seal responses were monitored using instantaneous scan sampling. Colony attendance (individuals remaining ashore) was found to be influenced by approach distance and time of day, but was not affected by environmental variables or season, whereas onshore resting behavior was influenced by approach distance, time of day, ambient temperature and wind direction. Onshore resting behavior decreased following experimental boat approaches to 75 m, but changes in abundance of individuals ashore were not observed at this distance. In contrast, approaches to 25 m elicited a strong response, with a steep decline in the number of individuals ashore. This response was strongest when approaches occurred in the morning, with a decline of approximately 47% of individuals, compared to a decline of 21% during afternoon approaches. With regard to onshore resting behavior, afternoon approaches to 75 m led to minimal response. The remaining three combinations of approach distance and time of day had a similar pattern of reductions in the proportion of individuals engaging in onshore resting behavior. The strongest response was again seen during approaches to 25 m conducted in the morning. These behavior changes suggest that unrestricted boat-based ecotourism at Australian fur seal colonies has the potential to increase energy expenditure and reduce the number of seals ashore. Increasing minimum approach distances to ≥75 m and/or restricting visits to afternoons may minimize these impacts at Kanowna Island during the post-molt and non-breeding seasons. As several studies have demonstrated considerable

Introduction
Despite their association with the marine environment, pinnipeds must haul-out on land or ice to rest, evade marine predators, and molt. In addition, they give birth and nurse their young ashore (Gentry and Kooyman 1986;Riedman 1990). While ashore, seals utilize sight, smell and hearing to detect potential above-water threats and to perceive the level of risk these threats present to their survival (Frid and Dill 2002;Nordstrom 2002). When a perceived above-water threat is detected, and individuals deem this threat to be significant, seals will respond by fleeing to the relative safety of the water (Cowling et al. 2015). In gregarious pinnipeds, individuals often detect threats and perceive risk based on the responses of their neighbors, leading to large-scale cascading responses in densely occupied areas that can result in injuries and have significant impacts on colony attendance (Barton et al. 1998;Jay et al. 1998;Stirling 1972).
Pinniped-based ecotourism activities make use of seal haul-out behavior to observe individuals in their natural habitat (Kirkwood et al. 2003). Although the purpose of ecotourism is to give patrons the opportunity to observe animals in the wild engaging in typical behaviors (Orams 1995), ecotourism-based human interactions may instead alter pinniped behavior by initiating responses indicative of predation risk (Frid and Dill 2002). Such responses can interrupt vital activities, increase energy expenditure, reduce breeding success and even cause injury or death (Boren et al. 2002;Marmion 1997;Shaughnessy et al. 2008).
Various environmental factors have the potential to influence a seal's ability to detect threat stimuli and, thus, may also affect responses to anthropogenic disturbance. Wind strength and direction can greatly affect both olfactory and auditory detection of a potential threat (Riedman 1990). Also, sea conditions may influence the willingness of individuals to return to the water, potentially confounding any response to boatbased disturbance. In addition, distance, speed and direction of anthropogenic approach can be important factors in perceived risk, further influencing a seal's response (Frid and Dill 2002;Shaughnessy et al. 2008).
The Australian fur seal (Arctocephalus pusillus doriferus) is one of the world's least abundant fur seals and is endemic to Bass Strait, southeastern Australia (Kirkwood et al. 2010). Ecotourism operators visit Australian fur seals at numerous locations, including haul-out sites (where few pup births occur) and breeding colonies (Kirkwood et al. 2003). Previous research on Australian fur seals at haul-outs has indicated that seals there are particularly sensitive to both land-and boat-based approaches (Burleigh et al. 2008;Shaughnessy et al. 2008). When on land, the seals perceive threats by sight, sound and smell, and generally respond by becoming alert (changing posture and looking toward the threat) and moving to the water (Kirkwood and Dickie 2006;Shaughnessy et al. 2008). Responses may be more pronounced at breeding colonies than at haul-outs, as haul-outs contain mainly adult/sub-adult male and juvenile seals, whereas colonies are densely populated with adult females and pups, groups that are particularly sensitive to disturbance (Barton et al. 1998;Boren et al. 2002;Holcomb et al. 2009). In addition, in colonies where seals have access to higher elevations individuals might perceive their distance to water as greater than seals on flat terrain, and thus sense a greater threat from the boat's approach. They may also detect the boat from a greater distance. Research on New Zealand fur seals has indicated that early detection of a threat tends to moderate responses (Boren 2001).
Despite a multi-million-dollar ecotourism industry based on visits to Australian fur seal breeding colonies (Kirkwood et al. 2003), descriptive and observational data have suggested a heightened response to above-water threats in this species. However, only limited conclusions can be drawn from such data. At the time of undertaking this study management guidelines at Australian fur seal breeding colonies varied by site and season, with minimum approach regulations ranging widely from 30 m (yearround, Seal Rocks) to between 50 m (non-breeding season, Kanowna Island) and 200 m (breeding season, Kanowna Island). These guidelines are largely based on anecdotal evidence. Needlessly strict guidelines may harm ecotour businesses, while unrestricted access could negatively affect seal populations. A recent study on the effects of boatbased approaches on New Zealand fur seals highlighted the ability of evidence-based studies to inform mitigation guidelines (Cowling et al. 2015). However, before evidence-based mitigation guidelines can be proposed for Australian fur seals, controlled studies, which account for potential confounding factors, must be conducted to assess species-specific responses to boat approaches. Therefore, the aims of the present study were to determine seal attendance and behavior in relation to controlled boat approaches at a breeding colony while accounting for potentially confounding environmental variables.

Study sites and field procedures
This study was conducted at an Australian fur seal breeding colony on Kanowna Island (39°10'S, 146°18'E; Fig. 1) in northern Bass Strait. The island is a granite outcrop with steep cliffs and tussock grass vegetation. The fur seal colony has an annual pup production of approximately 3000 (Gibbens and Arnould 2009;Kirkwood et al. 2005) and is dispersed around two main breeding areas, the larger Main Colony (north-western coast) and the smaller East Colony (east coast). Pups are born from early November to mid-December, molt their natal coat in March and are suckled through to October, with breeding females present at the colony year-round. Terrestrial access to Kanowna Island is restricted to researchers and National Park rangers (several visits a year) and the waters adjacent to the island are within a designated marine park. The then current management guidelines prohibited boats from approaching the island to less than 50 m in the post-molt period (March-October) and less than 200 m during the breeding/post-breeding period (November-February, Patkin 2005). A single commercial ecotour operator conducts 10 or less visits per year to the Main Colony, and the waters surrounding the island are infrequently used by recreational boats and fishing charters (both colonies may be exposed to these infrequent visits). Research-related boat landings occur around 10 times per year, along either the west coast or adjacent to the East Colony, with researchers occupying a plateau-area field camp away from the seals for several months each year.
Experimental boat approaches were conducted at the East Colony. The study area comprised a single stretch of smooth granite shelf where seals haul-out less than12 m from the water's edge. The shelf slopes upward toward the island's interior, with seals resting 1 to 10 m above sea level. The site had a concealed vantage point at the northern end, 14 m above the study area's highest point, from which an observer had an unobstructed view of all seals present. Two digital video cameras (GZ-MG630 and GZ-MG680, JVC, Yokosuka, Japan) were placed at this vantage point to record colony attendance and seal behavior for 30 min prior to, 15 min during and 60 min after each experimental boat approach. The number of seals descending from the hinterland of the East Colony (a portion of the island's elevated interior approximately 10-200 m from the shore that the seals access using two narrow paths) was counted by an observer during the same sampling periods.
Experimental approaches were conducted by three boats ranging in length from 5.4 to 10 m, with the majority (87%) conducted by the 10 m vessel ( Table 1). The experimental approach procedures were developed based on a mix of the then current regulatory guidelines and the approach technique at other non-naive colonies/ haulouts in the region. Logistic constraints limited this study to two approach distances, initially set at 25 m and 50 m from the shoreline, 50 m being the then current minimum approach protocol. A pilot study revealed strong responses at 50 m and, consequently, the distances were changed to 25m and 75 m. Thus, the experimental distances bracketed the current management guidelines. On each approach, the boat moved at a speed of 25 kn (50 km·h -1 ) from an out-of-sight point, came into sight 1.5 km offshore toward the study area, slowed to 5 kn (10 km·h -1 ) at 500 m and proceeded to one of the two approach distances (selected at random), where it remained for 15 min then departed on the same path. A handheld range finder (accuracy ± 0.91 m, Bushnell Yardage Pro Sport 450, Overland Park, KS USA) onboard each boat ensured approach distances were accurate. Time of day, meteorological conditions including ambient temperature, wind speed and wind direction (onshore or offshore) and sea state (i.e. calm, choppy, rough) were recorded for each approach. Wind speed and direction were subsequently transformed to a linear vector where values below 0 represented the speed of offshore wind and values above 0 represent the speed of onshore winds, hereafter referred to as direction-adjusted wind speed.
Seals in the study area were left for greater than 4 h to recover/redistribute between approaches and a maximum of two approaches was conducted per day (one in the morning between 08:00-11:00; one in the afternoon between 13:00-16:00). On average, there were 3.4 ± 0.9 days between experimental approach days. Approaches were conducted during two sampling periods: summer post-breeding (January-February), when most adult and sub-adult males have dispersed and the colony is occupied primarily by 2-3-month-old pups, adult females and juveniles; and winter post-molt (May-August), when the colony comprises mainly 6-7-month-old pups, adult females and juveniles (Warneke and Shaughnessy 1985).

Data collection
Digital video recording enabled observations of a large number of individuals simultaneously. Videos were analyzed using instantaneous scan sampling (Altmann 1974; Boren et al. 2002) at 5-min intervals pre-approach, at 1-min intervals while the boat was at the designated approach distance (to maximise our ability to detect rapid behavioural changes occurring during the peak disturbance times), and at 5-min intervals for 60 min post-approach (26 scans per boat approach). Boat arrival was defined as T= 00; thus, boat departure occurred at T= 15 and the end of the sampling period at T= 75.
During 11 trials, observations ended prematurely due to events such as inclement weather, an unrelated boat within sight of the colony or equipment failure. During each instantaneous scan, the age/sex, posture and behavior of each individual within the study area was recorded. Age/sex classes comprised adult males, adult females, sub-adult males, juveniles and pups, and were based upon Goldsworthy and Shaughnessy (1994), with classes differentiated based on size, head and snout shape, shoulder development and pelage ( Table 2). The ethogram used to classify behavior was modified from that used by Boren et al. (2002) for instantaneous scan samples at colonies of New Zealand fur seals (A. forsteri). Seals were classified into one of five behavior classes: At Rest, Comfort, Active, Mother-pup or Interaction (Table 3). Due to the rare occurrence of Interaction behaviors (less than 1 per scan), these were included in the Active category. To quantify colony-wide behavioral responses to vessel approaches, the above behavior classes were incorporated into one of two categories, Resting (including At Rest, Comfort, and Mother-pup) or Active (Boren et al. 2002).
Pre-approach age/sex composition, attendance and behavior were calculated from the averaged pre-approach scans for each trial. Due to the variation in attendance numbers between trials, for inter-trial comparisons all data were converted into a proportion of pre-approach attendance. For the same reason, the number of Resting seals ashore (versus Active) was also converted into a proportion of the number of seals ashore at baseline for inter-trial comparisons. Recovery time was measured at the colony level based upon the time it took for the study area to return to pre-approach numbers and activity levels (as a proportion of pre-approach numbers). The use of pre-approach attendance numbers to measure recovery assumes that the same number of seals that entered the water in response to approach returned to the study area however because we do not have information on individual identity of seals we only make our inferences at the colony level and do not assume that the same seals present pre-disturbance are those returning post-disturbance. Studies on other pinnipeds indicate that adult females exhibit extremely high site fidelity and generally do return to the same area following a benign disturbance (Kovacs and Innes 1990;Lidgard 1996). Accordingly, and as individuals that fled could not be tracked individually, it was assumed that they remained in the water until they returned to the study area. Additionally, the proportion of seals that remained ashore during each scan that were in the Resting category (versus the Active category) were analyzed. This provided activity level data that only took into account the seals that remained ashore during each scan as a further measure of colony recovery.  Goldsworthy and Shaughnessy (1994).

Age/sex class Description Adult males
Mature males, possessing well-developed chests, manes and shoulders. Sub-adult males Males, similar in size or slightly larger than adult females but distinguished from them by shoulder development, larger head and pointed snout. Adult females Mid-sized animals with smaller, sleeker heads than males and lacking shoulder development.

Juveniles
Smaller than adult females and sub-adult males but larger than molted pups, and with more defined muzzles and muscle tone. Pups Seals <1 year old, with black natal pelage until March/April, then molting to silver-grey juveniletype pelage. Table 3. Ethogram of seal behavior modified from Boren et al. (2002).
Lying with head down, lying or sitting with head arched in "skypointing" position, usually with eyes closed Comfort Grooming, scratching, shifting position, flipper-waving and other thermoregulatory behaviors Mother-Pup In pups, suckling; in females, nursing and/or sniffing, caressing a suckling pup Active Lying or sitting with head up and aware, alert or moving Interaction Interaction with another animal (excluding mother-pup pairs), noted if aggressive.
In addition to the seals within the study area, the movements of individuals between the hinterland above it and the shoreline were monitored to determine whether seals up to 100 m inland detected and were disturbed by experimental boat approaches. Seals accessed the hinterland via two narrow pathways that were easily visible from the observation point.

Statistical analysis
Behavioral observations were correlated temporally, while also being nested into individual boat approaches (i.e. multiple serial observations per boat approach). There was no expectation that the response of seals to boat approaches should show a linear response. To account for this nested structure and the expected non-linear response, data were analyzed using Generalized Additive Mixed-effect Models (GAMMs). GAMMs were fitted to investigate the response of the number of animals ashore (Poisson error distribution with log link) and the proportion of resting individuals ashore (binomial error distribution with log link) by treating each instantaneous scan sample as a single observation in time. Experimental (boat approach distance and time of day, time since disturbance) and environmental (breeding cycle period, ambient temperature, direction adjusted wind speed, sea state) fixed predictor variables were included in the models, with a random effect of unique boat approaches. Behavioural scans (attendance and activity level) were included as response variables.
Within these models, smoothing splines (thin plate regression splines) were fitted to the time since disturbance and ambient temperature at the time of observation. Degree of smoothness was calculated via Generalised Cross Validation following (Wood 2006). The interaction between time of day, approach distance and the time since disturbance was modelled using splines fitted to the time since disturbance split by experimental factors such that individual splines were fitted for each of the four unique factor combinations of time of day (AM or PM) and approach distance (25 m or 75 m). Additional environmental variables were treated as factors (sea state and season) or linear with second-degree polynomial (direction adjusted wind speed) covariates. Residual temporal correlations between observations were accounted for (to the best of our ability) using an auto-regressive correlation structure of the order 1 [corAR1(form= ~time|ApproachTrial. no)] within the GAMM models. The most appropriate model -with or without correlation structure -was selected for via AIC. The most parsimonious models were identified using a backwards step-wise AIC based model selection. All models were fitted using the mgcv package (v 1.8-7;Wood 2006) in the R statistical framework.
To investigate the differences in baseline attendance and activity level during periods exclusive of boat approaches paired and student t-tests were used, following assessments for normality and transformation of data where necessary. Non-parametric data were compared using Mann-Whitney U tests or Wilcoxon signed-ranked tests. Unless otherwise noted, all references to proportional change regarding behavioral data were quantified relative to pre-approach values. Data are presented as the mean ± one standard error (SE) and results are considered significant at P < 0.05.

The effects of boat approach on colony attendance
In the summer post-breeding period, a mean of 110.0 ± 12.7 seals were present in the study area prior to approaches, comprising 33.1 ± 1.9% pups, 38.8 ± 1.9% adult females, and 28.1 ± 1.4% other seals (primarily juveniles, occasional sub-adult males and rare adult males; Fig. 2). In the winter post-molt period, 146.9 ± 8.5 seals were present pre-approach, comprising 74.0 ± 2.3% pups, 23.5 ± 2.0% adult females, and 2.5 ± 0.4% other seals (Fig. 2). During the summer post-breeding period, the number of seals present pre-approach was lower during the afternoon than in the morning (t 6 = 2.627, P = 0.039). Time of day had no effect on pre-approach seal numbers during the winter post-molt period (t 6 = -2.194, P = 0.071).
Thirty-eight experimental boat approaches were conducted (18 to 25 m and 20 to 75 m). Approaches occurred between January and September with approaches ranging between 1 and 44 days apart (average: 6.2 d). Model selection on GAMMs fitted to assess the response of colony attendance to boat disturbance resulted in a final model showing the influence of approach distance and time of day to colony attendance, but no effect from the environmental variables or season (Table 4). Approach distance Table 4. Summary results for GAMMs assessing the number of individuals ashore and the proportion of individuals resting in response to boat approaches and weather conditions. Results shown here are for the most parsimonious models selected via AIC-based model selection.   provided the strongest influence on colony attendance to boat approaches (Fig. 3, Table 4). Experimental boat approaches to 75 m showed no change in the trend in abundance of individuals ashore, with there being either a slight but steady decrease (in the morning) or increase (in the afternoon) in numbers through the observation period (Fig. 3, Table 4). In contrast, approaches to 25 m elicited a strong response, with a steep decline in the number of individuals ashore beginning shortly before the boat arrived at the 25 m point and continuing until the boat had left (Fig. 3, Table 4). This response was strongest when approaches occurred in the morning, with a decline of approximately 47% of individuals from the beginning of the boat's arrival until it left, compared to a decline of 21% when the approach was conducted in the afternoon. There were sudden increases in the number of seals descending from the hinterland to the shore during eight (21.1%) experimental boat approaches, six in summer and two in winter (Descents by approach type and time of day: 25 m/AM = 1, 25 m/PM =2, 75 m/AM = 5). This suggests that seals distant from the shoreline study area occasionally detected the presence of the boat and perceived it as a threat.
Model selection on GAMMs fitted to investigate the influence of boat disturbance to the proportions of resting individuals ashore resulted in a final model demonstrating the influence of approach distance and time of day on colony attendance, as well as that of ambient temperature and direction adjusted wind speed (Table 4). Approaches conducted in the afternoon to 75 m showed an almost flat response with regard to the proportion of individuals resting during boat approaches (Fig. 4). The remaining three combinations of approach distance and time of day had a similar pattern in terms of the proportion of individuals resting, with reductions in the behaviours of individuals remaining ashore beginning between 5 and 10 min prior to the boat's arrival at the prescribed distance and continuing until approximately 5 min after arrival (Fig. 4). The strongest response, determined by the steepness of decline during the disturbance period, was again seen during approaches to 25 m conducted in the morning, where the proportion of individuals remaining ashore who were resting did not recover by the conclusion of the recording period, as evidenced by the numbers of resting individuals not returning to pre-disturbance levels during the monitoring period (Fig. 4).
Weather was also found to have an effect on the proportion of individuals resting with fewer individuals resting as the temperature increased (Fig. 5a, Table 4). There was also a slight linear reduction in the proportion of individuals resting as winds turned onshore and increased in speed (Fig. 5b, Table 4).

Discussion
Australian fur seals at the Kanowna Island colony responded to experimental boat approaches to both 75 and 25 m by becoming more active ashore and responded to 25 m approaches by fleeing to the water. Of the factors examined in the present study, proximity of approach had the greatest influence on the ability of seals to detect a boat and perceive it to be a threat. Time of day also influenced the strength of responses (determined by the steepness of decline during the disturbance period) to boat approaches, while environmental factors (temperature, wind speed/direction) were only observed to affect baseline resting behavior.
Australian fur seals, like other fur seals, are thought to rely on olfaction or auditory stimuli as principle means of detecting threats when they are on land (Riedman 1990;Shaughnessy et al. 1999). This is because they are short-sighted when out of water and cannot focus on objects in the distance (Hanke et al. 2006). The results of the present study, in which individuals were disturbed by boats up to 75 m away, suggest initial detection was by olfaction or audition. However, although individuals ashore displayed increased alert behaviors, when boats remained at 75 m this did not translate into animals leaving the colony. In contrast, when boats approached to 25 m, departures occurred. As such, it appears that while olfaction/auditory stimuli play a role in the initial detection of a potential threat, individuals may not respond further until a threat is perceived to be imminent (Tripovich et al. 2012).
The most severe response Australian fur seals exhibited to boat approaches was to flee toward the water; during approaches to 25 m, this caused dramatic changes in colony attendance. The periods fur seals spend ashore at colonies are particularly important for resting, evading predators, molting, breeding and rearing young (Gentry and Kooyman 1986;Riedman 1990). Fleeing behaviors in themselves expend energy, and time spent in the water as a result of flight responses can also be energetically costly due to active movement, being alert for predators and maintaining body temperature (Donohue et al. 2000). The actual energetic costs related to fleeing behaviors in this species are difficult to determine, but could result in the need for additional rest time ashore and/or greater caloric intake.
Fleeing responses to boat approaches may have even greater implications for adult females and pups, which are of particular concern because they are bound to breeding colonies and may be more energetically constrained than other age classes (Barton et al. 1998;Boren et al. 2002;Riedman 1990). Pups are also at risk of being trampled or falling from cliffs as a result of stampedes (Mattlin 1978). In addition, if they flee into the water they would incur high energetic costs and experience greater risk of predation, as they are small and have limited swimming and diving capabilities (Donohue et al. 2000;Shaughnessy et al. 1999). The separation of mother-pup pairs and resulting interruptions to suckling behavior also have serious implications for pup survival (Boren 2001;Kovacs and Innes 1990).
Although boat approaches to 75 m did not significantly impact seal attendance, the observation that approaches to 75 m resulted in significant behavioral changes still has substantial implications. Many seals remaining ashore changed posture and stayed alert for the entire duration of the boat's visit, investing more time and potentially more energy into vigilance behavior (Amo et al. 2006;Lidgard 1996) than they did before the boat arrived. Thus, events that appear to cause low-level disturbance without large-scale changes in attendance at Australian fur seal colonies could influence a seal's energy budget. Further investigation of the physiological impacts of disturbance on this species is warranted.
Time of day was also found to influence seal responses to boat disturbance, with individuals showing stronger responses to disturbance during morning approaches. It is possible that the weaker response observed in the afternoon was biased by afternoon trials being preceded by a morning approach, which may have habituated the seals present for the earlier approach or displaced sensitive seals to another location at the colony (Nordstrom 2002). Two approaches occurred on the same day for 80.0% of trials. Further investigation into decreased sensitivity with repeated disturbance, as well as potential temporal and density-dependent responses, is needed to fully explain these results. In relation to vessel approaches, greater caution or greater approach distances may be required at the first approach compared with latter approaches.
While colony-level responses to boat approaches were not found to differ between seasons, ambient temperature did have an effect on baseline resting behavior, with individuals more likely to be resting at cooler temperatures. However, no effect was identified with regard to the number of individuals remaining ashore during boat approaches, suggesting that temperature did not influence individual responses to threat stimuli. Cowling et al. (2015) identified a seasonal influence on seal response to such stimuli, with the proportion of New Zealand fur seals at rest decreasing during summer months. However, Cowling et al. (2015) did not measure ambient conditions directly but rather used time of year as a proxy for weather so it is not clear whether their results were temperature-related or due to other factors, such as breeding cycle. Further studies on the influence of temperature on fur seal behavior, particularly in terms of resting behavior and disturbance, is needed to clarify these relationships.
For many pinniped species, colonies are not restricted to the immediate shoreline and can often extend considerable distances and/or elevations inland (Stevens and Boness 2003). Seals resting inland may be less likely to detect an approaching boat and, because proximity to the water provides a level of security for some seals (van Polanen Petel 2005), they may perceive a higher threat level and, consequently, respond more severely than seals near the water's edge. Hence, disturbances that affect animals inland at colonies may result in extra energy expenditure due to the greater distances animals must travel to get to water, and may cause injury and even death if stampeding animals fall over cliffs or into crevasses (Mattlin 1978).
In the present study, the number of seals descending from the hinterland above the study area increased during 21.1% of experimental approaches. Three stampedes also occurred there as a result of boat approaches. During one such stampede (the final of the three) in the non-breeding season, two pups were observed to fall >10 m from a cliff and land on rock shelves below. Two of these stampedes occurred during afternoon approaches to 25 m, while one occurred during a morning approach to 75 m. Similar events were observed at a New Zealand fur seal colony due to land-based research activities, with stampedes causing several pups to be trampled to death and another pup dying as a result of falling 10 m off a cliff (Mattlin 1978). Direct injuries or mortalities to pups resulting from boat approaches have not been observed in Australian fur seals or other pinniped species (Boren et al. 2002;Burleigh et al. 2008;Nordstrom 2002;Shaughnessy et al. 2008). Due to the relative rarity of these events, we were unable to determine which approach distance or environmental factors (e.g. wind direction) influenced hinterland stampedes. Nonetheless, these findings do suggest that guidelines at the time of this study allowing 50-m approach distances at Kanowna Island would be unlikely to entirely preclude such events. It is possible that limiting approaches to ≥75 m and to afternoons could minimise hinterland stampedes but further research would be needed to confirm this. While the infrequency of these events suggests they are unlikely to have population-level effects, such disturbance impacts are in violation of state and federal regulations protecting marine mammals (IUCN 2007;Patkin 2005). Thus, the terrain of colony and haul-out sites, particularly seal use of hinterland/elevated regions, should be taken into consideration when developing management guidelines. Haul-out sites or low-lying colonies with limited access to such areas may not need to be managed as stringently as colonies with greater variation in terrain.
There were several limitations to this study. In addition to inter-specific variation among pinnipeds, several previous studies demonstrated a high degree of intra-species variation due to factors including age/sex class present, site type, density ashore and previous exposure to humans (Barton et al. 1998;Boren et al. 2002;Shaughnessy et al. 2008). Hence, developing appropriate management guidelines is challenging without extensive species-specific investigation. As this study was conducted at a single site, interpretation of data with regard to site structure, terrain, population density and previous exposure to vessel traffic must be made with caution, limiting the generalizability of these findings to other breeding sites. In addition, analysis based on boat size and vessel-related noise could not be conducted due to logistical constraints.
Furthermore, colony attendance and onshore behaviour were treated as independent predictor variables within two separate GAMMs. There exists the possibility that different behavioral phenotypes within the population meant that only a certain class of individuals remained onshore during boat approaches -the bold individuals -as such these data may not be truly independent. Model complexity (nested data structure with non-linear responses) excluded us from examining these data using multivariate methods. Therefore, data should be interpreted with the understanding that onshore behavioural responses represent the subset of the population that remained onshore and that this subset may not be a truly random sample of the population as a whole.
In summary, the findings of the present study reveal distinct gradients of response to the approach of boats at an Australian fur seal colony influenced by both approach distance and time of day. These results suggest that unrestricted boat-based ecotourism at colonies, particularly approaches to <75 m, may have implications for energy expenditure and reproductive success in the Australian fur seal. Current guidelines, implemented following this study, now limit boat approaches to 100 m at Kanowna Island from March through October now making them unlikely to affect behaviour at this colony. However, this study was limited to a single colony and further research will be needed to generalize these results to other sites.

Introduction
Many environmental threats, including climate change, biodiversity loss and degradation of land and freshwater are attributed to unsustainable agricultural practices (Foley et al. 2011). Agricultural intensification worldwide has had tremendous impacts on habitats, biodiversity and ecosystem functions and services (Zhang et al. 2007, Power 2010, Foley et al. 2011. The success of key environmental policy and legislation related to climate change targets (e.g. Roadmap to 2050) and meeting the 2020 biodiversity targets are heavily dependent on agricultural activities. In fact, Target 3 of the EU 2020 Biodiversity Strategy states explicitly that, by 2020, measures must be taken to "bring about a measurable improvement in the conservation status of species and habitats that depend on or are affected by agriculture and in the provision of ecosystem services as compared to the EU2010 Baseline" (http://www.alter-net.info/). Agroecosystems are often considered incompatible with conservation and protection. However, in Europe the growing recognition that biodiversity conservation depends on the persistence of low-intensity farming systems catalysed the formation of the concept of High Nature Value farmland (HNVf ) (Beaufoy and Jones 2012). HNVfs are parcels of farmland that are expected to support high levels of biodiversity or species and habitats of conservation concern because of certain features. HNVfs are classified into three categories, depending on the presence of semi-natural vegetation in the farmland, the diversity of low input farmland and the presence of species of conservation concern (Table 1) (Paracchini et al. 2008;Beaufoy and Jones 2012).
HNVfs are present all over Europe and include semi-natural pastures and hay meadows in Western Europe, Ireland and the UK and traditional cropping and orchard systems where livestock activities can also be present, such as olive groves, vineyards and cork oaks in the Mediterranean (Oppermann et al. 2012;Keenleyside et al. 2014). The largest extent of HNV farmland is found in Southern Europe with Spain having the highest percentage of HNVfs with respect to its utilised agricultural area (Keenleyside et al. 2014). HNVfs offer a great range of ecosystem services in addition to food provision, such as maintenance of crop genetic diversity, biocontrol services, erosion control, pollination and recreation services, while at the same time, they support rural livelihoods in remote areas (Keenleyside et al. 2014). HNVf landscapes are important for their ecological, cultural and socio-economic role. Some systems, such as the dehesas or montados in Spain and Portugal or the vineyards producing the sweet wine Commandaria on the south-western foothills of Cyprus, are historic landscapes. HNVfs have the potential to offer economic benefits and sustain rural populations as well as create benefits for the agrotourism industry and the wider public (O'Rourke et al. 2016).
The percentage of HNVf relative to the total utilised agricultural area (UAA) is acknowledged as an important indicator for monitoring landscape changes at spatial and temporal scales, as it is sensitive to anthropogenic impacts and it can be easily communicated in a policy-relevant context. The identification of HNVf in EU member states (MS) and monitoring its extent and condition is a policy instrument for agriculture and nature conservation purposes. HNVf is one of the 32 agri-environmental indicators developed by Eurostat to monitor the effects of agriculture on the environment and one of the impact indicators in the Common Monitoring and Evaluation Framework of CAP (http://epp.eurostat.ec.europa.eu/). At the same time, it is also a series of practical measures for biodiversity conservation (i.e. as habitats) and for providing valuable connections in the landscape for species movement (Doxa et al. 2012;Sutherland et al. 2010). For instance, vineyards and carob groves in Cyprus, provide feeding and breeding grounds to a diverse range of organisms, including many species listed in the Birds and Habitats Directives, such as the endemic birds, Cyprus warbler (Sylvia melanothorax) and Cyprus wheateater (Oenanthe cypriaca), as well as raptors, including the kestrel (Falco tinunculus) and the long-legged buzzard (Buteo rufinus). Along with bird species, many key insect species such as the critically endangered beetle Propomacrus cypriacus and a diverse range of pollinators, reside in agricultural fields and their margins (Sparrow and John 2016).
Regrettably, HNVf systems are currently threatened by two major opposing forces, intensification on the one hand (e.g. high inputs of pesticides or overgrazing) and land abandonment on the other. On the Mediterranean island of Cyprus, both forces have been acting simultaneously on the landscape for the last 50 years. The Utilized Agricultural Area (UAA) has seen a dramatic decline (Cyprus Statistical Service 2012) with extensive cereal crops and vineyards experiencing the largest decline. Intensification and land abandonment stem from the farmers' need to increase their income which is usually a response to political and economic imperatives. Especially in marginal areas, a farmer's profit can be discouragingly low for managing the land in a traditional, sustainable manner. The current challenge is to identify and apply strategies which will maintain the ecological benefits and create further socioeconomic turnover for the stewards of the HNV farming land across Europe. An essential first step towards this goal is defining what constitutes an HNVf system, describing its particular features and mapping HNVf land for each EU MS.
The methodology proposed by Andersen et al. (2003), revised by Paracchini et al. (2006), provides a European HNVf map which is not sufficiently accurate for applications at national and regional levels. Therefore MS are advised to use this indicative method (Beaufoy and Cooper 2008) and apply at national level, by adapting it to their specific conditions and contexts (EEA 2016). To this end, there have been various attempts to HNVf mapping (Hazeu et al. 2014;Halada et al. 2011;Galdenzi et al. 2012). In the case of Types I and II, two complementary approaches have been employed. The first is based on the use of land cover maps, which is suited for the localisation of HNV farmland areas. The second is the farm system typology, which combines agronomic and economic data derived from farms (e.g. FADN). The methodology employed at various national attempts follows or relies principally on the one proposed at the European level by Paracchini et al. (2006). Most attempts at national level use coarse data and therefore only approximate the spatial location of HNVfs. In a paper by Lomba et al. (2014), the caveats associated with the identification, mapping and monitoring of HNVfs in Europe are discussed, many of which stem from the very definition and concepts associated with HNVfs. The various attempts so far demonstrate that the methodology proposed is 'confusing' and is interpreted in different ways when applied (Lomba et al. 2014;Lomba et al. 2015).
Agriculture in Cyprus is characterised by a mosaic pattern of small plots with a diversity of permanent and arable crops. Two thirds of the farms in Cyprus are less than 20 ha in size (FSS 2007) while extensive farmlands with olives, carobs and vineyards have resulted in a diversity of cultural landscapes (Vogiatzakis and Manolaki 2017). This in turn provides an invaluable habitat to a rich biodiversity and, in particular, avifauna found within and around Natura2000 sites (Hellicar 2012). Yet and despite their importance for natural and cultural heritage, to date the extent of these High Nature Value farmlands on the island has not been sufficiently mapped. Therefore the objectives of the paper are to a) characterise and map HNVfs in Cyprus, adapting where necessary existing methodologies and b) identify potential synergies of HNVfs and the current spatial configuration of Natura 2000 areas.

Study area
Cyprus is the third-largest Mediterranean island and a biodiversity hot-spot with high levels of endemism (Tsintides et al. 2007). Agricultural fields cover ca. 50% of the land surface of the island in which the Government of the Republic of Cyprus exercises full control (CYGCA). Mapping was confined to these areas (CYGCA) since there were no data available for the northern part of the island.-In the last 40 years, the Utilized Agricultural Area (UAA) of Cyprus has seen a dramatic decline, from approximately 180,000 ha in 1975 to 115,000 in 2012 (Cyprus Statistical Service 2012). The largest acreage decrease was observed for cereals and vines, two extensive crops, which could potentially meet the definition of HNVf. Cereal crops (mostly rain-fed barley) decreased from a maximum of 73,000 ha in 2003 to 38,000 ha in 2012, while vineyards dropped from a maximum of 35,000 ha in 1975 to 7,000 ha in 2012. Carob groves, another extensive crop and an important feature of the island's landscapes, have also seen a dramatic reduction from 7,000 ha in 1985 (earliest available data) to 1,700 ha in 2012. Approximately 10% of the UAA of Cyprus is used for growing intensive crops, mainly vegetables, potatoes, citrus and other fruits.

Datasets and mapping process
Mapping of HNVfs usually relies on four sets of indicators (Lomba et al. 2014): landscape elements, extensive practices, crop diversity and indicator species. We followed a similar approach and employed high resolution data with experts' opinion and fine- Figure 1. Selection rules for mapping the three types of potentially HNVfs in Cyprus using the conservative storyline (25% threshold). The liberal storyline was applied with identical rules with the exception of the 10% threshold for % farmland within every 1 km 2 grid cell. tuning according to ground-truthing and local conditions. High resolution spatial data from the land parcel information system (LPIS) for 2013 (the latest available dataset), containing the delimited plots at national level, were combined with farm system data from the Cyprus Agricultural Payments Organization (CAPO). The CAPO dataset provides information on the type of farm and type of crops cultivated on each plot, including grasslands. Plots, including non-claimed farmland (no financial support requested), were excluded from the analysis, since such farmland was not eligible for payment under CAPO and there is no reporting on crop types or on whether they are active or abandoned. This information along with datasets on farming typology rules, agro-chemical inputs, water use intensity, Natura2000 sites, water bodies and Important Bird Areas (IBAs) were collated in a common spatial framework.
These data were complemented with CORINE Land Cover 2012 data (1:250,000) for the distribution of semi-natural vegetation. Seventeen CORINE level three types were used, namely the categories: 3.1. 1, 3.1.2, 3.1.3, 3.2.1, 3.2.2, 3.2.3, 3.2.4, 3.3.1, 3.3.2, 3.3.3, 3.3.4, 4.1.1., 4.2.1, 4.2.2, 4.2.3, 5.1.2, 5.2.2. CORINE land cover data provided the best proxy information on the distribution pattern of semi-natural vegetation, since national data at similar or finer scales do not include all semi-natural categories and were produced with different schemes and for different purpose (e.g. the state forests map of the island).
The spatial datasets were evaluated for the definition and mapping of HNVfs in Cyprus by combining the specific selection rules applied by JRC for Cyprus (Paracchini et al. 2008; Table 1) with principles of landscape ecology and field work (Fig. 1). A 1 km 2 grid was overlaid on the study area and two scenarios were used for delineating potential areas with HNVf within each grid cell. In the first storyline, farmland occupies at least 10% of total grid cell area (liberal storyline), while in the second, this threshold is raised to 25% (conservative storyline). The first threshold was considered reasonable given the size of agricultural plots on the island (average of ca. 0.5 ha), the fact that not all farmland is registered with LPIS and also considering that many areas are subject to agricultural abandonment. The 25% threshold was considered adequate to characterise grid cells which are predominantly agricultural (i.e. pure agricultural matrix). In each grid cell, the total area of farmland and semi-natural vegetation were tabulated, as well as the share of extensive farmland, the presence of water bodies, the diversity of seminatural vegetation types (used for Type II) and crop diversity. The mapping rules were derived using a dichotomous key with a set of decisive questions which were in line with the definitions of HNVfs while their quantification (threshold setting) was a result of individual consultation with experts and subsequent fine-tuning (Fig.1).
The potential HNV farmland was delimited by combining the three HNV farmland-groups described in the European HNV indicator study of Andersen et al. (2003) (see Table 1; Fig. 2). The spatial relationship of HNVFs with the Natura2000 network on the island was assessed using a simple overlay process.

Stakeholders' involvement
Fourteen organisations were invited to the two workshops. These organisations represented the major stakeholders in land management in Cyprus, i.e. Government, NGOs, Universities and farmers' associations. Twelve organisations attended two workshops held in Nicosia with the numbers of participants (24 in total with complementary expertise) shown in brackets: Departments of Environment (3), Agriculture (2), Forests (3), Water Development (1), the Game and Wildlife Service (2), Birdlife Cyprus (2), Terra Cypria (2), Table 1. HNVf definitions and mapping: A comparison between JRC's approach (Paracchini et al. 2008) and this study.  Figure 2. Extent of potentially HNVfs proposed by JRC (upper part) and potentially HNVfs according to the 10% threshold (middle part) and 25% thresholds in this study (lower part). (2), Management Authority for the Cyprus Rural Development Programme (1), Cyprus Agricultural Payments Organisation -CAPO (2), Cyprus University of Technology -Department of Agricultural Sciences (2), Open University of Cyprus -Terrestrial Ecosystem Management Lab (2).

Cyprus Federation of Environmental and Ecological Organisations
The first workshop brought together, early in the process, stakeholders to discuss the definition of HNVf and the mapping methodology for the Cypriot context. These people, coming from different disciplines and expertise, were introduced for the first time to the proposed rule-based methodology and engaged in a productive dialogue to improve the mapping rules. The rules for mapping were explained and the use of various thresholds of farmland coverage within each grid cell was discussed. Important decisions were made, such as, for example, classifying extensive as opposed to intensive farming activity. More specifically, one hundred crop types were classified into low and high intensity types, based largely on irrigation and agrochemical use intensity. In addition it was agreed that the term farmland was to be used in accordance with other studies (e.g. Plieninger and Bieling 2013), to denote land suitable or used for all kinds of agricultural activities, including grazing where data was available. Stakeholders pointed out that all mapped areas sensitive to nitrate pollution because of intensive fertiliser use / livestock farms, as well as areas where intensive irrigation takes place (mainly the east part of the island), should be excluded from mapping. These areas were also corroborated during selected field visits.
The outcome of the mapping process was presented during the second workshop to the same stakeholders and supplementary amendments were suggested and taken into consideration. During the second workshop, stakeholders pointed out that HNVf Types I and III were underestimated due to the under-representation of grazing lands in the mapping sources used (in both CORINE and CAPO). Conversely, HNVf II was considered overestimated for some areas since, according to stakeholders, the relative contribution of farmlands and natural/structural elements present in the estimation of 'mosaickness' should be somehow weighted with more emphasis placed on farmlands' presence. In addition, it was proposed that the most agriculturally intensive areas in the east of the island could be included in HNVf Type III only where there is confirmed presence of species of European importance (Type III). Stakeholders also commented that known occurrences of HNVfs in mountainous areas were omitted (i.e. areas with extensive presence of traditional stonewalls), particularly under the conservative storyline. The map was updated to reflect the changes proposed by the experts (see Results and Figure 4).

Identification and extent of HNVfs
The areal extent of the two HNVf storylines (25% and 10%) is given in Table 2. Areas containing potentially HNVf may extend from 22.5-34.5% of the island's surface, depending on the threshold employed. HNVf Type II dominates in both scenarios. HNVf Type I is the more restricted type for the 10% threshold (Fig. 3a), while for a b Figure 3. Potentially HNVf types according to a the 10% threshold and b the 25% threshold, in relation with existing Natura2000 sites. the 25% threshold HNVf Type III becomes the most limited type (Fig. 3b). Type II is the most extensive category while Type III is the more restricted in distribution. With a 10% threshold, ¾ of the study area is under one or more HNVf categories. There were also cases where cells could be assigned to more than one HNVf category ( Table  2). The highest overlap was between Type II and Type III for the liberal storyline and between Type I and II for the conservative storyline (Table 2).
During the second workshop, there was a general agreement amongst stakeholders that the map based on the 10% threshold better represented HNVf on the island. However, workshop participants pointed out that some amendments were required to improve HNVf representation; for instance, the exclusion of intensive agricultural lands on the east of the island, unless there was confirmed presence of species of European importance (see Materials and Methods). The incorporation of stakeholders' comments in the mapping process resulted in changes as explicitly shown in Figure 4. This resulted in the final map based on the 10% threshold (Fig. 4) which was adopted by the Ministry for reporting to the EU.

Spatial relation of HNVfs to Natura2000 network
When we applied a 10% threshold for farmland definition within every grid cell, there is a total of 431 km 2 of HNVfs in Cyprus already situated in the Natura2000 network, almost half belonging to cells characterised as both Type II & III. When we applied a 25% threshold, this figure is reduced to 228, although again half belong to cells characterised as both Type II & III. However, in this storyline, there are 43 cells within Natura2000 which are characterised as Type I, II and III (Table 3).   In the case of the liberal storyline (10% threshold), the spatial extent of HNVfs, together with the existing Natura2000 sites, cover 90% of the study area (Figure 3) while, with the conservative storyline (25% threshold), this figure drops to ca. 75% of the study area ( Figure 3). The 10% threshold results in substantial buffering around lowland Natura2000 sites in the central and east part of the island as, for example, Potamos Panagias Stazousas (CY6000007) -Kosiis Pallourokampou (CY6000009) -Alykes Larnakas (CY6000002) sites and more physical connections between Natura2000 in the landscape compared to the 25% threshold ( Figure 4). The only exception is the southeastern corner of the island where there is no connection or buffer between Cavo Gkreko (CY3000005), Limni Paralimniou (CY000008) and Fragma Achnas (CY000007). In the case of 25% threshold, there is no buffering effect around most of the Natura2000 sites of the central and east part of the island, while there is a lack of physical connection between most of the lowland Natura2000 areas, irrespective of the geographical location/district.

Discussion
The identification and mapping of High Nature Value farmlands is of utmost importance in preventing biodiversity loss but also in aiding with their preservation against intensification or abandonment (Andersen et al. 2003;EEA 2004;Paracchini et al. 2008). Although there is no doubt that the standard established procedures (Paracchini et al. 2006) provide guidance, it was obvious from the outset of this study that refinement is necessary at the national level, if HNVf delineation is to be of practical use for policy formulation and planning. On one hand, methodologies should not be prescriptive due to ecological historical and cultural differences across farming landscapes in Europe. On the other hand, the existing approaches are subject to a variety of interpretations and thus are problematic, while the dependency on local conditions and spatial datasets clearly influences the result (see also review by Lomba et al. 2014).
A previous study by JRC classified ca. 37% of the whole island as possible HNVf (Paracchini et al. 2008). However, the present study suggests that this percentage might be even higher given that mapping of the southern part of the island, which is under the effective control of the Republic of Cyprus government, resulted in 22.5% and 34.53% of the island's surface being potential HNVfs (depending on the storyline employed). This is the first study for Cyprus which takes into account expert derived information following a wide consultation and detailed GIS data on crop distribution, in addition to biodiversity datasets (IBAs) in order to classify HNVfs. Despite the availability of crop datasets in very high resolution, they had to be aggregated using a consistent spatial framework. The 1 km 2 grid chosen was small enough to accommodate the average size of farms in Cyprus and large enough to contain additional elements to account and measure heterogeneity in the wider landscapes (mosaickness). The approach is consistent with other HNVf studies (e.g. Kikas et al. 2017;Paracchini et al. 2008). Even so, there are some notable differences in the distribution of HNVfs between our approach and that of JRC's ( Fig. 2), particularly in the southeast and central part of the island, with a better match obtained when using the 10% threshold. As pointed out by the stakeholders involved in the process, the conservative storyline (25%) performs better in the lowlands where agriculture dominates but not in the uplands where agricultural plots are, on average, smaller in size and have fragmented distribution.
Scale is an important issue and the proposed HNVf methodology (Paracchini et al. 2008) should be revised to include different levels of data and different scales of analysis. At the European level, we need a nested hierarchical-based methodology (similar to CORINE), both in terms of spatial and thematic resolution, which countries can adapt according to the level and dataset availability. Traversing from coarse thematic resolution as derived from EU level CORINE, to fine thematic resolution as provided by CAPO in the case of Cyprus, proved problematic since there is no direct relationship between CORINE and CAPO classes and mapping schemes. This clearly demonstrated that the mapping result is directly related to the datasets employed, although conceptually the definition of HNVfs remains constant.
The need for identifying HNVf at the national level is as pressing as the need for clear definitions, e.g. in areas where semi-natural habitats form mosaics with more intensive agriculture. Perhaps from the three HNVf types, Type I is clearer and more consistently mapped across case studies (countries) while is also more easily identified when widespread, especially in marginal regions of Europe (Lomba et al. 2015;Beaufoy and Jones 2012). Type III definition is very clear, as it refers to farmland able to support populations of species of conservation concern, independently of farming intensity (Beaufoy and Jones 2012). Species data availability, particularly useful for characterising HNVf Type III, is greater in the north than in the south of Europe.
Determining Type II systems can be quite subjective as each case-scenario requires some value judgement as to what is the right proportion of semi-natural habitat needed to sustain wildlife populations. Tscharntke et al. (2005) refer to landscapes with more than 20% semi-natural habitats as complex landscapes while those with 2-20% seminatural land are considered as structurally simple. The landscapes with less than 2% of semi-natural land are referred to as cleared landscapes with limited conservation value. Le Roux et al. (2008) consider 20% to be the minimum threshold for maintaining significant biodiversity on farmland (Beaufoy and Jones 2012). Different approaches have also been employed, such as the one by Halada et al. (2011) who propose the inclusion of Annex I habitat types depending on agricultural management, while in Italy, Galdenzi et al. (2012) based their work on current and potential vegetation cover data obtainable from different integrated vegetation maps. The geographical context is evidently important in HNVf definition and mapping with fine grained mosaics characterising southern Europe landscapes compared to their counterparts in the north (Pinto Correia et al. 2018).
The area of agricultural land, classified as HNVf, has been proposed as a key indicator to monitor the impact of CAP of the environment in the 2007-2013 funding cycle and is a prominent metric for the next cycle of RDP (2014-2020). Within Axis 2 of the CAP, the agri-environmental measures are amongst the most important instruments for the support of traditional or HNVf farming (Pe'er et al. 2014). In Cyprus, 43% of the total public 2007-2013 RDP expenditure of 282 million Euros was paid for agri-environmental commitments. However, most of the funding allocated to Axis 2 did not target biodiversity conservation directly and included measures on financial support for certification schemes, such as GlobalCAP, in intensive citrus crops and potatoes. A measure specifically targeting HNVf conservation was included in the 2014-2020 RDP after the mid-term evaluation of 2017, just after the adoption of the current HNVf map. This measure targets cereals and tree crops: olive, carob and hazelnut.
Pressure on EU habitats and ecosystems is still prevalent. Only 17% of EU habitats and 11% of ecosystems are considered to be in a favourable state, nutrient surpluses persist in some water bodies (despite progress in others) and 45% of EU soils suffer from problems of quality (Pe'er et al. 2014). These challenges need to be addressed and the positive environmental contribution of farming and forestry should be strengthened. The Lawton Review (Lawton et al. 2010) advocates management of the countryside at landscape level to ease pressures on Protected Areas. HNVf areas have received increasing attention at the European level due to their importance for biodiversity conservation and their potential for establishing future connections amongst different protected areas (Strohbach et al. 2015). Even under the conservative storyline i.e. 25% threshold, in Cyprus more than 22.5% of the island's territory falls into one of the three HNVfs' categories, bearing in mind that this study did not map the whole island. In addition, in many areas, there is an overlap in HNVf types (Table 2) pointing to the difficulty in characterising these areas in Cyprus under the existing HNVf definitions. This highlights the importance for managing HNVfs for promoting nature conservation and safeguarding species and processes at regional and national level (Morelli 2018), in a complementary way to the existing Natura2000 network. Examples are already emerging in Europe where this is already taking place (Klimek et al. 2014;Sigura et al. 2010. There is a two-fold importance in this, since landscape level management has also been advocated as part of a wider adaptation policy to tackle climate change effects on biodiversity (Vos et al. 2008). These areas have the potential to facilitate spatial planning for nature conservation since establishing new protected areas will be more likely met with resistance due to competition with other land-uses. This is even more important in an island context where land availability is limited. Improving the conservation value of agricultural lands has been proposed as a complementary strategy (Troupin and Carmel 2014). However Kleijn et al. (2011) contest this and suggest that it is unknown how the extensive European agri-environmental budget for conservation on farmland contributes to the policy objectives to halt biodiversity decline.
Currently, the terrestrial part of the designated Natura2000 sites on Cyprus (CYG-CA) cover 752.6 km 2 and, although the percentage of habitats in a favourable state within the network is reported to be the highest in the EU (EEA 2016), a recent study on the effectiveness of the protected areas on the island suggests that this is far from adequate for conserving biodiversity (Christodoulou et al. 2018). These gaps more often than not relate to the site specific conservation measures implemented in Natura2000 sites without so much consideration of landscape level processes such as habitat connectivity or fragmentation (Opdam et al. 2006). Since establishing new protected areas to enhance the networks connectivity is unlikely, protecting and maintaining these HNVf areas can facilitate spatial planning and help address these gaps in nature conservation by designing appropriate ecological networks. As demonstrated by the results herein, HNVfs provide links amongst all Natura2000 sites which may improve the connectivity of the network and therefore its effectiveness for biodiversity conservation, although quantifying the degree of connectivity was beyond the scope of this study.
In particular, HNVfs can provide buffering from externalities to the lowland Nat-ura2000 sites which are under more pressure compared to their upland counterparts, due to their proximity to intensive agricultural and tourist activities. This is of utmost important for sites Potamos Panagias Stazousas (CY6000007) -Kosiis Pallourokampou (CY6000009) which host species such as Sylvia melanothorax, Oenanthe cypriaca, as well as Falco tinunculus and the long-legged buzzard (Buteo rufinus). The next step would be to identify and design ecological networks amongst Natura2000 sites with their constituent elements (buffers, corridors, stepping stones). A target site could be the Alykes Larnakas (CY6000002) where, physically, the site is split into two due to the construction of the international airport. Starting with the HNVfs' mapping, the challenge in Cyprus is to highlight synergies between nature conservation and agricultural production and to raise awareness about their importance and the financial opportunities they provide.
The high percentage of HNVfs in Cyprus may be good news for conservationists but disquieting for farmers unless they are somehow convinced that HNVfs have real benefit for them. Farmers, but also other stakeholder groups, as demonstrated during the workshops held for this case study, are apprehensive when it comes to lines on the map. When the discourse is about conservation elements, the assumption is that delineation will lead to designation and strict protection where restrictions to land management will apply. Already the erroneous perception on the concept of Natura2000 proves that management in the private lands within the Natura2000 network is problematic in Cyprus but also elsewhere in Europe (Hiedanpää 2002). A review by Plieninger and Bieling (2013) highlighted the vulnerability of HNVfs to socioeconomic changes and the need for a new dynamic and adaptive strategy for increasing their long term resilience. Such a strategy though can be effective only if supported by reliable HNVf mapping.