Application of a Quality-Specific Environmental Risk Index for the Location of Hives in Areas with Different Pollution Impacts

Abstract

Honeybees and hive products have long been recognized as excellent bio-indicators, as they provide valuable information on the state of the environments concerned. As yet, however, no tool is capable of contextually assessing the level of pollution of different environmental matrices in order to create maps of areas in which to place hives. In this paper, a possible method of creating a regional map to support the implementation of targeted monitoring plans on beekeeping is described. We obtained and identified related-activity pressure factors, which were subsequently combined by means of a multi-criteria approach through the analytic hierarchy process method (AHP). The different levels used were drawn up by first attributing qualitative values; these were converted into quantitative values through scoring elaborations and pair comparisons and then elaborated and standardized by means of different techniques in order to create an index with a spatial distribution of five risk classes throughout the region. To verify the correct execution of the procedure, a consistency ratio method was implemented on this index and validated the reliability of the application as the main source of information for sampling activities on beekeeping products. Creation of the specific environmental risk index enabled us to construct a map displaying the areas of greatest impact on beekeeping activities and a representation of the cumulative effects generated by the different pollutants in the air, water, soil and subsoil compartments. This index may, therefore, constitute an essential tool to support beekeepers in choosing sites for their apiaries.

1. Introduction

As economic well-being increases, global pollution levels are rising, to the detriment of the quality of human and animal health [1,2]. The development of population centers over the years in the Campania region has led to a spatial continuity between agricultural areas, urban areas and industrial areas. In the vicinity of large population centers, past urban planning has not provided for the proper division of space like in other parts of the world, resulting in the loss of natural areas and increased pollution levels. Some studies have shown that different levels of exposure to various pollutants cause localized damage to specific organs of the human body [3,4]. Severe environmental pollution also impacts ecosystems and the climate, leading to extreme climate variations and disastrous weather phenomena [5,6]. One of the possible methods of determining pollution by metals, organic pollutants and microplastics is to use bioindicators as tracers [7,8,9]. The use of hymenoptera, an order of insects encompassing more than 120,000 species worldwide, as bioindicators [10,11], or of beekeeping products [12,13,14] is well established. Indeed, bee products provide indications of the state of the environment and the chemical contamination present [15]. In particular, the honeybee (Apis mellifera Linnaeus, 1758) is a ubiquitous, easy-to-breed organism that displays great mobility. Moreover, as the bee’s body is covered in hairs, it picks up materials and particulates encountered in the environment. Bees can, therefore, bioaccumulate large amounts of materials from soil, vegetation, air and water. These characteristics make the honeybee a valuable bioindicator and an ideal agent for the easy monitoring of large areas, even in regions where infrastructure is scarce [16]. A study was carried out to calculate an index that would enable us to identify the most appropriate locations for hives in order to minimize the effect of pollution on bees and their products. Determining the most suitable areas for the placement of hives, and the mapping of these areas, is a very important task, as it can reduce colony losses [17] and, at the same time, ensure the necessary pollination services. Indeed, beekeeping contributes to the sustainability and conservation of natural resources through the pollination of a large number of native and cultivated plants [18]. In this context, the data present in the literature on the levels of pollutants in environmental matrices (water, air, soil and plants) in the Campania region of Southern Italy were processed [19]. In this way, a synthetic and specific environmental risk index with regard to beekeeping was produced in order to classify the environment of the Campania region. This index was created by using multi-criteria spatial-indicator construction techniques. These techniques are usually used to determine the choice of a specific site for the realization of industrial activities, as in the case of a photovoltaic field in Ecuador [20] or the determination of an optimal site for a waste-storage facility in Iraq [21]. The problem of understanding how multiple matrices (or levels) influence a certain activity has been analyzed in several studies using multi-criteria decision analysis (MCDA) techniques [22]. MCDA techniques are also used to evaluate spatially distributed data [23] and, especially in studies with a strong spatial component, in combination with GIS systems [24,25,26]. These approaches focus on a particular type of MCDA technique, the analytic hierarchy process (AHP) in a GIS environment. The synergy between the two techniques facilitates the creation of maps for specific activities, such as planning new coffee plantations in view of climate change [27], or as a support tool for soil-erosion risk assessments [28]. In the field of beekeeping, however, only one such study has been carried out, in the Calabria region of Southern Italy; this focused on the creation of suitability maps through the “fuzzy overlays” technique [29]. The present study aimed to draw up a practical approach to the location of hives in areas with different pollution impacts by using the AHP method in conjunction with spatial analysis techniques in a GIS environment. The novelties and innovations presented in this work are in the use of satellite imagery to assess air pollution levels, producing maps of pollutant concentrations in the lower troposphere and the use of modelist tools to create a regional map to support the implementation of targeted monitoring plans concerning beekeeping.

2. Materials and Methods

The index constructed was derived from the combination of four levels, chosen because they represent the main compartments where different types of pollution occur (Figure 1).

For this purpose, layers of different components were used, combined through an analytic hierarchy process (AHP), a type of multi-criteria decision analysis (MCDA), in a GIS environment. The AHP first decomposes the problem into a hierarchy of sub-problems. The decision-maker then evaluates the relative importance of the various elements by means of pairwise comparisons. The AHP converts these evaluations to numerical values (weights or priorities), which are used to calculate a score for each alternative. A consistency index measures the extent to which the decision-maker has been consistent in her/his responses. The use of decision-support techniques in conjunction with GIS systems is a topic that has been consolidated in recent years. Indeed, in the geospatial field, multi-criteria analysis techniques are increasingly used in synergy with GIS software to generate synthetic indices (Figure 2). In this study, the different layers were elaborated by using QGIS and ESRI^® ArcGIS^® 10.7 software with plug-ins for performing multi-criteria analyses.

In the first step, it is necessary to identify the possible areas of greatest impact on beekeeping activities. All potential sources of toxic substances that can interfere with beekeeping products were considered. Thus, four qualitative/quantitative variables, uniformly distributed throughout the territory of the Campania region, were identified: one level for the representation of land uses, one level for the representation of surface water bodies, an information level regarding air quality and an information level associated with the presence of metals in the soil. The different levels were processed in stages using a sequential methodology in order to obtain uniformity in the various workflows of each level used. For simplicity, the general scheme was divided into two stages: a first stage was used to process, correct and standardize the data to be used in the model, and a second stage was used to apply AHP analysis to the standardized data (Figure 3).

In the first stage, the geographical levels to be used in the model were identified, and they were processed and standardized by transforming them into four raster files (intermediate output) using various spatial geostatistics tools, such as inverse distance weighting (IDW) interpolation and statistical classifications using the “natural-breaks” method (Figure 4).

In the second stage, the actual AHP analysis took place, in which the intermediate levels obtained (output of the first stage and input of the second) were processed by applying weights and scores (Figure 5).

The data used were verified regarding correctness from the point of view of both numerical and spatial information. In the case of a qualitative variable, a method of scoring by means of a panel [30] was used to obtain a quantitative variable. All layers were processed in order to obtain a representation of coverage classes by using 9 numeric classes and the natural-breaks method [31]. This procedure is mandatory when the analytic hierarchy process (AHP) method is used in a GIS environment [32]. To assign weights to the individual levels, it is necessary to use the pairwise comparison method; this method is used in cross-comparisons between different pairwise-compared levels, in which numerical values of relative importance are assigned (

Table 1. Priority scale utilized in pairwise comparisons.

Value	Definition	Explanation
1	Equal importance	Two criteria contribute equally to an objective.
3	Moderate importance	The experience or judgment of a specialist slightly favors one variable over the other.
5	Strong importance	The experience or judgment of a specialist strongly favors one variable over the other.
7	Very strong or proven importance	One criterion is greatly favored over the other; the effect is demonstrable.
9	Extremely important	The evidence favors a criterion as highly as possible.
2, 4, 6, 8	Intermediate values between scale values (when a compromise is required).

).

Considering the scarcity of specific information in the literature, a method was devised to assign a weight value to each variable, thus, defining a priority scale among the identified levels. The database of apiaries (Banca Dati Apistica Regionale (BDA-R)) present in the regional territory was used to determine possible spatial interferences between levels (Figure 6).

Figure 6 shows the apiaries distributed throughout the Campania region, whose area is approximately 13,700 square kilometers. With these values, it was possible to determine the regional average density of the apiaries present (about 0.5 apiaries per square kilometer), which was relatively uniformly distributed. The priority scale among the levels was managed according to the following scheme:

land use > water bodies > air pollution > soil hazard.

The land use level was assigned a higher weight because it has a greater impact in terms of spatial interference, representing the territory in which the apiaries are located. The water bodies level was assigned a lower weight than the land use level; it represents the surface hydrographic network, which may or may not be located near the apiaries, contributing to a moderate interference compared to the apiary points. The air quality level was assigned a lower weight than the first two levels. This variable is represented as a column throughout the troposphere. Considering that the interference is limited to a few meters, a small portion of the column was assigned a moderate interference. The soil hazard level is considered to have a lower weight than all the other levels because the interference associated with this level comes solely from the analytical values of trace elements present in the soil and subsoil. To carry out the pairwise comparisons, a double-entry table (

Table 2. Weight attribution by pairwise comparisons.

	Land Use	Water Bodies	Soil Hazard	Air
Land Use	1	4	4	3
Water Bodies	1/4	1	3	4
Soil Hazard	1/4	1/3	1	1
Air	1/3	1/4	1	1

) was used; the values in the upper triangular matrix were entered in accordance with the priority scale; the values on the main diagonal were 1 for all cells; and the values in the lower triangular matrix were obtained from the reciprocal of each element in the upper triangular matrix.

Table 2 was used to obtain a “priority vector”, whereby a percentage weight was attributed to each individual variable by software output in order to determine how much each one contributed to obtaining the index (Figure 7).

The different levels were processed in the GIS working environment by means of QGIS software and ArcGIS desktop with plug-ins for performing multi-criteria analysis (extAhp 2.0—developed by Oswald Marinoni). In the following sections, the creation of the levels of the different variables is described in detail.

2.1. Quality Level of Surface Water Bodies

To assess the quality of surface water bodies, the analyses presented in the monitoring plan for surface waters [33], which were carried out by the ARPAC (Regional Agency for Environmental Protection in Campania region) in the years 2013–2020, were used. These data are the result of the processing of environmental surveys which qualitatively described the chemical and ecological status of the water at selected points along the watercourses. The classification of the chemical status of surface water bodies is based on the presence of chemicals defined as priority and hazardous substances, such as heavy metals, pesticides, industrial pollutants and endocrine disruptors, among others. These chemical substances are categorized into three groups based on their level of hazard, and distinct environmental quality standards (EQS) are established for each of them for the analysis matrices (water, sediments, biota) where they may be present or accumulate. Non-compliance with the established EQS for each of these substances implies an assessment of “not reaching good chemical status” for the water body. On the other hand, compliance with EQS leads to a “good chemical status” evaluation. The regulation in this regard has been established by Directive 2008/105/EC, updated later by Directive 2013/39/EU, and is implemented in Italy by legislative decree as of 13 October 2015, no. 172. The evaluation of the chemical status of surface water bodies is an important measure for the protection of environment and public health, as contaminated water can pose serious risks to human health and the ecosystem as a whole. Therefore, it is essential to maintain high environmental quality standards to ensure the protection of public health and the sustainability of the environment. The ecological status of inland surface waters, according to Italian legislation, is an indicator that describes the quality of the structure and functioning of aquatic ecosystems. The objective of ecological quality, established by Directive 2000/60/CE, refers to the capacity of the water body to support well-structured and balanced animal and plant communities, which are fundamental biological tools to support the self-purification processes of the water. The classification indicators are represented by biological quality elements, physical-chemical elements in support of biological elements, chemical elements in support of biological elements and hydro morphological elements. The classification is defined by calculating the average value for each parameter analyzed in each year of monitoring and using the worst status obtained in the reference period. The ecological status is defined by the element that is in the worst class according to the general principle, called “one-out, all-out”, of Directive 2000/60/CE. To support this level, it was necessary to use the database of surface water bodies reported on the Campania region website (https://www.arpacampania.it/web/guest/acque-superficiali, accessed on 9 January 2023). Processing of the levels was carried out in a GIS environment by means of the QGIS software. The ARPAC sampling points for the determination of the qualitative chemical and ecological status were used to identify the stretches of water to be evaluated. This involved using a multiplicative criterion and normalizing qualitative values to quantitative scores, which were assigned on a scale from 1 to 10 (

Table 3. Ecological and chemical quality scale.

Ecological	Ecological Score	Chemical	Chemical Score
Dry River/N.A.	1	Dry River/N.A.	1
Excellent	2	Good	5
Good	4	Bad	10
Average	6
Poor	8
Very Poor	10

N.A.: not available.

). On multiplying the score value by 100, these scores became the inputs for creating a buffer zone for each stretch of water analyzed. This calculation yielded the radius of the buffer zone in meters in order to create a sort of “interference band” around each stretch of water.

To determine the quality of the surface watercourse, the overall score was calculated by multiplying the ecological score by the chemical score. Subsequently, an interpolation technique (inverse-distance weighting) was implemented by inserting the previous score as an interpolation attribute. This elaboration yielded a raster file (Figure 8), which was used directly in the general model.

2.2. Quality Level of Land Use

Data on land cover, land use and the transition between different categories are frequently required for the formulation of sustainable land management and planning strategies. Indeed, they can provide information to support decision-making processes at the EU, national and local levels and enable the effectiveness of environmental policies to be verified. In this study, land-cover data derived from the 2018 edition of the Corine Land Cover (CLC) project (https://land.copernicus.eu/pan-european/corine-land-cover/clc2018, accessed on 19 December 2022) were used; this project was conceived at the European level, specifically for the detection and monitoring of the characteristics of coverage and use of the territory. The CLC outputs are based on the interpretation of satellite photographs by the national teams of the participating states (member states of the European Union and states that cooperate), according to a methodology and a standard nomenclature with the following characteristics: 44 classes at the 3rd hierarchical level of the Corine nomenclature, with some thematic insights at the 4th level regarding Italy. A panel test was distributed to technicians through the Delphi method [30] for the attribution of a quality judgment of the land-cover classes. For each CLC class, these judgments were based on the attribution of a judgment of interference with beekeeping activities on a 9-point qualitative scale, with 1 indicating minimum interference and 9 maximum interference (

Table 4. Corine Land Cover quality score attributions.

Impact	Score	Description
Lowest	1	3112—Woods with a prevalence of deciduous oaks (turkey oak and/or downy oak); 3113—Mixed woods mainly with mesophilous and mesothermophilous broad-leaved trees; 3115—Predominantly beech woods; 3122—Forests mainly of montane and oro-Mediterranean pines (black pine and larch, Scots pine, loricate pine); 3124—Woods with a prevalence of larch and/or stone pine; 3125—Woods and plantations predominantly of non-native conifers (Douglasia, insigne pine, white pine); 3131—Mixed woods with a prevalence of deciduous trees; 331—Beaches, dunes, sands; 332—Bare rocks; 333—Areas with sparse vegetation; 334—Areas covered by fires.
Very low	2	3212—Natural grassland with trees and shrubs; 3231—Sclerophyllous vegetation; 3232—Low scrub and garrigue; 324—Transitional woodland shrub; 3241—Young stands after cutting (and/or clear cuts).
Low	3	131—Mineral extraction sites; 3111—Woods with a prevalence of holm oak and/or cork oak; 3114—Predominantly chestnut woods; 3132—Mixed woods with a prevalence of conifers; 3211—Natural grassland prevailingly without trees and shrubs; 323—Sclerophyllous vegetation; 411—Inland marshes; 421—Salt marshes; 511—Water courses; 512—Water bodies; 521—Coastal lagoons; 523—Sea and ocean.
Quite low	4	3117—Woods and plantations mainly of non-native broad-leaved trees (robinia, eucalyptus, ailanthus, etc.); 3121—Woods mainly of Mediterranean pines (stone pine, maritime pine) and cypress trees.
Medium	5	241—Annual crops associated with permanent crops; 242—Complex cultivation patterns; 3116—Broad-leaved forests with a prevalence of hygrophilous species (forests with a prevalence of willows and/or poplars and/or alders, etc.).
High	6	141—Green urban areas; 224—Other permanent crops; 231—Pastures; 243—Land principally occupied by agriculture, with significant areas of natural vegetation.
Quite high	7	112—Discontinuous urban fabric; 123—Port areas; 124—Airports; 133—Construction sites; 142—Sport and leisure facilities; 244—Agro-forestry areas.
Very high	8	111—Continuous urban fabric; 221—Vineyards; 223—Olive groves.
Highest	9	121—Industrial or commercial units; 1211—Industrial areas; 122—Road and rail networks and associated land; 132—Dump sites; 2111—Arable land predominantly without dispersed (line and point) vegetation; 2112—Arable land with scattered (line and point) vegetation; 212—Permanently irrigated land; 222—Fruit trees and berry plantations.

).

This scoring yielded a qualitative land-cover map (Figure 9), which was then converted to raster-file format by means of GIS software to be used as input for the AHP model.

2.3. Level of Air Pollution

Air quality is essential to characterizing the state of the environment and the different pressures on regional territories. In addition to being a secondary source of contamination, air constitutes a possible route of migration of pollutants, owing to the mechanisms of volatilization and dispersion. The pollution of air by numerous pollutants may be of natural or anthropogenic origins (emissions from industrial activities, vehicular traffic, domestic heating and agricultural practices). To obtain a synthetic indicator of air quality, an elaboration of different bands from the Sentinel-5 Precursor satellite was used; Copernicus Sentinel-5 Precursor is a mission dedicated to monitoring our atmosphere and is the result of close collaboration among the European Space Agency, the European Commission, the Netherlands Space Office, industries, data users and scientists. The mission consists of one satellite carrying the TROPOspheric Monitoring Instrument (TROPOMI); this space-borne, nadir-viewing imaging spectrometer covers wavelength bands between ultraviolet and shortwave infrared and can detect daily concentrations of formaldehyde (CH₂O), carbon monoxide (CO), nitrogen dioxide (NO₂) and sulphur dioxide (SO₂) in the tropospheric column (https://s5phub.copernicus.eu/dhus/#/home, accessed on 12 October 2022). Today, few studies use TROPOMI in order to acquire information on air pollutants by using different bands of the on-board spectrometer [34,35]. To download the images and reduce the calculated values, a script was created by means of the Google Earth engine to elaborate all the images from the year 2021 and aggregate these by means of a value-reducer function in order to obtain minimum, median and maximum values at the 95th percentile (https://earthengine.google.com/, accessed on 12 October 2022); this approach yields a single image for each band of air pollutants, which represents the ratio between the moles of substance analyzed and the moles of air (

Table 5. The sentinel 5P range of values for each pollutant.

	$\frac{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{N}{\mathbf{O}}_2}{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{a}\mathbf{i}\mathbf{r}}$	$\frac{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{S}{\mathbf{O}}_2}{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{a}\mathbf{i}\mathbf{r}}$	$\frac{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{C}{\mathbf{H}}_2\mathbf{O}}{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{a}\mathbf{i}\mathbf{r}}$	$\frac{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{C}\mathbf{O}}{\mathbf{m}\mathbf{o}\mathbf{l} \mathbf{a}\mathbf{i}\mathbf{r}}$
minimum	0.0000814	0.0006046	0.0002453	0.0310240
maximum	0.0001271	0.0048747	0.0006101	0.0515890
average	0.0001042	0.0027396	0.0004277	0.0413065

).

A simplified AHP method, without weighting and pairwise comparisons, was implemented to obtain a single air pollution value; this was discretized by means of the natural-breaks classification, with a value of 1 indicating the lowest concentration and 9 the highest concentration (Figure 10).

2.4. Potential Hazard Level of Soil

Regarding soil and subsoil, an index modeled on the probability of exceeding the legal limits of pollutants was used. This indicator was constructed within the framework of an agreement between Federico II University and Istituto Zooprofilattico Sperimentale del Mezzogiorno [19]. This index was derived from the analysis of about 4000 soil samples, which were examined for the presence of 53 chemical elements. These data were analyzed by means of spatial statistics models (e.g., inverse distance weighting, kriging, geographically weighted regression) and enabled us to reconstruct continuous areas of concentration that covered the entire regional territory and to estimate the probability of exceeding the legal limits or reference values in areas not covered by sampling [36]. From this analysis, we created a map of the Potential Hazard, which was representative of the number of areas exceeding the CSC (contamination threshold concentrations). This enabled us to identify those areas of the regional territory where multiple inorganic contaminants were present; these were then assigned a relative score according to the level of contamination and were subdivided into classes by means of the “natural-breaks” method (Figure 11).

2.5. Statistical Verification of Consistency

To form a statistical evaluation of the correct execution of the AHP method, the consistency index (CI) was used [37]. In this index, the value of λ represents the maximum eigenvalue of the matrix and n the size of the matrix itself:

$$CI=\frac{({\lambda }_{max}-n)}{(n-1)}$$

(1)

If the value of the CI is zero, then the matrix is consistent; if it deviates from n, then the matrix is not perfectly consistent; however, the method used accepts a low degree of inconsistency, as this does not affect the validity of the result obtained. In the first approximation, the maximum eigenvalue of matrix A can be evaluated by referring to the average of the consistencies relative to the individual variables; the result is a maximum eigenvalue of 4.27271, which is close to the dimension n of the matrix A. By means of the CI, it is possible to define the Random Consistency Index (RI), the values of which are available as a function of the size of array A. In this case, for n = 4, the RI value is 0.9 (

Table 6. The Random Consistency Index.

N	3	4	5	6	7	8
RI	0.58	0.9	1.12	1.24	1.32	1.41

). At this point, the Consistency Ratio (CR) of matrix A can be calculated:

$$CR=\frac{CI}{RI}$$

(2)

In the present case, the ratio is 0.1, indicating that matrix A is consistent.

3. Results

Elaboration of the specific environmental risk index by means of the AHP method yielded a map of the points of greatest impact on beekeeping activities. The representation of the raster map was split into five risk classes (using the natural-breaks subdivision method) in order to obtain the best separation between the classes (Figure 12). The five risk classes represent different territorial realities. Classes five and flour, which represent “very-high risk” and “high risk”, respectively, correspond to areas of high population density, major industrial centers and widespread pollution of water and air. Risk class three, corresponding to “average risk”, represents buffer zones, corresponding to areas with high population density but few heavy industries and mostly agricultural enterprises. Risk classes two and one represent “low risk” and “very low risk” values, respectively; these are areas with limited industrial or agricultural activity, low population density and very low or zero air and water pollution levels.

Figure 12 reveals a geographical pattern; in the southern and eastern areas of the Campania region, there are more class one and class two areas, while in the northwest and center there are more class four and class five zones. Overall, however, class one and class two zones account for most of the area of the Campania region, and classes four and five together constitute less than 19% (Figure 13).

By making a map overlay between the calculated index and the database of apiaries, it was possible to observe a uniform spatial distribution of all the apiaries surveyed in the regional territory (Figure 14).

On those points, a territorial statistic was made by counting the points falling into the different levels of risk, and the relative percentage was calculated (

Table 7. Distribution of apiaries in the different risk levels in the Campania region.

Risk Level	Number of Apiaries
Lower risk	18.5%
Low risk	36.8%
Medium risk	26.8%
High risk	14.0%
Higher risk	3.9%

).

By overlapping the layers, it can be observed that the number of beehives in high-impact areas is significantly lower compared to other areas. These data could be explained because highly polluting anthropogenic activities could coexist in high-impact zones, resulting in concentrations of pollutants to which bees are particularly sensitive. The index helps us evaluate the cumulative effect of pollutants from different environmental compartments. An example of a cumulative effect that may arise is as follows: in the land-cover level compartment, there may be the presence of industrial areas, roads and railways, intensive agricultural areas and a lack of vegetation that is functional to bee nutrition in the interference area corresponding to 3 km². In the water quality compartment, there may be the presence of potentially toxic elements whose legal limits are exceeded. In the soil compartment, there may be the presence of inorganic pollutants in the soil and subsoil. In the air quality compartment, there may be the highest concentrations of gaseous pollutants, such as SO₂ and NO₂, due to the geomorphological features of the territory that can lead to photochemical smog effects in particular climatic conditions [38]. This complex condition that arises could create imbalances with significant repercussions on the fitness of honeybees. To validate this hypothesis, studies will be necessary to evaluate the fitness of bees using precision sensors.

4. Discussion

Different papers cited in the introduction part use multi-criteria techniques in applications for geographic information. These techniques use multiple factors in the same analysis with the goal of generating a single indication to support decision-making processes. The difference between multi-criteria techniques is in the way that the criteria are chosen and processed. The aim of this study was to construct an environmental risk index through which to map the Campania region in order to provide a tool to support the implementation of beekeeping-related environmental monitoring plans. Among the papers with the greatest similarities to the topic of beekeeping is Zoccali et al. In fact, the similarity to our work is to generate sensitivity maps for beekeeping activity, whose objective is to evaluate the most suitable areas to carry out beekeeping activity. The application of geostatistical methods aimed at creating beekeeping-specific maps is becoming an increasingly used technique; previous studies have attempted to identify the most suitable areas in which to undertake beekeeping activities. Zoccali et al., for example, argue that honeybees are of fundamental importance to the environment and the economy; therefore, they developed a rapid method of identifying areas suitable for beekeeping, with a view to maximizing productivity and reducing the risk of colony loss. To achieve this goal, they superimposed different levels, such as road networks, temperature, hydrographic networks and CLC, using a fuzzy overlay analysis that did not use the assignment of weights on the different factors used. In this study, data were acquired on all environmental matrices involved in bee production (water, air and soil) in order to gain an overall view of the environment in which bees operate, with a major focus on pollution, taking into account the cumulative effect generated by the different pollutants in the air, water, soil and subsoil compartments. One limitation of the study stems from the regional nature of some of the data used. Indeed, the soil hazard level was derived from work carried out only in Campania, and the quality of surface water was assessed on the basis of surveys carried out by the Campania region’s Environmental Protection Agency (ARPA CAMPANIA). To minimize this local approach, the CLC and Sentinel 5P were taken into account. In the future, a standard technique of territorial mapping should be drawn up using nationwide data that are easily retrievable. Further applications should try to use as many global coverage tools as possible so that the method can be easily replicated.

5. Conclusions

The specific risk index can be a very valid tool, both for the construction of bee maps from a production point of view and for the identification of those areas where environmental contamination can be a serious problem for bees and beekeeping products. Environmental pollution is one of the leading causes of mortality in the beekeeping sector and is, therefore, one of the main problems that must be addressed by the scientific community and policymakers. In recent decades, many agricultural practices have been designed to improve agricultural production; however, these have resulted in long-term pressures on the natural environment, negatively affecting native pollinators and reducing pollination services. In addition to developing general governmental agricultural policies, local and regional authorities could take measures to develop beekeeping on a local scale. Such environmental planning approaches are currently gaining increasing attention; there is a need for collaboration between stakeholders and researchers to create incentives so that decisions taken by individuals, communities, companies and governments can promote widely shared values that are compatible with ecosystem services. With this in mind, the environmental index can be a suitable predictive model for assessing environmental impacts and, thus, an essential tool to support beekeepers in choosing sites for their apiaries. In bio-monitoring work, this index tool facilitates the exact identification of points and locations to be considered according to the level of risk. Additionally, it could also provide indications regarding the wellbeing of bees related to pollution issues by predictively monitoring the fitness of honeybees.