Pollination services mapping and economic valuation from insect communities: a case study in the Azores (Terceira Island)

Insect pollinators provide vital ecosystem services through its maintenance of plant biological diversity and its role in food production. Indeed, adequate pollination services can increase the production and quality of fruit and vegetable crops. This service is currently challenged by land use intensification and expanding human population growth. Hence, this study aims: (1) to assess the pollination services in different land uses with different levels of disturbance through GIS mapping technique using insect pollinators abundance and richness as indicators, and (2) estimate the economic value of pollination by insects in agricultural crops. Our study takes place in a small oceanic island, Terceira (Azores, Portugal). Our results showed, remarkably, that not only the pristine vegetation areas, but also the orchards and agricultural areas have relatively high values of pollination services, even though both land uses have opposite disturbance levels. For the economic valuation, we analyzed 24 crops in the island and found that 18 depend on pollinators with one-third of these crops having 65% or 95% dependence on pollinators. The economic contribution of pollinators totals 36.2% of the total mean annual agricultural income of the dependent crops, highlighting the importance of insect pollinators in agricultural production and consequent economic gain productions. Nature Conservation 18: 1–25 (2017) doi: 10.3897/natureconservation.18.11523 http://natureconservation.pensoft.net Copyright Ana Picanço 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
Research at the interface of ecology and economics to characterize, value, and manage ecosystem services (henceforth ES) has supported a paradigm shift in how society thinks about biodiversity, ecosystems and human relationships to them (MEA 2005;TEEB 2010;Garbach et al. 2014). This awareness of the ES started with classical papers of Daily (1997) and Constanza et al. (1997); and in 2005, the Millennium Ecosystem Assessment (MEA) promoted and defined the concept of ES as "the benefits that humans recognize as obtained from ecosystems that support, directly or indirectly, their survival and quality of life". MEA suggests to group ES into four categories: (1) provisioning services, such as food, water, timber, and fiber; (2) regulating services that affect climate, floods, disease, wastes, and water quality; (3) cultural services that provide recreational, aesthetic, and spiritual benefits; and (4) supporting services such as soil formation, photosynthesis, and nutrient cycling (MEA 2005).
The valuation and mapping of ES constitutes a continuous and very complex work for several national governments and organizations, and this process is only currently available for few countries (e.g. Portugal, Pereira et al. 2009;UK, Maresca et al. 2011;France, Watson et al. 2011). ES assessment aims usually to estimate of the marginal values of these services to inform decisions and to evaluate how trade-offs in ES provision will affect human well-being. Therefore, researchers are interested in developing methods for quantifying the provision and value of ES so this information can be incorporated into mapping, planning and decision-making at different scales and in different public and private sectors (see e.g., Losey and Vaughan 2006;Allsopp et al. 2008;Nelson et al. 2009;Tallis and Pollaski 2009;Villa et al. 2009;Maes et al. 2012;Nemec and Raudsepp-Hearne 2013;Nahuelhual et al. 2013;2015).
Pollination together with seed dispersal is considered as one of the key ES, classified by the Common International Classification of Ecosystem Services (CICES) coding system (Haines-Young and Potschin 2013) as a "Regulation & Maintenance ES" with the code 2.3.1.1. Among other studies, Klein et al. (2007), Aizen et al. (2009), ), Calderone (2012 and Giannini et al. (2015) show that pollination services contribute significantly to the agricultural production and subsequently assures 75% of food production worldwide ) (as well as to other flowering plants) by ensuring plant reproduction, fruit set development and dispersion (e.g. Ollerton et al. 2011;Altieri et al., 2015). Notably, the pollination of some vegetable crops (e.g. cabbage and other brassicas, carrots, turnips, lettuce, chicory and onions) increases the quality of the seed production (Gallai and Vassière 2009). In addition, insect pollinators enhance fruit and seed quality (Garibaldi et al 2013;Bartomeus et al 2014;Garratt et al. 2014;Marini et al. 2015;Saeed et al. 2016) and reinforces pest management (Cross et al. 2015) which constitutes an indirect and difficult benefit to measure, but extremely important for the agricultural market. Also, a recent study on pollination by wild insect pollinators has showed their capacity to increase the seed production in 41 agricultural systems globally, regardless of the abundance of honey bees (Garibaldi et al. 2013). Additionally, it was also documented that wild insect pollinators can buffer the impact of climate change on crop production ), most likely due to their high biological diversity that can in turn stabilize ES against habitat disturbances (Cardinale et al. 2012).
Besides these findings, there is also a general consensus that native pollinators abundance and richness are declining throughout the world (Ghazoul 2005;Biesmeijer et al. 2006;Winfree et al. 2009). This global decline has sparked the formation of a global policy framework for pollinators, primarily through the International Pollinator Initiative within the Convention of Biological Diversity (CBD) and several other programs (e.g. Food and Agricultural Organization (FAO) Global Action on Pollination Services for Sustainable Agriculture; Bee Life European Beekeeping Coordination). All of these initiatives emphasize the need to assess and monitor the pollinators in different regions in order to better plan their conservation, restoration and to preserve the ES they supply for humans.
The Intergovernmental Platform on Biodiversity and Ecosystem Services (IPBES) from United Nations Environment Programme (UNEP; Zisenis 2015; Schmeller and Bridgewater 2016) was recently created as a Knowledge-Policy interface (Díaz et al. 2015;Schmeller and Bridgewater 2016). In the fourth plenary of IPBES (IPBES-4) the agenda's item 5 (work programme of the Platform) included the development of works towards the approval of the thematic assessment on pollinators, pollination and food production ("Deliverable 3a" -see http://www.ipbes.net/workprogramme/ pollination). This "Deliverable 3a" highlighted substantial knowledge gaps in different regions on the status and trend of pollinators and pollination, making the global assessment of insect pollinators (henceforth IP) not possible due to lack of data, although regional and national assessments indicated that more than 40 % of insect pollinators are threatened locally (Schmeller and Bridgewater 2016).
These knowledge gaps unveil how the interactions between plants and insects are numerous and complex. So, the understanding of how plant-insect species' interactions affect ecological functions and are affected by land management ) is central to maintain and enhance associated ES. As a vital and increasingly threatened ES, pollination Potts et al. 2010) has become an often-cited example of how the ES are economically valuable . Two additional recent examples of studies about ES pollination in Europe (EU) that complement each other are from Leonhardt et al. (2013) and Schulp et al. (2014), both showing results that provide an overview of ES importance, variation and influence throughout European regions.
In this work, we assess the ES provision and values provided by insect pollinators in the Azores archipelagic region (Portugal) where few studies on ES assessment (e.g. Cruz et al. 2011;Mendonça et al. 2013;Vergílio et al 2016) or related to pollination and seed dispersal services have been undertaken (e.g. Pereira 2008;Heleno et al. 2009;Olesen et al. 2002Olesen et al. , 2012. We use a database on the spatial distribution of insect pollination in Terceira Island (Azores) recently collected (Picanço et al. 2017) to provide the first insight of the bees and other IP contribution to the pollination services and for assessing pollination-related ES in a small oceanic island. With this purpose, we applied two types of methodological approaches: (1) mapping pollination services with geographic information systems (GIS; e.g. Nemec and Raudsepp 2013) using bees and other IP abundance and richness numerical values as indicators; and (2) economic valuation -through the production function approach -by using crops production estimates and crops dependence ratio ; Gallai and Vassière 2009;. Our goals were to determine: (I) the spatial variations of the pollination services; (II) whether the variations of the pollination services were influenced by the different land-uses and/or level of disturbance; (III) the number of crops for which production has a certain level of dependence on IP (or vulnerability ratio); and (IV) estimation of the island's IP economic value.

Study area and sampling sites
Terceira Island, with an area of approximately 402 km 2 (length=29 km and width =17 km) is a small island of the central group of the Azores archipelago (Portugal), located in the North Atlantic Ocean (38°37'N, 38°48'N, 27°02'W, 27°23'W). Like the other islands of the archipelago, Terceira is of volcanic origin and the third oldest island after Santa Maria and São Miguel, with an age of about 3,52 million years (Forjaz et al. 2004). The island is formed by four main volcanic complexes namely Cinco Picos, Guilherme Moniz, Pico Alto and Serra de Santa Bárbara, the latter corresponding the highest point of the island (1023 meters).
Terceira climate is temperate oceanic, characterized by both high levels of relative atmospheric humidity and low temperature fluctuations throughout the year. Particularly, winter and autumn are marked by heavy and regular precipitations often associated with strong winds. The average annual precipitation exceed 3400 mm in "Serra de Santa Bárbara" summit, and reaches almost 1000 mm per year in all the island. The average annual temperature varies between 9°C in "Serra de Santa Bárbara", to 17°C on the coast. Minimum temperature in the winter varies between 4°and 12°C while maximum temperature in the summer varies between 14°and 26°C .
The insects (Suppl. material 1: Table S1) were observed from five relevant habitat types, corresponding to an increasing gradient of disturbance, namely natural forests (NatFor) mainly characterized by Juniperus-Ilex montane forests and Juniperus woodlands, naturalized vegetation areas (NatVeg) composed by Pittosporum spp. and Rubus spp., exotic forests (ExoFor) with Criptomeria japonica and Eucalyptus sp., semi-natural pastures (SemiPast) with Lotus sp., Holcus sp., Rumex sp. and intensively managed pastures (IntPast) with Lolium sp. and Trifolium spp.. These habitat types were previously selected according to landscape disturbance index from  see supporting information), with the aim to assess the impact of land-use change on flower-visiting insect species community structure in Terceira Island (for further details see Picanço et al. 2017). In each habitat type, 10 sites were selected. In each site, 10 meters' linear transects with 1 meter width were set up (Pollard and Yates 1993), making a total of 50 transects located across the entire island ( Fig. 1) (for further details on the sampling protocol see Picanço et al. 2017).Ecosystem service mapping Digital Elevation Models (DEM) are a numerical representation of topography, made up of squared equal-sized grid cells (pixels) with an elevation value associated to each pixel. DEM constitute the most widely used data structure to store and analyze topographic information in GIS (Rishikeshan et al. 2014). The pollination service mapping was performed with the ArcGIS10© software, by applying the "Topo to Raster" interpolation technique, which was designed for the creation of hydrologically correct DEMs. This method uses an iterative finite difference interpolation technique. It is optimized to have the computational efficiency of local interpolation methods, such as inverse distance weighted (IDW) interpolation, without losing the surface continuity of global interpolation methods, such as Kriging and Spline. It is essentially a discretized thin plate spline technique for which the roughness penalty has been modified to allow the fitted DEM to follow abrupt changes in terrain. Furthermore, the quantity of input data can be up to an order of magnitude less than that normally required to adequately describe a surface with digitized contours, further minimizing the expense of obtaining reliable DEMs (Wahba, 1990, Hutchinson 1988, 1993ESRI 2016).In this work, DEM were generated using respectively as elevation data the bees and insect pollinators' abundance and richness quantitative information collected from field surveys, of the 10 transects of each habitat type (or land use). We've chosen to separate the bees and total insect pollinators data, because many studies about pollination services are more related to bees than to the insect pollinators in general, and also, to analyze if there would be differences between the DEM of the possible pollination services contribution from these two groups of data. This latter also applies relating to the abundance (i.e. number of individuals) and richness (i.e. number of species) information on both groups (Suppl. material 1: Table S1). In this way, by applying all the fieldwork data, we intend to be more accurate as possible while developing DEM that deliver information on pollination services.
To complement this spatial analysis, we applied the formerly mentioned index of landscape disturbance metric based on the attributes of the landscape matrix ). This index, ranging from 0 to 100, corresponds to a local index of disturbance by taking into account the level of disturbance in the surrounding areas. Values of the disturbance index (D) was obtained by ranking the different land uses attributing a value of "local disturbance" (L) on a land use map of 100 × 100 m resolution built from aerial photography and fieldwork, and for each 100 × 100 m cell the D was calculated (see  and Suppl. material 1: Fig. S1).
For each analysis, we overlaid the respective pollination services' interpolation maps delivered by the fieldwork data on bees and other insect pollinators from Picanço et al. (2017) with the land use and the disturbance index D. We've created thresholds to analyze disturbance index D influence on the amount and diversity of bees and other insect pollinators and mapped these categories in eight classes for bees' abundance (N) and richness (S); and in 12 classes for insect pollinators' abundance (N) and richness (S). The created thresholds values for the different classes are specified in Table 1. The numbers of classes established follow the minimum and maximum abundance and richness values (Suppl. material 1: Table S1) obtained by Picanço et al. (2017) for the different habitat types -natural forest, naturalized vegetation areas, exotic forest, semi-natural pasture and intensively managed pasture. The exceptions are urban, agriculture and orchard areas due to unavailable technical resources. Bees (Hymenoptera) a very important functional group, are constituted by the following most abundant species Apis mellifera, Bombus ruderatus and Lasioglossum spp., while the other wild insect pollinators groups are constituted by Coleoptera, Diptera and Lepidoptera, being the most abundant species Anaspis proteus, Meligethes aeneus, Stomorhina lunata, Rhinia apicalis, Episyrphus balteatus, Eristalis tenax, Hipparchia azorina azorina and Pieris brassicae (for further information related to the species list see Suppl. material 1: Table S1 in supporting information).
The disturbance level was organized in four classes, including a first one with very low disturbance level typical of high altitude native forests (D<20), two intermediate classes and finally a class with high levels of disturbance (D>40). The number of individuals of bees was divided in two classes in a logarithm scale (less than ten and more than ten individuals). The number of species of bees was divided in two classes with one species and two or more species. For insect pollinator abundance and richness three classes were prepared: for abundance, we created one for the rarest species, one for intermediate and one for the most abundant; for species richness we divided the classes arbitrarily in less than 10 species, 10 to 15 and more than 15 (see Table 1). These created classes were evaluated through a quantitative analysis of the area covered by each class in Terceira Island.Economic valuation Terceira Island's main economic activity is agriculture, with the production of dairy products and raising livestock. Many small farmers practice subsistence agriculture or produce in small quantities to cooperatives. The island consumer is relatively similar to the southern Europe consumers, when comparing the GDP per capita of Azores region and Portugal to other countries of Europe (Suppl. material 1: Tables S2, S3), with Azorean economy comprising a conventional interval of prices elasticities -1.2 and -0.8, as in Gallai and Vassière (2009). FRUTER/Frutercoop is the "Association of Producers of Fruit, Vegetables and Flowers' in Terceira Island". Using their data from 2011 to 2015, we calculated the mean annual productions of 24 common fruits and vegetables in this island. Five-year means were used instead of the latest yearly production figures, in order to smooth out annual variations in crop output.
We estimated the value of pollination gain in agricultural crops and its respective vulnerability by using the crop production amount (Kasina et al. 2009), market and producer prices for each crop. This method was adapted to a regional rating scale, according to the methodology of FAO  previously developed by . The data on crops were derived from multiple sources: Klein et al. (2007; only for crops grown in Terceira Island), FAO  ,  Table 1. Distribution of disturbance index (D) for bees' and insect pollinators' abundance (N) and richness (S) per classes. D  N  S  IP class  D  N  S  1  D<20  >10  >2  1  D<20  >73  >15  2  D<20  <10  <2  2  D<20  25<S<73  10<S<15  3  20<D<30  >10  >2  3  D<20  <25  <10  4  20<D<30  <10  <2  4  20<D<30  >73  >15  5  30<D<40  >10  >2  5  20<D<30 25<S<73  10<S<15  6  30<D<40  <10  <2  6 Klein et al. (2007), and posteriorly adapted by , into the following classes: essential, great, modest, little, increase seed production, increase breeding and no increase. We also corresponded the dependence ratio (DR) to these classes according to : essential, DR = 0.95 (meaning that the value of pollination-driven yield lies between 90 and 100%); great, DR = 0.65 (40-90% of yield is dependent on pollination); modest, DR = 0.25 (10-40% of yield is dependent on pollination) and little, DR = 0.05 (0-10% of yield is dependent on pollination). We multiplied this ratio by the economic value of the mean annual crop production to obtain the pollination services' economic value . The production value was obtained through the market prices and producer prices provided by the regional authority -"Serviços de Desenvolvimento Agrário da Terceira" (2016). For the current assessment we did not consider currency values, regional or seasonal variations in the crop labour costs and food prices.

Ecosystem service mapping
By analyzing together both the land use map of Terceira Island (Fig. 1) and the four pollination services' interpolation maps (Fig. 2) we can observe that: (i) bees abundance (N) comprised by some abundant species like Bombus ruderatus and Lasioglossum morio (Suppl. material 1: Table S1) presented higher density values around the northwest,  east, south-eastern coast and also at north, near the centre of Terceira island, matching especially with the current areas occupied by orchards and agriculture; (ii) bees richness (S) high density values also correspond mostly to orchards and agricultural areas, namely in the north, along the west to the southern coast and in-between the centre and the eastern side of Terceira island; (iii) insect pollinators (IP) abundance (N) with the most abundant species being Anaspis proteus, Stomorhina lunata, Eupeodes corollae, Sepsis neocynipsea and Pieris brassicae azorensis (Suppl. material 1: Table S1) presented higher density values around the north-western coast till near the center, and also in the eastern and central parts of the island, corresponding these higher density spots to the main Terceira island's biodiversity hotspots (pristine vegetation forests): "Serra de Santa Bárbara" and "Pico Alto" (that are both classified as protected areas). In the south-eastern coast of the island some orchards and agricultural areas also presented high IP abundance; finally, (iv) insect pollinators (IP) richness (S) comprised by many hoverfly species (Diptera, Syrphidae) when compared to the other insect pollinators groups Coleoptera and Lepidoptera (Suppl. material 1: Table S1) followed a very similar spatial pattern to that of IP abundance. Nevertheless, orchards and agricultural areas in the north-western coast also presented high density values of IP richness.
In order to strengthen the previous analysis, we assessed the influence of the disturbance index (D), as calculated by , in the pollination services and also assessed the area covered by bees and IP classes within the island (Tables 2-5).
As a result of overlaying each previous pollination service output with the matching disturbance index D spatial data (see full description of classes in Table 1), we observed that "Class 1" spatial distribution (areas with disturbance index D lower than 20 and high values for both abundance (N) and richness (S) of bees and IP) corresponded in every output to the small areas of pristine vegetation (biodiversity hotspots) at high altitudes and consequent most difficult human access, namely "Serra de Santa Bárbara" and "Pico Alto" protected areas (Fig. 3), which corresponds to the smallest % of island area (between 0.06 -0.56%) occupied (Tables 2-5).
According to the same Fig. 3 and to Table 1, classes 4 and 6 for bees' abundance (N) and richness (S) ( Table 2 and 3), as well as classes 5 and 8 for IP's abundance (N) and richness (S) ( Table 4 and 5), respectively, are the predominant spatial patterns around class 1's areas.
Moreover, both bees-related maps (abundance -N and richness -S) in Fig. 3, Tables 2 and 3 show that the whole island is predominantly covered by highly disturbed areas (disturbance index D higher than 40) that seriously affect these pollination services, resulting in low abundance (N) and richness (S) for bees (classes 7 and 8). In fact, for the bees' abundance (N), class 8 covers the north to north-eastern coast, passing through the centre until the western coast. Class 7 is predominant from east to the southwestern coast. Regarding the bees richness (S), the class 8 occupies the centre and the area from north to the south-eastern coast, as class 7 covers the areas from northwest to south and the territory between the centre and the eastern coast of the island. Both classes 7 and 8 mostly occur in orchards/agricultural areas, and in IntPast land use, respectively (see Fig.1 and Fig. 3).    In the case of IP-related maps (Fig.3, Tables 4 and 5), the spatial pattern of disturbance versus pollination services is quite similar to that of bees' pollination services. Highly disturbed Class 11's areas (see Table 1) are predominant in the whole island for both IP's abundance (N) and richness (S), occupying around 39% and 20% respectively (Tables 4 and 5). In the case of IP abundance (N), this class covers relevant areas in the north-western, eastern, south-eastern and southern territories of Terceira Island. For IP richness (S), class 11 covers large areas in the west, south and east of Terceira Island (Fig. 3). Most disturbed areas with lower IP-related services (abundance and richness) performance mostly occur in orchards, agricultural areas and other land uses strongly affected by human activity (Fig. 1).

Economic valuation
According to the data provided by Frutercoop for the period between 2011 and 2015, the total value of production for the 24 referred crops in Table 6 represented an amount of €874,925.51, from which only 29% of the production is from crops with known pollinator dependence ratio (see Tables 6 and 8).
In terms of welfare, an assessment of the social cost to Terceira Island consumers resulting from pollinator decline estimated that the consumer surplus (economic measure of consumers benefit) loss was from €156K to €231K, which reflects the impact on the price of the crop on the market, based upon average price elasticities of −1.2 to −0.8, respectively (Table 7). When considering these values, we must also take into account that the production from Frutercoop represent approximately 54% of the entire island's production.
Among the 18 crops relatively dependent to IP, the greatest economic value generated by the IP was originated by the class "little" or DR = 0.05, with 46.9% (€119,833), as well as the one originated by the class "great" or DR = 0.065, with 29.5% (€75,465) (Tables 6 and 8). In each class "little" and "great", the most representative crop productions were respectively "tomatoes" and "apples" (Table 6).
On average, in recent years (2011)(2012)(2013)(2014)(2015), IP contributed to pollination service in crop production with about €91,957 (total economic value of IP, EVIP), representing 10.5% crops ratio of vulnerability (VR) ( Table 7), according to Frutercoop dataset. Extrapolating the IP contribution estimation from Frutercoop production data to the entire island, we can consider the IP value to be of approximately €170,291. This value represents about 36.2% of VR from the mean annual agricultural income (€469,867) resulting from the dependent crops.

Discussion and conclusions
Under the same thematic as "Deliverable 3a" from IPBES (Schmeller and Bridgewater 2016) -i.e., assessment of the contribution of insect pollinators to the pollination and Table 6. Array of crops used directly for human food following FAOSTAT and listed by common names of crops.  food production -our findings highlight the great importance of insect pollinators on a small oceanic island economy. Our results are relevant since they are based both in field and economic data with the aim of providing quantitative information as in Leonhardt et al. (2013) and Schulp et al. (2014), but by using a completely different approach to evaluate insect pollinators distribution in comparison to Londorsf et al.

Crop common name
and Polce et al. (2013), which have used other biological indicators and modeling techniques. Concerning the field-based mapping of pollination-related ES, similar spatial patterns were revealed for both bees and overall insect pollinators (IP): (i) high values of abundance (N) and/or species richness (S) are directly associated to the pristine native forest areas with lower disturbance (D), on one side with low percentage of island area covered (Tables 2-5); and (ii) on the other side these same high values of pollination services are also observed in orchards and agricultural areas with high level of disturbance (D) covering large island areas (Tables 2-5). These results show that Azorean native pollinators (e.g. Pieris brassicae azorensis, Anaspis proteus, Lasioglossum spp., Eupeodes corollae, Stomorhina lunata -see Suppl. material 1: Table S1) are provid-ing key pollination services not only in native habitats for which they are originally adapted, but also in low altitude agro-ecosystems in which they expended their range. This finding call for the need of a whole island integrated management strategy for pollinators in Terceira in order to decrease the 32.6% VR of crops production. However, intensive managed pastures, the most dominant land use in the island with highest disturbance index D (see classes 8 and 12 in Tables 2-5), showed low abundance (N) and/or richness (S) classes for both bees and IP (Fig. 3), evidencing therefore a low performance of pollination services, as observed in previous studies (e.g. Batary 2010;Sjödin 2007). Indeed, this land-use, subject to frequent and intense grazing events does not foster the occurrence of abundant pollinator populations. Based on the results obtained for low altitude agricultural areas, the disturbance index D variable, in contrast to other studies (e.g. Boieiro et al. 2013;Cardoso et al. , 2014Florencio et al., 2013), do not fully and adequately explain the spatial abundance of native pollinators in this island. Unmeasured variables associated to current and past land uses that reflect specific agro-ecosystems management regimes in Terceira Island may have driven the current spatial heterogeneity of the pollinators' abundance and diversity. The numerous resources available for pollinators at low altitude (e.g. private gardens, abandoned orchards) together to a low input of pesticides in abandoned orchards are possibly fostering an ideal situation for the spread of native pollinators across the landscape (see also Picanço et al. 2017).
This study also highlights the fact that about one-third of Terceira Island crops have an essential or great dependence on pollinators, therefore complementing the above information on high values of insect pollinator abundance and richness in low altitude agro-ecosystems. The economic contribution of pollinators totalizes 36.2% (€170K) of the mean total annual agricultural income of the dependent crops (€469K). This EVIP percentage represents also the VR of agricultural production. Moreover, the consumer surplus loss was estimated between €156K and €231K based upon average price elasticities of −1.2 to −0.8 respectively. This interval of prices on the consumer surplus loss represents the difference between what island consumer are willing or able to pay for the ES relatively to its market price, in case of pollination services loss. These values referred to Frutercoop production only represents 54% of the island's total crop productions (Tables 6 and 7). However, the presented estimates are underestimated values, since not all agricultural production is officially declared (family production, production in backyards, urban gardens, etc.).
Our study also indicates the high socioeconomic relevance of pollination-related ES in a small oceanic islands' context. Nevertheless, bio-economics based valuation studies have been inherently and generally unable to provide thorough and consistent results, due to frequent changes in currency values, labor costs and food prices. This type of approach has also failed to consider and propose realistic and cost-effective mitigation efforts that might reduce the impact of a pollination crisis. In general, the costs are still being strongly dependent on the local agro-ecological setting, namely the crops phenology, the local insect populations, and the existing ecological relationships between farmland and surrounding natural or semi-natural areas.
Some crops, despite their modest or little dependence, showed very high values of mean annual production and, therefore, even in these cases, the contribution of pollinators is significant (Gallai and Vassière 2009; Tables 6 and 7). Moreover, there is no available information on pollinator dependence for some relevant crops, showing the urgent need to address this issue through basic research on reproductive biology and pollination ecology.
As a result, these pollinator-dependent crops are crucial for maintaining the agricultural food balance of the increasing population-growth of Terceira Island's consumers. Meanwhile, at the world scale, IP are becoming increasingly more vulnerable to (i) land-use intensification (Power et al. 2012); (ii) use of pesticides (Kevan 1999;Suchail et al. 2001;Dos Santos et al. 2016;Geslin et al. 2016 (Gill et al. 2016;Ferreira et al 2016); and (viii) the interactions of these ecological stressors (Potts et al. 2010;Vanbergen 2013). Nevertheless, it seems that intensive pastures aside, IP populations in Terceira Island are abundant and diverse in several agro-ecosystems (Fig. 3), and performing adequate pollination services to crops.
With the expected need for an increased production of vegetables and fruit in Terceira Island in the coming years, integrated mitigation measures (e.g. biological pest control, wild flowering plants production areas, promotion of organic farming), as well as (cost-) effective, innovative and attractive (for farmers) agri-environmental schemes are required in order to adequately promote pollination services and to compensate for some eventual crops' failing production (e.g. Wilson and Hart 2001; Power et al 2012; Andersson et al. 2014). It appears to be increasingly consensual that organic farming regimes benefit biodiversity, zoophilous wildflowers and IP abundance on a local scale (Gabriel and Tscharntke 2007). As such, if strategically and effectively promoted and applied, this management practice may have the potential to benefit crop pollination and to increase IP abundance across the whole island. This needs to be taken into account for the sustainable long-term management and conservation of pollinator communities and insect-pollinated plants in Terceira Island (e.g. Power et al. 2012).
Agri-environmental schemes aiming to foster and to pay/compensate farmers for a more sustainable management of low-intensity pasture systems and to implement integrated farm management and organic agriculture practices should be especially encouraged in the north-western, eastern and south-eastern agro-ecosystem areas of Terceira Island.
Finally, this broad, straightforward and cost-effective methodological approach may be able to be applied in further small oceanic islands with the aim of improving the capacity of effectively assessing and monitoring pollination-related ecosystem services, in order to improve the existing decision support systems for land use planning/ management policies, especially those related to agriculture and nature conservation.

Introduction
Species distribution models (SDMs) use environmental data from known locations of a species to predict places where that species could potentially occur within landscapes or regions (Booth et al. 2014). SDMs have been used to identify critical habitats for species with greatly reduced distributions (Hamilton et al. 2015;Jetz and Freckleton 2015;Manthey et al. 2015), provide potential translocation sites for species based on known habitat requirements (Adhikari et al. 2012) and predict the movement of invasive species across landscapes under different scenarios (Kearney et al. 2008;Elith et al. 2010). Spatially explicit probability of presence, or prediction of occurrence maps, generated using SDM algorithms, have been used to inform conservation planning and habitat management at both coarse and fine scales. Consequently, they can guide or prioritise future survey efforts and aid in assessing the conservation status of target species. SDMs relate known occurrences of a species with various environmental variables and predict a probability that a species will occur in areas where no data on its occurrence is available. Thus, they help to identify potentially suitable habitat Guisan and Thuiller 2005). The accuracy of a SDM depends on such factors as the quality and appropriateness (in regard to sample size and representativeness) of the presence and/or absence data for the target species or community, the expertise of the modeller, the selection of an appropriate modelling tool (or software package), the selection of an appropriate suite of predictive/independent variables, the quality of the variable data used, and an acknowledgement of the strengths and limitations of the SDM (Elith et al. 2010;Guillera-Arroita et al. 2015). However, where there are shortfalls in appropriate data, modellers often compensate by focussing their efforts on the development and evaluation of novel methods to improve the performance of their models, and thus, to better predict the environmental suitability for species in applied studies (Barbosa and Schneck 2015;Guillera-Arroita et al. 2015).
In this study we set up a SDM to determine the potential distribution (PD) of the northern quoll Dasyurus hallucatus population found in the Pilbara biogeographic region of Western Australia (Thackway and Cresswell 1997). This is a unique population of a threatened species for which conservation management to date has been limited by significant knowledge gaps (Cramer et al. 2015). Limited and largely opportunistic sampling has resulted in constricted and potentially biased presence data and this in turn has resulted in problems in determining appropriate predictive or independent variables. Furthermore, the literature associated with these SDMs shows no evidence that the models have been tested for covariance and sample bias as described by Hijmans (2012) and Fithian et al. (2015), nor have differences between these SDMs been evaluated or explained.
Northern quolls are a suitable subject for distribution modelling as they have a strong habitat affiliation with complex rocky areas, often in close association with permanent water (Begg 1981;Braithwaite and Griffiths 1994;Oakwood 1997;Pollock 1999;Schmitt et al. 1989). In the Pilbara, quoll habitat affiliation aligns with mesas or ranges which are often the focus of iron-ore extraction (channel-iron deposits and banded iron stone formations) and granite outcrops which are often quarried for road and rail bed materials (Ramanaidou and Morris 2010).
Once widely distributed from the Pilbara region of Western Australian (WA) across northern Australia to southern Queensland (Figure 1), the mainland distribution of the northern quoll has now contracted to several disjunct populations (Burbidge et al. 2009;Oakwood 2008). This collapse has largely been linked to invasion by the toxic cane toad Rhinella marina (Braithwaite and Griffiths 1994;Doody et al. 2015;How et al. 2009;Oakwood 2004;Woinarski 2010). Other impacts currently causing rapid and severe declines in northern Australia's critical weight range mammal fauna (i.e. terrestrial species within the weight range 35 -4200 g are considered to be particularly vulnerable to introduced predators) are also likely to also be impacting on the carnivorous northern quoll (Burbidge and McKenzie 1989;Cramer et al. 2016). These include altered fire regimes, the grazing impacts of introduced herbivores, climate change and both enhanced mortality and competition for resources (including prey animals) from introduced predators, in particular the feral cat Felis catus (Burbidge et al. 2009;Cook 2010;Woinarski et al. 2015). As a consequence of all these recent declines, the northern quoll is listed as endangered under both the Commonwealth's Environment Protection and Biodiversity Conservation Act 1999 (EPBC Act 1999) and the Western Australian Biodiversity Conservation Act 2016.
The main WA populations of northern quoll occur in two discrete mainland regions, the Kimberley and Pilbara, separated by the arid Great Sandy Desert. Both mitochondrial DNA sequences and nuclear microsatellite loci reveal clear differentiation between Kimberley and Pilbara populations and a greater distinction between these populations than those in the Northern Territory and Queensland (Spencer et al. 2013;Westerman and Woolley 2016). These WA populations also differ from those remaining in Queensland and the Northern Territory in regard to both genetic structure and demographic parameters and represent the last intact populations in Australia that have not experienced major declines subsequent to the introduction of the cane toad and consequently display the highest levels of genetic integrity (How et al. 2009;Spencer et al. 2013;Spencer 2010).
Given that the Pilbara population of the northern quoll is genetically and demographically distinct from all other populations, retains its pre-European genetic diversity, is currently outside of the cane toad's distribution, and has much of its habitat still intact, this population has been assigned a high conservation, research and management priority (Cramer et al. 2015).
The available presence data for this population is clustered around areas of development interest to the mining industry, or where targeted surveys have occurred. This begged the question: was this an example of sample bias or a true representation of northern quoll distribution? Sample bias, where sampling has not been uniform over the project area, e.g. where only easily accessed areas, or known populations have been sampled, has the potential to distort a SDM ). Lacking true absence data for this exercise and being aware of the limited capacity of pseudo-absences to compensate for high levels of sample bias ), we sought to find a method to eliminate or minimise sample bias in our SDM.
Our objective was to construct a high resolution SDM for the Pilbara northern quoll by applying an innovative form of bias compensation to a well proven modelling method, MaxEnt, and testing this SDM with an ensemble model.

Study area
The area modelled for this exercise is the Pilbara biogeographic region (Fig. 1) (Thackway and Cresswell 1997). This selection encompasses all the known occurrences of the unique Pilbara population at the time of the study, and satisfies the requirements and priorities of this project's stakeholders and funding bodies.

Presence data
All presence data, both for northern quoll and surrogate species, was supplied by the West Australian Department of Parks and Wildlife NatureMap database (Naturemap 2015), and comprised of (~2000) records for this species within the Pilbara up to, and inclusive of, 2014. These species records are both targeted and opportunistic, sourced from museum, and Parks and Wildlife fauna surveys as well as compulsory returns from biological consultants and industry. The NatureMap threatened species database is continually being verified and updated by cross-referencing from grey literature, fauna returns and reporting required under the Western Australian Biodiversity Conservation Act 2016.

Variable selection
To obtain optimum efficiency, minimize multicollinearity and prevent overfitting, the suite of variables used should be kept compact (preferably ≤10 in number) and be comprised of those variables best able to define the PD of the target species or community (Beaumont et al. 2005;Elith et al. 2011;Hijmans 2012;Van Gils et al. 2012). To accomplish this, we reviewed the literature on the Northern Quoll in general and the Pilbara population in particular, to identify independent variables likely to be in- fluential or linked to its distribution and therefore suitable for producing a SDM (see Suppl. material 1). All variable data sets were downloaded at, or converted to, a pixel resolution of 30 seconds (~1km 2 ) using the WGS 1984 datum and clipped to the Pilbara bioregion. Map projection (WGS 84) was consistent in all data sets, which were aligned by importing raster data onto a common grid and converting the outputs to ASCII grid file format using ArcGIS 10.4.
To avoid using unnecessary time and resources in identifying an appropriate suite of predictive variables, a two-stage process was adopted. The first stage used a series of statistical tests (described below) to halve the number of potential variables so as to quickly and effectively limit multicollinearity in the bioclimatic variables and to remove those variables for which their contribution to the model was low or counterproductive. The second stage was to use a stepwise elimination process to identify a final suite of predictive variables suitable for use in all further modelling.
Firstly, to reduce multicollinearity between scalar variables, we calculated both the Pearson and Spearman rank correlation coefficient between each pair of variables using the northern quoll presence data. This was done using the pairs function in the 'psych' R package (Revelle 2014). For each pair of highly correlated variables (r > 0.70), we selected only the single variable deemed to be the most relevant for identifying northern quoll presences based on ecological relevance and expert opinion (Phillips and Dudík 2008). Categorical variables were tested for association with each other using Pearson's chi-square test for significance. Similarly, relationships between categorical and scalar variables were tested using linear models with a 0.05 level of significance (Agresti and Kateri 2011).
The final cut was undertaken through a step-wise elimination process using Max-Ent (Phillips et al. 2006) looking at both the contribution of each remaining variable and the consequences of its omission. This was done by building a MaxEnt SDM using all remaining predictive variables against the Pilbara northern quoll presence-only point data set and then iteratively removing the worst performing variable in regard to variable contribution and jack-knife tests. The model was then re-run and if sufficiently robust, the process was repeated. The model at the start of this process was robust, with acceptable Area Under Curve (AUC) and regularised training gain values, so if the revised model was not sufficiently robust the variable was reintroduced to the model and the process repeated with the second worst performing variable deleted. This process was repeated until we had a minimum number of variables capable of producing and SDM with a minimum AUC value of 0.9. This minimum threshold value was set to ensure a statistically very strong model (Elith et al. 2011).

MaxEnt modelling
For our primary modelling tool we selected MaxEnt for its capacity to produce effective SDMs using presence-only data (Booth et al. 2014;Yackulic et al. 2013). MaxEnt, or maximum entropy, modelling is a machine learning modelling tool which seeks to estimate a target probability distribution by finding the probability distribution of maximum entropy, i.e. where variable parameters are closest to homogenous, subject to the limitations of the data used (Guillera-Arroita et al. 2015). This is achieved by applying the predictive or independent variables against training data which are a subset of presences randomly selected and assumed to be representative for the modelled distribution (Radosavljevic and Anderson 2014).
Some limitations have been recognised with MaxEnt, notably a tendency for it to underperform where there is a biased sample, poorly chosen predictive variables or inadequate testing of results (Bystriakova et al. 2012;Elith et al. 2011;Kramer-Schadt et al. 2013). However, where these limitations are addressed, it remains a well-supported modelling tool because it is relatively easy to use and has a capacity to link fine-scale bioclimatic data to species distributions to produce accurate probability-based outputs suitable for informing conservation management actions (Guillera-Arroita et al. 2015; Phillips and Dudík 2008;Syfert et al. 2013;Williams et al. 2012). Consequently, we constructed our MaxEnt SDM with the following criteria: • Withholding a random 30% of presences for testing purposes over 10 bootstrapped repetitions (Barbet-Massin et al. 2012). • Combining all presences within a pixel (~1km 2 ) as a single presence. This resulted in a reduction in the number of presences from 1984 to 324.

Bias compensation
To compensate for our limited presence data we followed Fithian et al. (2015) in constructing a bias layer using substitute species which we applied only to the MaxEnt model. Therefore, we generated a bias grid by substituting records of other non-volant critical weight range (CWR) mammals of 35-4200g (Burbidge and McKenzie 1989), obtained throughout the Pilbara in lieu of northern quoll presence/absence data. These data were used as a presence/absence point data set that would reflect a sampling effort likely to detect northern quoll. The reasoning behind this is: 1) CWR presence records reflect sampling effort; 2) sampling for non-volant CWR mammals would most likely result in northern quoll presence records (e.g., capture, sighting, tracks, scats, or other physical evidence such as remains) if indeed they were present; 3) therefore, presence records for non-volant CWR mammals are suitable for use as pseudo-absence data; and 4) a point density analysis of non-volant CWR mammal presences for the whole of the Pilbara would indicate the degree of bias present in the northern quoll presence records and could therefore be used as a bias grid in the MaxEnt model. To construct this bias grid, presence records for all non-volant CWR mammals (including northern quoll) in the Pilbara were gathered from the Department of Parks and Wildlife Nature Map data base (Department of Parks and Wildlife 2007-) and categorised into northern quoll presences and pseudo-absences. We note that, although this was a separate data set to that of the original presence data set, many presences were replicated. All records were then used to conduct a Point Density Analysis (PDA) using the Point Density function of the Spatial Analyst toolbox in ArcGIS 10.3. This function counted the total number of records for each cell within a 44 cell radius (the default radius). The resulting shapefile was then directly incorporated as a bias grid into the MaxEnt SDM.

Testing the SDM with an ensemble package
As the above SDM was compiled using just one modelling tool and as different algorithms and methodologies can yield very different and often contradictory results, we opted to test the rigour of the preferred (MaxEnt) SDM using ensemble modelling techniques. This involved compiling a suite of different algorithms to construct multiple SDMs for the target species within a single platform and then combining these SDMs to produce a single ensemble, or composite, SDM (Crimmins et al. 2013;Grenouillet et al. 2011). This approach enabled us to compare the MaxEnt SDM with individual and ensemble model outputs and differences between modelling algorithms to be compared.
The ensemble modelling was undertaken using the biomod2 package in the R platform (Thuiller et al. 2013). This package allowed the use of the same variable, presence and pseudo-absence data used to develop the preferred MaxEnt SDM.
We selected the five best performing modelling algorithms for our ensemble model. These were Generalised Linear Model (GLM), Generalised Additive Model (GAM), Generalised Boosted Model (GBM), Flexible Discriminant Analysis (FDA) and Multiple Adaptive Regression Splines (MARS). In running these a random 30% of presences would be used to calibrate the model and 70% of presences could be withheld for testing. This process was then repeated 10 times to add rigour to the results. Unlike Max-Ent, the biomod2 package does not provide an option to use a bias layer to compensate for sample bias. Therefore we applied the surrogate presence data set, from which we constructed the bias layer, as a substitute for true absences in our model imputs.
All outputs of all algorithms were evaluated with a True Skill Statistic (TSS), and Receiver Operator Characteristic (ROC -a test comparable with the MaxEnt's AUC statistic) and combined. A weighting was given to each algorithm based on ROC performance and all model outputs were combined to produce a weighted mean SDM which we used as our biomod2 output.
Comparisons between the MaxEnt and biomod2 SDM were again made using the pairs function in the psych R package. This was done by compiling a point data set of 10,000 random points across the study area. This point data set was then used to extract values from both SDMS and the two resulting data sets compared through the pairs analysis.
The individual modelling packages used, their results and the results of the ensemble modelling process are given in Suppl. material 2.

Variable selection
From the broad suite of variables tested (Suppl. material 1) we derrived a final suite of seven variables with acceptable levels of covariance for use in all models ( Figure 2). Spearman rank correlation coefficients are well under the 0.7 threshold. Although Pearson correlation coefficients were also under this value the Spearman value is used as the most appropriate given that not all distributions appear normal. The MaxEnt analysis provided the contribution and importance values for each variable (Table 1).

MaxEnt SDM
The bias file (Figure 3) demonstrates that sampling has not been uniform with most sampling being undertaken in the north east of the region in areas subject to relatively heavy mining activity. On further examination more instensive sampling along infrastructure corridors, e.g railways, major roads and powerline routes, became obvious.
The MaxEnt SDM (Figure 4a) appears to be a robust model with a high average AUC value of 89.5 (the full MaxEnt readout including AUC plots, variable responses,  sensitivity, threshold diagnostic data is given in Suppl. material 2). This model identifies a high probability of occurrence for many areas already known to be northern quoll habitat such as the rocky habitats on the western edge of the Hamersley Ranges, the rugged Chichester Ranges and the granite outcrops of the Abydos Plains (a map of these areas is given in Suppl. material 1). However, it projects beyond known presences to predict a low probability of occurrence in the Fortescue catchment, the sandy coastal regions of the Pilbara and in the southern areas of the Hamersley Ranges, and to predict a high potential for northern quoll habitat in many areas where this species has not been previously identified, particularly in the central west and far eastern parts of the region. When compared to the MaxEnt SDM, constructed without the use of the bias layer (Figure 4b), it appears that the use of the bias layer extends predicted habitat further beyond those highly sampled areas where northern quoll are frequently encountered identifying further areas, where this species has not been, or has rarely been, previously encountered.

Comparison with the biomod 2 ensemble model
The full outputs for the biomod2 modelling process are given in Suppl. material 3. This is a robust model with ROC values for individual SDMs varying between 0.88 and 0.97. To facillitate a visual comparison, the ensemble model has been projected at the same extent, symbology and resolution as the MaxEnt model ( Figure 5). A comparison between this (emsemble) SDM and the MaxEnt model ( Figure 4) shows  a strong general similarity between both models in terms of the areas selected as likely habitat, with the main differences being variations in habitat value of the areas selected rather than which areas were actually selected. The pairs analysis ( Figure 6) shows a strong similarity between these SDMs with an overall Spearman correlation coefficient of 0.86. Histograms show very similar distributions across the full range of both models and like the MaxEnt SDM the biomod2 ensemble model also highlights similar, but little sampled, areas as potential northern quoll habitat.

Discussion
The independent variables that we used were suitably diverse in nature with an acceptable level of covariance within the variables used in the final selection. The use of the jack-knife test to determine the final suite showed that the removal of any of the variables selected in the final suite would have compromised what was a strong model. The literature informs us that the number of variables used was appropriate for a modelling exercise of this nature (Phillips and Dudík 2008;Stockwell and Peterson 2002). Furthermore, the fact that all runs of all models applied these variables differently and, in most instances, allotted acceptable levels of significance to all or most of these variables (while consistently producing high test values), supports this finding. It should  be acknowledged that the bias grid was derived from pseudo-absence data for which suitability was determined through a series of assumptions. Consideration of sampling bias was necessary because of a potential shortfall in suitable presence/absence data and the fact that most of the targeted sampling for the northern quoll has been undertaken over a relatively limited area and timeframe. Specifically, many surveys have been conducted in mining areas and associated infrastructure (either historical, current or proposed) which tend to confine the survey effort to particular types of geomorphology. Surveys of non-volant mammals conducted in the region tend to be for environmental impact assessments associated with mining or part of comprehensive regional surveys (e.g. Bamford Consulting Ecologists (2013); Biota Environmental Sciences (2012); Eco Logical (2012)) or part of comprehensive regional surveys ) and hence absence of northern quoll in these surveys could be reasonably assumed to be true absences. As demonstrated (Table 1, Figure 2, Suppl. materials 2 and 3) northern quolls were found to conform strongly to ecological habitat associated with the vegetation, and slope, topography of the rocky areas of the Pilbara bioregion. Primary areas of northern quoll occupation in the Pilbara included the western edge of the Hamersley Ranges, in the granite outcrops of the Abydos Plain, and in the more rugged areas of the Chichester Ranges. Our models identified a low likelihood of occurrence in the Fortescue River floodplain and its upper catchment, the sandy coastal regions of the Pilbara and in the central and southern parts of the Hamersley Ranges. We also identified several areas with a high probablity of presence which have not, as yet, been adequately sampled and which have not been identified as habitat in previous modelling efforts. This is particuarly so in the eastern Pilbara IBRA region. Here, areas are given a high probability of northern quoll although few to no records of the species exist. This area and others identified as of high probability would be logical areas for further survey effort (as well as providing a ready field-based validation of the model), and demonstrate the utility value of SDMs.
The conservation of a threatened species requires good information about population locations and ecological requirements within its geographic range, particularly when threatening processes are ongoing. Applying appropriately selected surrogate species data to both the MaxEnt and biomod2 software packages has enabled us to develop SDMs which identify remarkably similar core areas of likely quoll habitat, as well as less optimal habitat that may only be occupied in favourable conditions. These apparently less favourable habitats may however be of high conservation value as current information suggests that all Pilbara northern quoll populations are genetically linked, and high level of dispersal occurs between geographically distant populations (Spencer et al 2013, Woolley 2015. Consequently, smaller populations of northern quolls in less-preferred habitat may be important in maintaining gene-flow throughout the region. The SDM pairs analysis ( Figure 6) shows a very strong overall correlation between the preferred MaxEnt SDM and the ensemble model with minor differences in predicted habitat being more about conflicts of scale rather than the actual areas nominated. The correlation tests also indicate a strong similarity between the preferred MaxEnt SDM and the ensemble model overall, and a very strong similarity between these models in identifying areas of high habitat value, i.e. when comparing the highest quartiles. This supports our earlier observation that the difference between the two models is largely one of resolution rather than different predictions.
In comparison with previous SDMs on this northern quoll population (Biologic 2012;CliMAS 2014;Eco Logical 2012) our SDM (particularly the MaxEnt model) picks up many of the same areas identified as having a high habitat value, but with an apparently greater level of definition. This is particuarly so in those areas which have not been heavily and specifically sampled for northern quoll. In short, there is a tendency for these models not to project far beyond those areas where northern quoll are known to exist. As such these models tend to tell us little more than what we already know, thereby limiting their value for conservation planning. This could also be said for the trial MaxEnt and biomod2 SDMs constructed without the use of surrogate species, either in the form of a bias file or as pseudo-absences. In both cases, areas of high quality habtat appeared to be much more limited to areas already known to be habitat and less likely to projet into areas where northern quoll records are either absent or less frequent.
In this study we have demonstrated a methodology capable of addressing three of the more common problems associated with SDMs, specifically how to: 1) address bias in a high resolution SDM over a large and diverse landscape with a limited, and potentially biased, presence only data set; 2) selecting an appropriate suite of predictive variables for the construction of such a model; and 3) establish a means by which the suitability and outputs of a modelling tool can be verified. We have developed an innovative approach to constructing an SDM by pre-emptively identifying problems likely to arise due to data limitations and addressing these issues by reviewing the options available and selecting a combination of responses which minimise bias effects and meet the needs and constraints of the modeller.
We note that a true comparison between a model with randomly selected psuedoabsences and a bias-corrected model using a bias layer created from surrogate presences, remains the preferred way to demonstrate the usefulness of this form of bias compensation. However, in the absence of broad-scale sampling and ground truthing to test truth versus prediction (as is proposed) the demonstrated approach remains the most feasible form of bias compensation under the given circumstances.
Being aware of the resource and skill limitations preventing many conservation managers from constructing SDMs, this methodology was deliberately selected to meet these limitations. All data used was freely available and all software used was freeware, other than the use of a commonly used GIS package that could also be substituted with freeware. Skill levels were limited to what the authors considered average for an ecological project team, i.e. no modelling, statistical, GIS or programming specialists were required for this study.

Conclusion
In this exercise we produced a SDM that predicted areas where new populations and sub-populations of the northern quoll might be found outside of areas currently known to be habitat. By identifying areas of high habitat value, this SDM facilitates the identification and conservation of high priority habitat areas, potential translocation sites and potential movement corridors. As a consequence of having identified a sound suite of predictive variables, our understanding of the habitat requirements of the Pilbara population of the northern quoll has been increased. Finally, by comparing the distributions identified through this exercise with proposed mining and infrastructure projects in the environmental impact assessment process, this SDM can be used to minimise impacts on this unique and important northern quoll population.
The ensemble modelling process validates our choice of the MaxEnt model with bias file as the preferred SDM. However, this is a desktop exercise derived from a relatively small and uniform sample. This model should be validated and refined through on ground sampling and research. Our study exemplifies a preferred practice in the use of SDMs by corroborating our findings with a control, in this case one derived through an ensemble modelling approach. We commend this practice to modellers and caution against outcomes derived from only a single modelling approach. Our study has revealed the most comprehensive known refined distribution map for the endangered Northern Quoll. It has confirmed their reliance on rocky upland habitats and their limited distribution within the Pilbara region. These outcomes will be of great importance to land managers when considering the impacts of planned developments within the region.

Introduction
The increasing amount of infrastructure and road networks in the landscape is a major cause of habitat fragmentation (e.g. Forman andAlexander 1998, Jackson andFahrig 2011). Roads may act as barriers to animal movements, resulting in animals avoiding or not being able to cross roads, but also leading to animals suffering a high mortality when moving between habitats intersected by roads (Forman andAlexander 1998, Trombulak andFrissel 2000). Most studies on barrier effects of roads on animals have so far mainly focused on mammals, birds and amphibians (c.f. Trombulak and Frissel 2000). Considerably fewer studies have focused on insects (Muñoz et al. 2015), despite the fact that insects are the most diverse and species rich animal group on earth. Some studies on the effects of roads on insects exist, indicating that large roads can cause disruption of movements between habitats in butterflies (Askling et al. 2006) and bumblebees (Bhattacharya et al. 2003). There is also some recent evidence that the mortality of insects can be substantially high along large roads (McKenna et al. 2001, Soluk et al. 2011, Skórka et al. 2013, Baxter-Gilbert et al. 2015, Skórka et al. 2015. Barrier effects may also ultimately have the consequence that populations on different sides of roads become genetically isolated (Keller and Largiader 2003). It is likely that barrier effects of roads are most pronounced in species with poor dispersal ability, and studies of butterflies indicate that it is mainly smaller species that avoid crossing roads (Askling et al. 2006) and also suffer from the highest road mortality (Skórka et al. 2013). There is thus a possibility that roads can affect the community composition of insects, particularly with respect to species that are poor dispersers. However, to our knowledge this has not been investigated before.
One additional reason that barrier effects of roads on insects deserve attention is that road verges have been highlighted as important habitat for many species (Way 1977), and may as such be considered for conservation actions. In Sweden and many other parts of Europe, changes in land use have caused dramatic areal declines in grassland and meadow habitats during the last century (Critchley et al. 2004, Cousins 2009), and road verges may today be important surrogates of these grasslands in modified landscapes (Ruiz-Capillas et al. 2013). Road verges can be species rich in vascular plants (e.g. Way 1977) and may therefore function as habitats for pollinating insects (Saarinen et al. 2005, Hopwood 2008). Furthermore, roads can also serve as linear elements in the landscape, thereby facilitating dispersal of insects alongside the roads (Sjödin et al. 2008). However, if roads act as barriers and obstruct dispersal between habitats on different sides of roads for insects inhabiting the road verges, there is a potential conflict between the conservation value of the road verges and a fragmentation effect from the road (Skórka et al. 2013). Aculeata (i.e. bees and wasps) is a species rich group that may provide several ecosystem services (e.g. Harris 1994, Tscharntke et al. 1998. Especially bees, but also to some extent wasps, are flower-visiting insects important for pollination of crops and native plants (Tscharntke et al. 1998). Wasps are also predatory insects, potentially reducing populations of pest insects (Harris 1994). Moreover, these insects have been proposed as important bioindicators for ecological change and habitat quality (Tscharntke et al. 1998). It is known that many bees and wasps utilize road verges for foraging and nesting (Hopwood 2008, Hanley andWilkins 2015), but very little is known about how large roads may act as barriers to movements between habitats and thereby affecting community composition.
The aim of this study was to investigate if large roads may act as barriers on the movements of bees and wasps (Aculeata). Specifically, we compared the species composition of bees and wasps on different sides of a large road, both for the whole community and for only the smallest species. We expected the community composition to differ between the two sides, given that there was a barrier effect. As a complement, we also analyzed general differences in foraging and nesting variables, and the abundance and number of species of bees and wasps, between the two sides of the road.

Field inventory
We performed an inventory of bees and wasps (Aculeata) using yellow pan traps at five sites in Sollentuna municipality, north of Stockholm, during late July-early August 2015. The sites were situated along the highway E4, which runs approximately in a north-south direction through Sollentuna municipality (Figure 1). This segment of E4 highway has a maximum speeding limit of 100 km/h and has a traffic flow of approximately 90 000 vehicles/day. We selected the sites as to be more or less evenly distributed along the highway. However, exact sampling locations had to be determined in the field after considering logistic factors, as large parts of the highway were not always easily accessible. At the sites, the distance from the road to the vegetation was approximately 50 cm, separated by a stripe of gravel.
At each site, we placed four yellow pan traps, two on each side of the road, i.e. 20 traps in total. At four sites the distance between the traps across the road was approximately 50 m and at one site the distance was approximately 200 m. The distance between traps at the same side of the road varied somewhat between the sites, but the distance between the traps on the same side of the road were at each site approximately the same as the distance to the nearest trap on the opposite side of the road (with a square-shaped setup, were each trap was a corner). The vegetation at the sites had previously been mowed during the summer. The exact dates for mowing are not known, but the vegetation height was between 15-40 cm at the sites. The traps consisted of 0.8 l aluminum boxes (with the lid removed) sprayed with yellow paint. The traps were filled with water and a drop of detergent with the purpose to reduce surface tension of the water surface. The traps were placed at the sites during early morning at days with suitable weather conditions (sunny weather, 20-25°C), and the traps were emptied and removed from the sites four days later, at late afternoon/ evening. The collected material was stored in ethanol. All bees and wasps where identified to species-level.
At each collection plot (i.e. four plots per site), we made a rough estimation of three variables of importance for foraging and nesting of bees and wasps. These variables were based on four photographs of the vegetation in a 1x1 m wooden frame on the ground placed at four different positions at each collection plot; one at the position of the trap, two positions ten meters from the trap in each direction parallel to the road, and one position ten meters from the trap perpendicular to the road. These photographs where then used to estimate three variables: 1) the number of flowering plant species within the test square (as a qualitative measure of the foraging habitat), 2) the proportion of vegetation (in %) covered by flowering plants within the test square (as a quantitative measure of foraging habitat), and 3) the proportion of the area (in %) covered by bare soil within the test square (as a measure of nesting habitat availability). Based on the estimations from the four photographs, we then calculated an average of the three variables for each collection plot, which was used in the analysis. At each site we also subjectively assessed the inclination of the road verges on the eastern and western side of the road based on photographs taken at each plot. This was done to detect possible patterns in inclination between the two sides of the road that potentially could affect our results.

Statistical analyses
To test whether the species composition of the bee and wasp community differed between the eastern and the western side of the E4 highway, we used Canonical Correspondence Analysis (CCA; Oksanen et al. 2015). Trap catches for traps placed on the same side of the road at each site were pooled, as these trap catches may not be regarded as independent replicates. First, we performed a CCA on the entire community of bee and wasp species. We included two explanatory variables, (1) side of the E4 (eastern/western) to test for the barrier effect of the road and (2) the north coordinate to control for any potential north-south gradient in the data. Second, we performed a CCA only on the smallest species, as they could be expected to have a lower dispersal ability compared to the larger species (e.g. Greenleaf et al. 2007), and thus may be more sensitive to barrier effects. For this purpose, we categorized all species as either small or large. We used the size of a honeybee Apis mellifera as the criterion for a large species following Samnegård et al. (2015), as this species usually is considered as being mobile with a relatively high dispersal capacity (Schneider and Hall 1997;Beekman and Ratnieks 2000). We categorized all species larger than the honeybee, which has a maximum body length of approximately 16 mm (Amiet and Krebs 2014), as being large (mobile, 9 species), and all other species as being small (less mobile, 30 species) and performed a CCA on only the small species, with side of E4 (eastern/western) and the north coordinate as the explanatory variables. The data on size of the species were collected from identification literature and reference collections (See Suppl. material 1: S1 Aculeata species list and body lengths), and we consistently used the maximum body length reported.
As a complementary analysis, we investigated if the road verges on both sides of the E4 differed in vegetation characteristics, as these variables may have consequences for the interpretation of the results from the CCAs. For this purpose, we analyzed the number of flowering plant species, cover of flowering plants (in %), and cover of bare soil (in %) with the side of E4 (eastern/western) as the explanatory variable, using generalized linear mixed models (GLMMs) with site as a random factor.
Furthermore, we also examined if there were any general differences in abundance and number of species in bees and wasps between the eastern and the western side of the E4, and if the three vegetation variables could explain these two response variables. For this purpose, we used GLMMs with Poisson distributions and site identity as a random factor. However, initial analyses indicated that the number of flowering plant species was tightly correlated with the cover of flowering plants (Pearson correlation: r=0.73, p<0.001), and as the first provided a slightly better model fit (when tested separately) we decided to only include the number of flowering species, the cover of bare soil and the side of the road as explanatory variables in the final models.
All analyses were performed using R 3.2.2 (R Development Core Team 2015) with add-on libraries VEGAN (version 2.3-2) for the CCAs and glmmADMB for the GLMMs.

Results
In total, 111 individuals of 39 species of bees and wasps (Aculeata) were collected in the 20 yellow pan traps. The trap catch comprised 20 bee species (Apoidea), while the remainder consisted of various wasp species, such as Crabronidae (8 species), Pompilidae (6 species), Sphecidae (3 species), Tiphidae (1 species) and Vespidae (1 species) (see Suppl. material 1). We found several relatively rare species such as Lestica clypeata, Panurgus calcaratus and Psenulus brevitarsis, of which the first two previously were redlisted in Sweden (Gärdenfors 2005, Gärdenfors 2010. Panurgus calcaratus and Philanthus triangulum have been suggested to be indicators for species rich sandy habitats and dry meadows in Sweden (Swedish Board of Agriculture 2003, Karlsson 2008).
The canonical correspondence analyses (CCAs) showed a significant difference in species composition between the eastern and the western side of E4 (Table 1, Figure 2) when analyzing the whole community, and this pattern became even stronger (lower p-value) when we excluded the largest species (i.e. with a maximum body size >16 mm) (Table 1, Figure 2). In both analyses, there were also significant relationships between the species composition and the north-south gradient ( Table 1). The variation explained by the two explanatory variables was 37% for the whole community and 47% when excluding the largest species.  We found no significant difference in the number of flowering plants (p = 0.10), the cover of flowering plants (p = 0.16), or the cover of bare soil (p = 0.81) between the two sides of the E4. The abundance and number of species of bees and wasps did not differ significantly between the two sides either (Table 2). However, the abundance of bees and wasps increased significantly with an increasing number of flowering plant species, while the number of bee and wasp species did not show such a relationship ( Table 2). The cover of bare soil did not explain any of the two response variables. There was no clear pattern in the inclination of the road verges between the eastern and western side of the road: at site 1 both verges were inclining rather strongly towards the road; at site 2 the eastern road verge was flat while the western road verge inclined weakly towards the road; at site 3 both verges inclined weakly from the road; at site 4 the eastern road verge inclined towards the road while the western inclined from the road; and at site 5 both verges were rather flat.

Discussion
By analyzing differences in species composition along two sides of a large road we found indications that roads may act as barriers on the movement of bees and wasps. Even if the yellow pan traps were in most cases only separated with ~50 m across the road, the two analyses of the community composition both pointed in the same direction -the community composition differed between the western and eastern side of the road. Also when considering the fact that the trap catches where rather small (~10 individuals/site), the results show that species composition at sites situated on the same side of the highway were more similar to each other compared to the sites at the opposite side ( Figure 2). There were also differences in species composition in the north-south direction, which is more expected as this gradient spanned almost 9.5 km. The barrier effect seemed stronger when excluding the large species from the community (based on the lower p-value and the increased explained variation). We suggest this pattern to be a result of the fact that for flying insects large species in general have a higher mobility, thereby constituting better dispersers than small species Tscharntke 2002, Greenleaf et al. 2007). This suggests that larger species do not perceive large roads as barriers as much as smaller species. For example, among the large species in this study almost 45% are bumblebees, which are known to be strong fliers with the ability of quickly covering long distances (Walther-Hellwig and Frankl 2000, Steffan-Dewenter et al. 2001, Hagen et al. 2011) and possibly also flying relatively high above the ground compared to the smaller species. Even though we did not directly study insect movements and mortality, our results agree with previous studies on insects (Askling et al. 2006, Skórka et al. 2013. Askling et al. (2006) showed that smaller butterfly species (Aphantopus hyperantus and Maniola jurtina) were less likely to move between pastures on different sides of a large road compared to larger species (Pieris rapae, Gonepteryx rhamni and Anthocharis cardamines). The suggested reason for this was that the first two species have lower dispersal abilities and often tend to fly close to the ground, while the latter three species often move at higher levels above ground and are capable of quickly covering great distances across a landscape (Askling et al. 2006). Flying close to the ground also increases the susceptibility to motor vehicles (Soluk et al. 2011), which may explain why small butterfly species have been shown to be overrepresented among road-killed butterflies (Skórka et al. 2013). We are currently not able to determine whether our results are best explained by an increased mortality in smaller species, or if smaller species to a greater extent avoid crossing large roads or highways. To investigate this would require a study design with several roads, preferably of the same width, along a gradient in traffic volume (McKenna et al. 2001). However, even though we did not analyze the actual mechanism, our result suggests that large roads could have a fragmentation effect in the landscape by acting as barriers to insects, with consequences for community compositions that is critical enough to justify detailed research. Road verges have been highlighted as important grassland habitats, often rich in species of vascular plants and pollinators (Way 1977, Saarinen et al. 2005. Moreover, these habitats are often important for red-listed and rare species (Helldin et al. 2015). This was also confirmed by this study, where some rare or interesting species from a conservation perspective were found. However, the number of flowering plant species along the road verges had a relatively weak effect on the richness and abundance of bees and wasps, and only for abundance the relationship was statistically significant. Somewhat surprising, we found no effect of the cover of bare soil, which may indicate that the soil along the road verges is not optimal as nesting habitat for these species. However, it may also be due to the method used to quantify this variable. The occurrence of bare soil has a rather patchy distribution along roads, and it may therefore be more difficult to get a representative quantification of this variable based on four 1x1 m squares compared to e.g. the cover of flowers which usually is more evenly distributed. There was no clear pattern in the inclination of the road verges between the eastern and western side of the road, and it is therefore unlikely that the differences in community composition we found was caused by differences in inclination. However, inclination can most likely affect the microclimate of the road verges and it would therefore be interesting to investigate its effect on insect communities in future studies.

The influence of motivational factors on the frequency
of participation in citizen science activities ment and Perceived Competence was obtained by the group of respondents that participate a lot and the

Introduction
Citizen Science can be defined as the general public involvement in scientific research activities and has recently become a mainstream approach to collect information and data on many different scientific subjects (Miller-Rushing et al. 2012). The huge number of data collectors engaged in citizen science allows scientists to tackle questions that were previously out of their reach (Silvertown et al. 2011). With traditional scientific methods, the cost of such data collection would become a limitation due to budget and time constraints. Therefore, an increasing number of researchers have started to work with citizens, realizing that those directly involved in research activities exhibit a rapid increase in scientific literacy (Bonney et al. 2009;Lowman et al. 2008;Silvertown 2009). As such, citizen science has been recognized not only as an instrument for a given research experiment, but also as an education and outreach tool for researchers. The level of participation in citizen science studies is however remarkably different between regions and countries (Dierkes and von Grote, 2000;Forte and Lampe 2013). For citizen science projects to become successful, it is therefore essential to understand the motivations behind the different levels of participation of citizens. These motivations may be different, depending on the local historical and cultural background and among different societal groups.
Some studies aimed to identify the main motivations for people to participate in citizen science projects and have identified several reasons. The desire to learn more about scientific issues behind the project, the feeling that they are helping the environment and the enjoyment of developing activities in nature were recognize as important motivations to participate (Bell et al. 2008;Van den Berg et al. 2009;Raddick et al. 2010, Rotman et al. 2012). It was also described that getting to know other people with similar interests, making new friends, having the feeling that they are an active participant and co-owner of the project and gain recognition for their input and achievements were also reasons that encourage people to participate in citizen science projects (Bell et al. 2008;Van den Berg et al. 2009;Raddick et al. 2010, Rotman et al. 2012. In this study, we aim to analyze differences in motivations concerning gender, age, level of education and level of participation in one of the largest and longest running citizen science project in Portugal, the biodiversity web portal Biodiversity4all (www.biodiversity4all.org). The BioDiversity4All is a nationwide project that aims to increase citizens' biodiversity knowledge. Currently BioDiversity4All has nearly 2500 registered users, a network of 50 partners representing different citizen groups and other stakeholders and a validation panel already encompassing 49 taxonomic experts. The project has currently over 400000 observations of 7000 species, and includes nearly 98000 pictures associated to sightings. Users can add to the database either point species observations (sightings) or polygon areas for species occurrence which are later validated by taxonomic specialists (invited scientists or non-academic experts) and through this validation process, users progressively learn to identify and recognize local and national biodiversity. In order to understand the level of intrinsic motivation of Portuguese participants in this citizen science project, we tested the self-determination theory (SDT). SDT is grounded in the assumption that people have basic psychological needs to feel competent, autonomous and have a sense of belonging or relatedness to others (Ryan and Deci 2008). Autonomy involves feelings of willingness and choice in regards to activities undertaken; relatedness refers to feelings of closeness to other people; and competence involves feeling able to master challenges and having effective interactions with the environment (Katz and Assor 2007) (Figure 1). SDT predicts that, as a result of developmental experiences that engender competence, autonomy, and relatedness, individuals will advance towards more autonomous motivational orientations (in other words, the amount of self-determined motivation increases) (Katartzi and Vlachopoulos 2011). The most self-determined form of motivation is intrinsic motivation, representing the motivation to engage in an activity purely for the sake of the activity itself and because it is inherently pleasurable (Deci and Ryan 1985;Lepper et al. 1973). Intrinsic Motivation Inventory (IMI) is a multidimensional measurement instrument intended to assess participants' intrinsic motivations related to a target activity's subjective experience. It has been used in several experiments related to intrinsic motivation and self-regulation (e.g. Ryan 1982;Ryan et al. 1983;Plant and Ryan 1985;Ryan et al. 1990;Ryan et al. 1991;Deci et al. 1994). It assesses participants' Interest/Enjoyment, Perceived Choice, Perceived Competence, Pressure/Tension, Effort, Value/Usefulness and Relatedness. The category Interest/Enjoyment is the most direct measure (self-report) of intrinsic motivation. This category assesses the interest and inherent pleasure when doing a specific activity. Perceived Choice and Perceived Competence are theorized as positive predictors of intrinsic motivation and are related to the SDT innate psychological needs of autonomy and competence. Perceived Choice evaluates how individuals feel they engage in one activity because they choose to do it, and Perceived Competence measures how effective individuals feel when they are performing a task. Pressure/Tension, conceived as a negative predictor of intrinsic motivation, evaluates if participants feel pressure to succeed in an activity. Effort is a separate variable, which is important when taking into account motivation in specific issues and contexts. It assesses the person's Figure 1. Self-determination theory, illustrating basic psychological needs, defining features of the types of motivation and position in the relative autonomy continuum (adapted from Ryan and Deci, 2007). investment of his/her capacities in what he/she is doing. The Value/Usefulness category embodies the idea that people internalize and develop more self-regulatory activities when experience is considered as valuable and useful for them. Finally, Relatedness refers to the degree of a person's feelings connected to others and is used in studies where interpersonal interactions are relevant (Monteiro et al. 2015). The IMI statements are often slightly modified to fit specific activities. Thus, for example, a statement such as "I tried very hard to do well at this activity" can be changed to "I tried very hard to do well on these puzzles" or "…in learning this material" without effecting its reliability or validity. Concerning redundancy there are statements within the categories that can overlap. Making a randomization of the presentation of the statements makes these categories less evident to the respondents.

Survey instrument
We prepared a web survey that was sent to citizens registered in the BioDiversity4All project through the project's Newsletter's (Suppl. material 1).
The survey was composed of three sections. The first introduced the research and addressed survey ethics and data security. The second section asked about respondents' demographic and professional characteristics like gender, age, self-reporting level of participation in the project (from never participated to participate a lot), nationality, profession, and level of education. The third section (see Table 1 for all the questions in this section of the survey) was an adaptation of Fonseca and Brito's (2001) version of the IMI (McAuley et al. 1989).
The seven categories employed, (Table 1) although derived from the Intrinsic Motivation Inventory (IMI), were generated by reviewing the theoretical literature and relevant published instruments. In the present case, these were modified to refer to citizen science activities connected with biodiversity assessments. In the analysis seven categories were considered: Interest/Enjoyment (eight statements), Perceived Competence (nine statements), Effort/Importance (five statements), Perceived Choice (seven statements), Value/Usefulness (seven statements), Project Relatedness (six statements), and Group Relatedness (five statements). All motivational statements were rated on a seven-point Likert scale ranging from one (strongly disagree) to seven (strongly agree), with an intermediate score of four (moderately agree) (Munshi 2014). From the original IMI we excluded the category Pressure/Tension once is not expected to be felt by participants that do this activity in a volunteer basis and included the category Group Relatedness created according to the features of the project. Ryan and Deci (2000) describe relatedness as a sense of belonging and connectedness to the persons, group or culture disseminating a goal. Although the IMI analysis is designed to tap into individual motivation for doing a certain activity, the statements on the Group Relatedness category lend themselves readily to the assessment of the degree to which a person feels connected to other persons that do the same activities.
Because BioDiversity4All is a project developed in Portuguese language, the survey was only available in Portuguese even if the participants were from other nationalities. It was assumed that, if they had registered in the Portuguese platform, they could read in Portuguese.
The link to the web survey was sent in May 2015 to all the Newsletter's receivers of BioDiversity4All Project (N=1450), independently of their participation or not in the project. Five answering reminders were sent till October 2015.

Data analysis
The results from the survey were ranked and analyzed considering the questions referring to the participants' socio-demography and the IMI-related statements. All results describing the characteristics of participants and their motivation to participate were reported as a percentage of total responses.
To analyze differences between gender, of the average scores of the statements ranked on Likert-scales, we did a Mann-Whitney-Wilcoxon test. After calculating the medians with an interquartile interval (Q3-Q1) for age classes, levels of education and levels of participation, a multiple comparisons analysis was performed with the Kruskal-Wallis test (multiple comparisons and unbalanced sample sizes). Significant Table 1. IMI categories used in the survey with corresponding statements. The (R) after a statement is just a reminder that the score attributed is the reverse of the participant's response on that particular statement. differences between average scores were determined for α≤0.05. All the statistical analysis was performed using R 3.1 (R Development Core Team 2014).
Finally, we performed a cluster analysis to group participants according to similarities in the answers they provided. We used hierarchical agglomerative clustering with Ward method (Murtagh and Legendre 2014) and performed the cluster analysis using the package Cluster of R 3.1.

Results
We received 149 survey responses corresponding to 10.3 % of the Newsletter's receivers. Most of the responses were given by Portuguese citizens 92.6% with the remaining representing six other nationalities: Brazilian, Spanish, British, French, Dutch, and Swiss.
From the total amount of responses 77 were given by males (51.7%) and 72 by females (48.3%) and participants' ages varied between 19 and 71 years old with an average of 43.5 ± 11.4 (Figure 2). Concerning the level of education 83.1% had higher education (44.6% bachelor degree, 25.7% MSc, 12.8% PhD) and 16.9% high school (Figure 2).
Respondents that had registered in the project and only occasionally participate were responsible for largest (55.7%) fraction of the survey's responses, followed by those that had registered in the project but never participated (30.2%). Of the remain-  ing a small fraction (12.1%) regularly participate and very few (2.0%) showed a high degree of participation ( Figure 2). Concerning their professional activity, 28.9% of the respondents to the survey have education related jobs and 25.5% have environmental related jobs. Considering all survey participants, the highest IMI scale-score was obtained by the category Project Relatedness, with an average of 5.8 out of 7, followed by Perceived choice and Value/Usefulness with an average of 5.7. Interest/Enjoyment had an average of 5.3, Group Relatedness an average of 4.7 and Perceived competence an average of 4.5. The lowest average obtained referred to Effort/Importance with 3.8. In the correlation analysis of the different IMI scores, Interest/Enjoyment and Value/Usefulness, Interest/Enjoyment and Project Relatedness, and Value/Usefulness and Project Relatedness were strongly correlated (Table 2).
No statistical differences were found between genders, age classes and levels of education for the averages in any category of analysis. However, levels of participation were significantly different for all categories except Interest/Enjoyment (Table 3 and Figure 3).
The cluster analysis of the answers given by the participants supports the differences of motivations of the respondents with different levels of participation (Figure 4). The first cluster group was composed of people that never participated or that participate only occasionally. The third group included people that participate a lot and most of the people that participate regularly.
The highest value of Interest/Enjoyment and Perceived Competence was obtained by the group of respondents that participate a lot and the lowest by the ones that never participated. For Effort/Importance, the lowest value was obtained by the group that participates occasionally and the highest by those who never participated. For Value/ Usefulness, Project Relatedness and Group Relatedness, the highest value was obtained by the ones who show high participation levels and the lowest by the ones that never participated. For Perceived Choice the highest value was obtained by the ones that participate regularly and the lowest by those that never participated.

Figure 3.
Percentage of answers of the IMI categories, rated from 1 (strongly disagree) to 7 (strongly agree), with an intermediate score of 4 (moderately agree), for the four groups of people with different levels of participation in the project (never participated, participate occasionally, participate regularly, participate a lot).

Figure 4.
Hierarchical agglomerative clustering of the answers given by the participants, with Ward method. The symbol represents the level of participation in the project (never participated -, participate occasionally - participate regularly -, participate a lot -).
Concerning the group of people that never participated, the lowest IMI was Perceived Competence and the highest was Perceived Choice. For all the other groups, the lowest IMI was Effort/Importance and the highest was Project Relatedness.

Discussion
In this study, we wanted to assess citizens' engagement and meaningful experience in citizen science projects, using motivational measures. This study revealed lessons of interest for citizen science projects when participants' motivations is concerned, in a country with limited culture of public participation. Assessment of intrinsic motivations in countries with higher levels of engagement with biodiversity and participation in citizen science, could present different results and a comparative analysis would be an interesting approach.
Analyzing survey respondents, the majority of participants have higher education, a fact which is not representative of the Portuguese reality (only 16.5% of Portuguese people have or are undertaking higher education, PORDATA 2015). Moreover, the age groups <25 and ≥65 were the ones with less answers to the survey (5% each); one reason might be that these are the groups with less participants in the project, or that these are the groups showing less willingness to answer to web surveys. For a general characterization of respondents, we also included questions about nationality and professional activity. The survey was developed for Portuguese speakers and this may have hampered people from other nationalities to participate. Several participants from other nationalities collaborate with the project either through the Portuguese project or through the international platform. Some of these participants are residents in Portugal and presumably speak Portuguese however, less than 8% of survey respondents were from other nationalities. Although the professional activities of respondents are diverse, 54.4% of respondents have education or environmental related jobs. The demographic factors of nationality and profession were just used to characterize respondents and not to test the motivational differences. Nationality because the project has an inherently national scope and answers to profession because they were too generic to allow any conclusions.
A high percentage of respondents had registered in the project BioDiversity4All but never participated. When we analyzed the responses to IMI categories given by groups with different levels of participation, we found that people who never participated were the ones responding more differently compared to other groups. This group shows the lowest levels in all categories except Effort/Importance. This might indicate that those people do not have intrinsic motivations for participating in such a project. Of these people, some registered after a project presentation, a media news or a launch of a contest but did not pursue with their participation. A possible lesson to draw from these results is that extrinsic motivations may be needed to foster participation in these cases, while creating mechanisms to increase competence, autonomy and relatedness on participants, to drive more autonomous (self-determined) motivations.
Frequently, citizen science projects use extrinsic motivation instruments to induce citizens' participation, such as incentives, certificates of recognition and challenges, which stimulate people's interest in the project (Dickinson et al. 2012). Nevertheless, it is important to include mechanisms to foster intrinsic motivations in order to create continued support and involvement in citizen science initiatives after these initial extrinsic motivations erode (Cialdini 2008). For example, one could use contests and prizes that include educational material, feedback on the effort already invested, group activities, interacting with a similar community and different ways of participating, increasing perceived choice.
In contrast with the respondents that never participated a small percentage (2%) participate a lot. This is not unexpected regarding results from other citizen science projects. In the Wikipedia project, with one million registered users, about 10% contribute with ten or more entries and about 0.5% contribute to a large number of tasks to keep Wikipedia running (Tapscott and Williams, 2008). The group of respondents that participate a lot had the highest levels of intrinsic motivation, scoring highest in all categories except Effort/Importance and Perceived Choice.
These findings are aligned with past research on intrinsic motivation which has focused on identifying and examining the activity-level psychological factors that promote or inhibit the development of intrinsic motivation. This approach has yielded important insights, some of which that (1) enjoyment is positively related to competence valuation (i.e. the degree to which one cares about performing well at a given activity; Elliot et al. 2000;Goudas et al. 1995;Harackiewicz and Manderlink 1984;Reeve and Deci 1996;Sansone 1989;Tauer and Harackiewicz 1999), and (2) enjoyment is positively related to the degree to which activities are perceived to be ''optimally challenging''-not too easy and not too difficult (e.g. Harter 1978;Keller and Bless 2008;Moneta and Csikszentmihalyi 1996). Stated more generally, the degree to which the potential rewards of ongoing activity engagement are realized would seem to be dependent on the degree to which attentional resources are devoted towards these potential rewards.
Early experiments showed that positive feedback enhanced intrinsic motivation relative to no feedback (Boggiano and Ruble 1979;Deci 1971) and that negative feedback decreased intrinsic motivation (Deci and Cascio 1972). Deci and Ryan (1985) linked these results to the need for competence (White, 1959), suggesting that events such as positive feedback provide satisfaction of the feeling of competence, thus enhancing intrinsic motivation, whereas events such as negative feedback tend to thwart the feeling of competence and thus undermine intrinsic motivation. That is why it is understandable that people who participate a lot in the Biodiversity4All project had the highest levels of Perceived Competence. The feeling of competence leads them wanting to participate more.
These results indicate that citizen science projects should nurture participants with positive feedback and adapted participation modes to their level of competence. This may yield higher levels of motivation to participate and foster intrinsic motivation.
Project Relatedness and Value/Usefulness were the highest scoring IMI categories for all groups except those who never participated. People tend to value the feeling of relationship and trust in the project, moreover since they feel that the project has an important mission to accomplish.
A note should be given about the category of Perceived Choice. Most respondents feel they had a high level of Perceived Choice which is in line with the voluntary nature of the project. However, we have students participating in the project and some contests with schools and scouts which may explain why some respondents may feel that they had no choice in their participation.
With the cluster analysis we wanted to confirm similarities in the answers given by different respondents to find, if people with the same level of participation, have comparable intrinsic motivations and in fact, we detected the expected result.
In conclusion, in recent years much has been written on communication and recruiting participants for citizen science projects (Dickinson et al. 2012;Roy et al. 2012;Silvertown et al. 2013). However the results from our work show that working deeply on people's involvement is fundamental to increase and maintain their participation on citizen science projects. If, for initial recruitment and in countries with low participation culture, mechanisms of external motivation may be necessary, to guarantee higher levels of long term participation, citizen science projects should foster intrinsic motivations which can be done by incorporating in project design experiences of relatedness, capacity building, positive feedback and adapted participation modes.