Pollination by wild insects and other animals is an intermediate regulation ecosystem service (ES), which is to say a service rendered by the ecosystem to itself, necessary for the provision of final ES, those from which the anthropogenic system directly benefits. In this case, the ecosystem performs functions of intermediate ecological regulation in support of the final service of crop biomass provision, on which fertilisation and agricultural productivity depend eventually (Cory et al. 2016). It is known that pollination can increase yields, quality and stability of fruit and seed crops (Bommarco et al. 2012). In fact, it was found that 87 out of 115 crops of global importance (~ 70%) benefit from animal pollination (Klein et al. 2007). The Italian bee fauna is one of the richest in the world in relation to the surface of our country. The latest official list includes 944 species (Pagliano 1995) belonging to six of the seven worldwide known families. In 2018, the IUCN Italian Committee, on behalf of the Ministry of the Environment, drew up the Red List for Apoidea. In total, there are 151 species of native bees in Italy for which there are sufficient data and indications of decline: out of 151 species assessed according to the rigorous criteria of the IUCN, 34 are those with different levels of threat. Nevertheless, there is a need for better spatial assessments of pollination as few studies are based on Species Distribution Models (SDMs) and these are at relatively coarse spatial resolution. Especially in highly fragmented landscapes like the ones in Italy, this resolution may fail to cover important habitats, such as hedgerows, small pastures and forests, with the risk to obtain biased pollinator distribution maps. Three most important elements are needed for pollinators to live in a habitat: suitable places to nest, sufficient food (provided by flowers) and availability of water near their nesting sites (Vallecillo et al. 2018). For this reason, the pollination model focuses on resource needs, flight behaviour of wild bees (the most relevant pollinator group for most crops (Free 1993) and, finally, on pollination service associated with certain crops (Capriolo et al. 2020). We rely on biophysical modelling to identify the areas in which pollination is most likely to occur, based on the overlapping presence of pollination-dependent crops (Klein et al. 2007) and areas with high pollinator habitat suitability (Martínez-López et al. 2019). The model quantifies the supply of the pollination service, on the basis of a relationship between the suitability of the habitat and the activity of pollinating insects. Habitat suitability expresses in synthesis key characteristics of the ecosystem: its ecological suitability for nesting, floral availability (representative parameter of the insect's foraging function) and its proximity to water sources such as rivers, lakes and waterways in general (CCN 2019). The notion of habitat suitability can be represented with an index with values ranging from 0 (habitat with no suitability) to 1 (very high habitat suitability). The activity of pollinators is simulated as a black spherical model representing a function of the average annual temperature (°C) and the average annual solar irradiation (W× m-2), as they can greatly affect insects foraging and, in turn, the number of active individuals (Martínez-López et al. 2019 for complete equations). For this reason, we use the black spherical model as a proxy of the abundance of pollinators. Modelling the relationship between active number of pollinators and habitat suitability describes the activity/microclimate relationship, from which we obtain a representation of the ES supply which refers to a potential flow (Capriolo et al. 2020). In order to represent the demand side of the ES, the model considers the production in metric tonnes of pollination dependent crops, multiplied by the pollination dependency rate (Klein et al. 2007). Additionally, only crops with a dependency rate greater than 0.2 have been considered for the modelling exercise. As a result, when ES supply (occurrence of insect activity and habitat suitability) and demand (production of pollinated-dependent crops according to the crop-specific dependency rate) are overlapped, we identify first the met demand and then the ES use. Due to a lack of national maps on agricultural production spatially distributed, we used the most recent and accessible spatialised data on crops extent and yield (Monfreda et al. 2008) which include 173 types of crops, but only six crops were taken into account: almond, apple, apricot, melon, pear and watermelon. Production data of the six crops for the years 2012 and 2018 were obtained from the Italian agricultural yearbooks and are illustrated in Table 1 (INEA 2013, CREA 2020). The economic valuation was carried out by applying a market-based method using the “farm gate”*5 price at the producer site. The increase in agricultural production can be measured as a share of agricultural production attributable to pollination flow and, in monetary terms, by multiplying the output of the use of the service with the ‘farm gate’ basic price, for each of the different crops (CREA 2020, INEA 2013) and for both years. This component of production would not exist in the absence of ecosystem service and, therefore, represents as a whole the additional value deriving from the presence of pollinators (Vallecillo et al. 2018). These quantifications converge in the construction of the SEEA-EA accounts and can be collected in dedicated tables, both in biophysical and monetary terms (UN 2017, UN 2014a, UN 2014b, UN 2021a, UN 2021b).

The flood risk regulation model maps and evaluates the service by identifying areas at risk of flooding (FHP: flood hazard probability) through an index consisting of a first climatic-weather parameter (CCN 2019) which gives an indication of presence, intensity and volume of precipitation, a second parameter that contains a topographical humidity index (Kirby and Beven 1979, Manfreda et al. 2011) which combines precipitation with the flow of water accumulation descending along the slope obtained from digital elevation model (DEM) and, lastly, a temperature parameter of the wettest quarter for the years (Trenberth et al. 2003, Hijmans et al. 2005, Utsumi et al. 2011). A further parameter needed to determine the supply component of the model is the Curve Number (CN), an indicator of potential run-off, which estimates the ability of vegetation and soils to retain excess run-off from precipitation. CN is a function of land cover, soil hydrological groups and slopes in some contexts (SCS 1985, Zeng et al. 2017) and is formed by combining each CN soil groups value with Corine Land Cover (CLC) classes for both 2012 and 2018 years (Suppl. material 1). CN scores range from 0 to 100, with higher scores indicating higher run-off. In order to evaluate the potential run-off over the six years, we consider the variation between the CN values associated with the actual CLC soil groups and CN values associated with the bare soil (or bare areas), where the run-off is greater, in line with a hypothetical scenario of substantial loss of service. This change is then multiplied by the FHP to obtain the biophysical indicator of the Flood Regulation Supply (FRS), which identifies the areas subject to flooding (equation 1). This relation allows us to compare the current probability of flood hazard (FHP) with that of a possible loss of service (bare areas) and to obtain an estimate of the actual flow of the service.

\(\text{FRS} = \text{FHP} \times (CN_{bareareas}-CN_{soilsgroup})\) (eq.1)

Eq.1 represents the supply (potential flow/capacity) of the ecosystem service and is then spatially described (Fig. 1) by the areas on which the ES has an impact. Overlapping these areas on the land-cover and use data, which instead constitute the service demand and identify the assets and values benefiting from the service (or potentially at risk in case of service loss), it is then possible to estimate the use of the service. This model is the result of a simplification of those already published on a global or continental scale (Stürck et al. 2014, Ward et al. 2015), but has the advantage of being easily replicable even in contexts with little data available. The methodology adopted for the economic evaluation uses different monetary damage functions estimated for different land-cover and use classes (Huizinga et al. 2017). Damage values were initially elaborated from existing studies across some EU countries and the average damage value per land-cover and use class was applied to other EU Member States scaled to GDP per capita. The damage functions are built on observations from nine countries and, being nationally homogenous, they do not account for regional differences. The direct economic impacts (equation 2) are derived from depth-damage functions that express the damage cost in EUR/m² as a function of the flood water depth (in metres) for different classes of land-cover and use. Specifically, it is a function of the land-use type (damage value per each land use), the level of the damage (damage factor, based on water depth) and the extension of the flooded area by land-use type.

\(\sum_{i=1}^{n} \text {damage value}_i \times \text{damage factor}_i \times \text{extension of flooded area}_i\) (eq.2)

where i = land use type.

Monetary damage values (Suppl. material 2) are available for almost each class (few classes in Suppl. material 2 are missing with respect to Table 6) as they have been adjusted to the CLC classes identified as economic assets currently protected by the service. The economic values derived from this direct economic impact analysis represent the potential avoided damage arising from the presence of the ecosystem service.

Nature-based tourism-related services are defined as the ecosystem contributions, in particular through the biophysical characteristics and qualities of ecosystems, that enable people to use and enjoy the environment through direct, in-situ, physical and experiential interactions with the environment. Nature-based tourism-related services are a part of the recreation-related services. According to the SEEA EA (UN 2021c), recreation-related services include services to both locals and non-locals (i.e. visitors, including tourists)

The model is inspired by previous works mapping outdoor recreational activities in Europe (Paracchini 2014) and is aimed to identify nature-based tourism in areas with a high naturalistic value. Data available on this kind of tourism are very scarce, both globally and nationally: in this purpose, one of the most useful sources of information for this process has been identified in a previous work (Paracchini et al. 2014), which has the merit of having involved an in-depth survey on visits to protected areas in more than 50 countries.

The model presents four fundamental assumptions:

1. only foreign tourism is considered (i.e. incoming or inbound) (Balmford et al. 2015) and, as a consequence, the role of domestic tourism is ignored in this specific study, although already addressed for Italy in a previous and recent work (Capriolo et al. 2020);
2. all the “naturalist” tourists go initially to protected areas as represented in the Balmford dataset (Balmford et al. 2015), while other tourist destinations, based on nature, are not considered – this, in particular, is probably a strong assumption considered that most of nature-related attractive spots are not covered by PA in the Italian context;
3. this approach disregards other country-specific touristic attractive elements, such as the weight of cultural tourism in the Italian context;
4. an incomplete dataset is used and, therefore, results should be interpreted with a high level of caution.

Since Italy was not significantly represented in the Balmford's dataset and, therefore, there was no chance to carry out a regression analysis that could include data for our country, it was decided to use the dataset only to run the econometric analysis and estimate the portion of nature-based tourists and then proceed with a distribution of nature-based visitors over all areas of potential naturalistic interest, going beyond the perimeter of the formally-protected areas. This has been done considering the conformation of our territory which, unlike other countries, where fewer large naturalistic areas coincide with protected areas, often shows landscape and naturalistic value even in non-protected areas or in areas that do not belong to the Natura 2000 network. However, it is desirable to proceed towards an improvement of this estimate as soon as more accurate local data on the number of visits to protected areas become available.

The model assigns an indexed and normalised value to each grid cell on the basis of specific features in order to obtain the supply:

• the naturalness of different land cover types – every landcover type is associated with a level of hemeroby*6 to assess the degree of human influence in that area;
• presence of protected areas;
• distance from mountains and water bodies - to identify attractive areas for outdoor recreation;
• indicator of fauna richness linked to the presence of vertebrate species.

Potential supply is then weighted considering accessibility, intended as the distance from cities and road infrastructures, which means that the most accessible hot-spots for outdoor recreation are also the most likely to be visited by “naturalist” tourists. The result is a map of areas identified by weighing the attractive nature of each cell within the grid with its accessibility and this can be used to spatially spread the number of tourists at national level over the attractive areas (Fig. 2).

We have built a simple univariate regression model to quantify the relationship between visits to protected areas and natural parks (dependent variable) with data on tourists arriving for leisure and recreational purposes from the United Nations World Tourism Organization*7 (World Tourism Organization 2019). The univariate regression model, which intends to provide a starting analysis tool capable of estimating annual visits to protected areas as a function of incoming tourists associated with leisure activities, indicates that about 50% (precisely 50.05%) of inbound tourism is involved in nature-based tourism. The estimate obtained has a coefficient of determination R² equal to 0.835, which is an index of excellent confidence for regression diagnostics. The addition of further variables in a multivariate model was explored, including the percentage of land area covered by protected areas within each country and the percentage of the country economy represented by tourism, but these additional variables did not add predictive value nor were they statistically significant. Further details on the regression results and the test run on the econometric model are summarised in Suppl. material 3.

We use the proportion obtained from the regression analysis described above to predict the number of tourists at national level attracted by naturalistic elements and then to calibrate the relative attractiveness map. Once the number of inbound visitors for Italy, linked to the enjoyment of nature has been estimated, this has been spatially distributed on the basis of the characteristics of the landscape and of the operating dynamics of the model described above for the supply.

In order to obtain the share of expenditure relevant for nature-based tourism, the total expenditures for tourism in the economy is multiplied first by the percentage of travelling for holidays, leisure and recreation purposes and then by the percentage of this group of travellers engaging in nature-based tourism. The monetary value associated with nature-based tourism (MVNT) for the year 2018 is presented here below (equation 3):

\(\text {MVNT} = \text {Total expenditure on inbound tourism} = \frac{\text {Leisure tourism and leisure time}} {\text{Tourism made for personal reasons} } \times 50.05\) (eq. 3)

We have considered data on total inbound foreign tourism expenditure, available from the United Nations tourism database (World Tourism Organization 2019). We have extracted the share attributable to leisure tourism and leisure time by applying the rate of visits with this specific destination preference to the total of visits made for personal reasons (i.e. those inclusive of both holiday tourism or linked to weddings funerals, etc., but net of travel for business). Then, the share of recreational tourism associated with experiences of enjoyment of nature is computed using the percentage estimated in the regression analysis (50.05%).

For the year 2012, in light of the absence of WTO data on tourism for Italy, the monetary value associated with the nature-based tourism (MVNT) is calculated as follows:

\(\text{MVNT} = \text{TE} \times \text{HT} \times \text{VNT}\) (eq. 4)

Where,

TE = total expenditure for inbound foreign tourism;

HT = percentage of holiday tourism expenditure: 62% (CISET 2014);

VNT = percentage of visitors for nature-based tourism: 24.3% (CISET 2014).

The pollination service supply table for the selected six crops (Table 1) describes the spatial extension of land-cover areas that provide the service for the two years considered jointly to the supply variation over this period (Table 2). The map in Fig. 3 reports, as an example, the use value spatialised for the “apple” crop, while Table 3 reports the biophysical use values expressed in metric tonnes for the six crops considered. Observing the production of pollinated crops and comparing it to the total production for the respective years, we have an indication, in percentage, of the contribution of pollination to these fruit crops production that remains substantially unchanged (6.10% in 2012 and 5.97% in 2018) (Table 4). Table 5 highlights the monetary use value that estimates the economic benefits (prices 2018) which arise from the pollination service. Fig. 4 provides, again as an example, the monetary use value of the apple pollination service. The monetary values are placed in the areas where the service is present with different intensity: pollination use index values range from 0 to 1 and are marked with different colours.

The use table represents the accounting format of ES and can be constructed in physical or monetary terms (UN 2014a, UN 2014b, UN 2017, UN 2021a, UN 2021b). In particular the first column of Table 6 shows the assets that were considered. Then the next columns show the area of each asset that potentially benefits from the flood risk regulation service, other than the change over time in absolute and percentage terms. The economic sectors that benefit from the ecosystem service are considered in Table 7. Table 6 and Table 7 show the change over time between 2012 and 2018. The monetary value index in Fig. 5 highlights areas of greatest potential economic risk in case of service loss.

Table 8 represents the use table of ES and reports the tourism rate in percentage for the two years, in the first column and data in monetary terms, in the other two columns, related respectively to the total value of inbound tourism expenditure and the total value of inbound tourism, based on nature. Fig. 6 shows the spatial distribution of the monetary value index for nature-based recreational tourism service (price 2018) for the two years taken into consideration.