**An Investigation of Takagi-Sugeno Fuzzy Modeling for Spatial Prediction with Sparsely Distributed Geospatial Data**

**Robert Thomas 1, Usman T. Khan 2, Caterina Valeo 1,\* and Fatima Talebzadeh <sup>1</sup>**


**Abstract:** Fuzzy set theory has shown potential for reducing uncertainty as a result of data sparsity and also provides advantages for quantifying gradational changes like those of pollutant concentrations through fuzzy clustering based approaches. The ability to lower the sampling frequency and perform laboratory analyses on fewer samples, yet still produce an adequate pollutant distribution map, would reduce the initial cost of new remediation projects. To assess the ability of fuzzy modeling to make spatial predictions using fewer sample points, its predictive ability was compared with the ordinary kriging (OK) and inverse distance weighting (IDW) methods under increasingly sparse data conditions. This research used a Takagi–Sugeno (TS) fuzzy modelling approach with fuzzy c-means (FCM) clustering to make spatial predictions of the lead concentrations in soil. The performance of the TS model was very dependent on the number of outliers in the respective validation set. For modeling under sparse data conditions, the TS fuzzy modeling approach using FCM clustering and constant width Gaussian shaped membership functions did not show any advantages over IDW and OK for the type of data tested. Therefore, it was not possible to speculate on a possible reduction in sampling frequency for delineating the extent of contamination for new remediation projects.

**Keywords:** fuzzy modelling; marine sediment; Takagi–Sugeno; ordinary kriging (OK); inverse distance weighting (IDW); spatial predictions

#### **1. Introduction**

The release of pollutants into the natural environment has been a problem of global concern since the beginning of the industrial revolution. Highly toxic persistent environmental pollutants often occur in marine harbour sediments as a result of industrial practices around the world and pose a significant risk to human health [1]. The major contributor to exposure of humans to contaminated sediment is through the ingestion of contaminated food as a result of bioaccumulation through the food chain. Chronic exposure to contamination can lead to infertility, birth defects, impaired child development, diabetes, damage to the immune system, disruption of hormonal function, and cancer. Seafood in the most common exposure pathway of sediment contamination to humans; therefore, for the protection of human health, the remediation of these contaminants in aquatic environments is of the upmost importance [2–5].

The first step required for remediation is an accurate assessment of the spatial distribution of contamination in order to ensure the most effective and least costly remediation. Spatially continuous data are required to delineate the boundaries of unsafe levels of contamination and to determine the volume of contaminated material to be removed. However, in aquatic environments, point samples are generally collected on a predetermined grid spacing, and the contaminant concentrations are spatially interpolated to provide a continuous surface that may introduce uncertainty. Spatial interpolation methods are traditionally grouped into deterministic and geostatistical methods. The most commonly

**Citation:** Thomas, R.; Khan, U.T.; Valeo, C.; Talebzadeh, F. An Investigation of Takagi-Sugeno Fuzzy Modeling for Spatial Prediction with Sparsely Distributed Geospatial Data. *Environments* **2021**, *8*, 50. https:// doi.org/10.3390/environments8060050

Academic Editor: Sílvia C. Gonçalves

Received: 23 April 2021 Accepted: 28 May 2021 Published: 29 May 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

used deterministic method is inverse distance weighting (IDW), which uses a function that estimates values at un-sampled points through the linear combination of known sample values weighted by the distance between them [6,7]. IDW is a simple method requiring minimal modeler input [8]; however, because it is based solely on distance, it often performs very poorly with sparsely distributed geospatial data. The most commonly used geostatistical method for spatial interpolation is kriging, which employs a semi-variogram that plots the semi-variance between points against the distance between them [6,9]. From the semi-variogram, it is possible to determine the range of spatial dependence of sampled points that can be used to determine a value at an unknown location [10]. However, the accurate estimation of a semi-variogram is complicated, computationally expensive, and can introduce modeler bias. Furthermore, kriging has been shown to have a significant smoothing affect, where areas of high pollutant concentration could be missed. This is not ideal as an underestimation of the pollutant concentrations could lead to an increased risk to human health. Both deterministic and geostatistical methods have one major commonality, that is, the greater the sample density, the greater the accuracy of the spatial interpolation [11–13]. However, the sampling cost of sediment in marine environments and the analytical assessment cost for dioxins are extremely high. Therefore, obtaining an adequate number of samples to achieve an acceptable resolution during interpolation may not be possible, and this high cost may be prohibitive to remediation projects.

There is a separate family of data-driven predictive methods that utilize fuzzy set theory [14], which have been proven to be a suitable method for the prediction of sparse nonlinear data and have been used for many applications of spatial estimation in geoscience [11,15–17]. Fuzzy set theory [14] provides a convenient way of describing the degree of belonging (membership *μ*) of an element to a set, between 0 (no belonging) and 1 (complete belonging). For example, let U be an ordinary set with elements {*x1*, *x2*, ... , *xn*} and *A* be a fuzzy subset of U, in which the elements *xi* have degrees of membership (belonging to *A*-) given by a membership function *μ<sup>A</sup>* (*xi*) <sup>=</sup> *<sup>α</sup>*, which dictates that an element *xi* has a degree of membership *α* to fuzzy set *A* , where 0 ≤ *α* ≤ 1 (see Figure 1).

**Figure 1.** Example of a (**left side**) crisp and (**right side**) fuzzy set, and their respective membership functions.

#### *The Takagi-Sugeno Method*

System modeling techniques that employ fuzzy set theory are commonly referred to as fuzzy modeling [15,18,19]. One of the most commonly used fuzzy modeling techniques for spatial estimation is the Takagi–Sugeno (TS) method [15,18,20]. TS fuzzy modeling breaks

down the input data space into a number of fuzzy regions and creates a linear function for each region. It is advantageous for spatial modeling because of its transparency and interpretability [21]. The degree of belonging that an un-sampled point has to each region in the data is used to predict a value at that location. The partitioning of the input space into fuzzy regions is achieved though fuzzy clustering, which is the foundation for the fuzzy spatial modeling using the TS method. Fuzzy clustering differs from traditional crisp data clustering in that in fuzzy methods each element (sample point) can have a degree of membership to multiple clusters within the data. Each cluster is defined by a cluster center, which has a value calculated from the membership-weighted average of the members of that cluster [22,23]. For spatial modeling, data are clustered in the three-dimensional product space defined by the Cartesian map coordinates (*x*, *y*) and the magnitude of a pollutant concentration (*p*) (see Figure 2).

**Figure 2.** Theoretical visualization of fuzzy clustering in the three-dimensional product space.

After clustering, each data point has membership to one or more cluster centers, depending on the pollutant concentration. The clusters are then partitioned onto the *x* and *y* Cartesian product space and a membership function is generated for the *x* and *y* coordinate axis of each cluster. The membership functions from each cluster are then used to determine the degree of membership an un-sampled location has to the different clusters within the data (see Figure 3).

Ultimately, the combination of these degrees of membership is used in solving a pollutant concentration at an unknown location. Once the data have been clustered, they are subjected to a rule-based fuzzy inference system (FIS) that makes inferences about un-sampled geographic locations based on their membership to the clusters within the data. A single rule is introduced for each cluster using conditional "IF-THEN" statements. FIS uses input variables referred to as antecedents for each rule; in the TS fuzzy model, the antecedents are the fuzzy sets generated from the clustered data. The antecedents are subjected to IF-THEN statements to produce a weighting for the prediction from each rule based on membership to that rule (cluster). The output of each rule is referred to as a rule consequent. When a rule in the FIS is executed, if the antecedent is unaffected by the IF-THEN condition, that rule is skipped and the next rule is executed. If the IF-THEN condition produces a consequent, then that rule is deemed to have fired or been executed.

**Figure 3.** Example of partitioning membership to cluster centers to the map coordinate axes and applying a membership function (represented by solid lines) for three fuzzy clusters to determine membership based on geographic location.

Each rule in the TS fuzzy model uses *x* and *y* Cartesian coordinates as inputs to solve the consequent (output variable) for each rule in the form of a linear equation to produce a crisp value [23]. The form of the general first order TS model is shown in Equation (1).

$$\begin{array}{l} \mathcal{R}\_{i} = \text{if } \mathbf{x}\_{1} \text{is } A\_{i1} \text{ and } \mathbf{x}\_{2} \text{ is } A\_{in} \text{ then } p\_{i} = a\_{i1}\mathbf{x}\_{1} + a\_{i2}\mathbf{x}\_{2} + b\_{i} \\ \qquad \qquad i = 1, \ldots, K \end{array} \tag{1}$$

where *Ri* is the *i*th rule; *x*<sup>1</sup> and *x*<sup>2</sup> are the antecedent variables (Cartesian coordinates *x*, *y*); *Ai*<sup>1</sup> and *Ai*<sup>2</sup> are fuzzy sets for the *i*th rule and the respective coordinate axes; *pi* is the *i*th rule output; *K* is the number of rules (clusters); and *ai*1, *ai*2, and *bi* are the unknown model parameters that must be solved. This is accomplished by least squares regression using the data in each cluster. The model output for a given input is obtained through the aggregation of all utilized rule consequents weighted by their membership to the respective rules [15,19,20]. The TS method's simplistic nature and interpretability make it advantageous compared to other FIS [24].

Applications using TS fuzzy modeling for spatial estimation have shown its predictive capacity to outperform kriging [19,25]. The advantage of TS modeling is that it allows a complex data surface to be broken down into more easily modeled individual fuzzy surfaces. Clustering and rule-based methods are then used to estimate a value at unknown locations by solving the intersection of the contributing surfaces based on the degree of belonging an unknown location has to each surface.

Fuzzy set theory has been cited as a useful tool for modeling under sparse data conditions [26]. Fuzzy modeling can take advantage of this additional information to produce an adequate contaminant distribution map for remediation planning using fewer point samples. Lowering the number of samples required would reduce the initial cost of new remediation projects and may lead to more remediation projects being undertaken in the future, eventually leading to a cleaner environment and to lower anthropogenic health risks [27].

This research employs a spatially distributed marine, soil geochemical dataset to compare the predictive abilities of the TS fuzzy model, IDW method, and ordinary kriging method (OK) using a sparse dataset. The geochemical signatures present in the data are the result of naturally occurring mineralization and are not related to anthropogenic sources. The data are useful, because the indicator geochemical signatures associated with the mineralization are analogous with several pollutants of concern (POC) in harbour settings. The spatial distribution of the contaminates in the data are very similar to that from anthropogenic sources, as both generally originate from point sources and spread out into the surrounding environment.

#### **2. Materials and Data**

The dataset used in this research is comprised of 1535 spatially distributed sample points (illustrated in Figure 4), each consisting of a 51 element spectral array with notable POCs of arsenic, mercury, and lead. Although no two spatial datasets would have identical statistical distributions, the challenges faced during spatial interpolation are similar regardless of the data type. Therefore, the use of this data still provides an adequate initial test of the relative predicative ability of the TS fuzzy model under increasingly sparse spatial data density.

**Figure 4.** Spatial expression of the data used for this research, Datum: NAD83, UTM Zone 8, Eastern Yukon Territory.

The 1535 soil samples were collected on an approximately 50 m by 50 m grid spacing at depths ranging from 0.1–1.2 m below the ground surface. The samples had compositions ranging from well to poorly developed soils to glacial sediment. The sample area was approximately 5000 m × 4000 m in size and was in the Yukon Territory, Canada. The samples were collected between 17 August 2015 to 22 August 2015, and between 16 June 2016 to 25 June 2016 by a sampling team from Archer, Cathro, and Associates Limited under contract for ATAC Resources Limited. ALS Minerals in Vancouver BC analyzed the samples using metallurgical assay with inductively coupled plasma and atomic emission spectrometry. The spatial coordinates of each sample were recorded in the spatial datum NAD83 UTM Zone 8. Each data point was characterized using an easting and northing UTM coordinate in meters. Because of its relevance to human health, its high rank as a POC in urban settings and its availability in the soil data, lead was selected to test the efficacy of TS fuzzy modeling compared with OK and IDW in this work.

As OK and IDW are well-developed tools, this research employed the GIS platform ArcMap 10.4.1 (ESRI, 2011) and used the Geostatistical Analyst Package to perform the OK and IDW methods using optimized parameters. Using the spatial predictions from ArcMap, model validation was performed using MATLAB R2017a (MathWorks, 2017). The fuzzy model was also developed and validated, and its results were compared with OK and IDW using MATLAB.

In order to assess the predictive ability of the TS fuzzy, kriging, and IDW methods, the sample data were split into training and validation sets, where the training data were used to predict the known pollutant concentrations at the locations in the validation set. To assess the predictive ability of the models under increasingly sparse spatial data conditions that simulated fewer samples, the amount of training data used were incrementally reduced from 75%, to 66%, to 50%, to 33%, to 25%, to 12.5% and labelled as sets A–F, respectively. At each increment, the data not used for training were utilized for validation.

To determine the optimal number of clusters, the fuzziness parameter (*m*) for the FCM clustering and the membership function width (*σ*) for the FIS, a pseudo-optimization was performed. This was accomplished by iteratively increasing the number of clusters and at each iteration varying *m* through a reasonable range, and at each iteration of *m* varying σ through a reasonable range. The model performance results from all combinations were then assessed and the combination of *m* and *σ* that yielded the most accurate prediction for the respective validation sets were retained as optimal.

The key things that determined the competency of a spatial interpolator were its predictive ability, its ability to handle data of different types and variance, and the smoothness or abruptness of the surface generated [6]. The performance metrics used to assess the model performance of TS fuzzy, OK, and IDW modeling were the coefficient of determination (R2), root mean squared error (RMSE), and the mean absolute error (MAE). R2 is commonly used and provides a metric by which to assess the variance of the model prediction. A higher R2 value indicates that a higher proportion of the total variation of the prediction is explained by the model [19]. An R2 value equal to 1 would indicate a perfect prediction. Each model result for MAE, R2, and RMSE received a score between 1 and 10. The mean score then determined the overall performance of that model result, with 10 being the best and 1 being the worst. The ranges for the scoring system were chosen arbitrarily based on the performance ranges observed from each metric. This approach aimed to produce enough separation between scores to draw realistic inferences on the differences in performance.

#### **3. Results**

When determining the optimal number of clusters, it was observed that the performance of the of TS fuzzy model would reach an initial peak or plateau and subsequent model runs using a higher number of clusters would not produce any further increase in score. In the interests of maintaining minimal complexity within the model and reducing the risk of over clustering, the lowest number of clusters that yielded the highest model performance was retained as optimal. In cases where more than one TS fuzzy model iteration with the optimal number of clusters but different parameters produced the same mean score, the result with the lowest mean absolute error (MAE) was deemed optimal. The MAE was selected because of its prevalence in assessing the accuracy of spatial predictions [28]. Figure 5 displays the results from the selection of the optimal number of clusters for data reduction increments A and B for subset 1. In each case, the highest mean score from each number of clusters tested is displayed. Therefore, for each cluster, the results displayed are the optimal performance with the optimal *m* and *σ* for that number of clusters.

**Figure 5.** Example results from the determination of the optimal number of clusters for data reduction for (**a**) increment A and (**b**) increment B, for subset 1.

In all cases, a clear initial peak was reached and dictated the optimal number of clusters for that particular training set. After identifying the number of clusters that yielded the highest model performance, the combination of fuzziness (*m*) and membership function width (*σ*) that contributed to the highest mean score were identified. Figure 6 displays an example of the results from determining the optimal *m* value. In each case, the highest mean score for each *m* value tested with the optimal number of clusters is displayed.

**Figure 6.** Example results from the determination of *m* using the optimal number of clusters for data reduction for (**a**) increment A and (**b**) increment B.

In all instances, a single *m* value produced the highest model performance for each increment of the data reduction and the *m* value had a visibly large effect on the performance of the TS fuzzy model. The optimal *σ*, which produced the highest model performance, was determined to occur over a range and was not overly sensitive. In general, a range of 80 m to 150 m produced similar mean scores. Figure 7 displays a 2D visual representation of the clustered data, the cluster centres, and the spatial distribution of membership to the cluster centres. The concentration of lead at a validation point was solved using its membership to each of the clusters in the data.

**Figure 7.** (**a**) Visual representation of clustered data for training increment A and (**b**) training increment F.

Finally, using the optimal model parameters, the fuzzy model was trained using each increment's training data (A–F); then, the lead concentrations at the spatial locations of each increment's validation sets were solved using the model. The predictions were then compared with the actual values and the performance metrics applied. This same methodology was also applied to the OK and IDW methods, such that the performance of the three models could be compared. Figure 8 displays the performance results for each model at all of the data reduction increments.

**Figure 8.** *Cont.*

**Figure 8.** Ordinary kriging (OK), inverse distance weighting (IDW), and the Takagi–Sugeno (TS) fuzzy model from each increment of data reduction performance metrics (**a**) mean score, (**b**) mean absolute error (MAE), (**c**) coefficient of determination (R2), and (**d**) root mean squared error (RMSE).

#### **4. Discussion**

In general, the three models tested here performed very similarly. For all three methods, data reduction increment B produced the highest performance scores. This is because this increment's training set has a much higher range, variance, and kurtosis than the respective validation set, making for simpler predictions closer to the mean. The discrepancies identified between the training and validation sets appear to have a large effect on the performance of the models—specifically, between increment A and B, where, for increment A, the kurtosis and range are higher in the validation set and for increment B, the kurtosis and range are higher in the training set. The effect of statistical differences between the training and validation sets does appear to be significant. However, in a real-world situation, if fewer samples are collected and analyzed, the results could have a lower range and kurtosis; so, for the purposes of simulating spatial data sparsity, the training and validation sets used here are useful because they illuminate weaknesses in the models tested. This proves that if the samples collected fail to capture the overall variance of a system, their spatial predictions will be poor, which Li and Heap (2011) found when reviewing 18 studies using different spatial interpolators.

Throughout the data reduction increments, on average, all of the methods' performances decline, which is the desired result of the evaluation. For increments A–F, IDW consistently outperformed the other methods, and the TS fuzzy model outperformed OK, but only by a small margin. To further quantify these observations, the performance of the individual metrics from the optimal TS fuzzy model results and for the OK and IDW methods from increments A–F are also plotted against the amount of training data used in each increment (shown in Figure 8). The differences between the three methods are so subtle that it is difficult to draw any significant conclusions from the results. As is consistent with the mean scores, based on these performance metrics, IDW appears to produce the best results, with the TS fuzzy model outperforming OK by a small amount. For the lead predictions made in this analysis, OK, IDW, and the TS fuzzy model all appear to have a significant smoothing effect at all data reduction increments, which agrees with previous research [12,29]. Previous research comparing IDW and OK found that in general, the two methods perform very similarly, and that IDW often produces a lower error, and this research agrees with that finding [29–31]. Furthermore, previous research has also shown that TS fuzzy modeling has the ability to outperform kriging [25,32,33] as it has in this study, albeit by a small amount. Although not clear in this research, the TS fuzzy model

did produce a higher mean score than OK for data reduction increments C–F. Although the three models are not using identical distributions in their predictions of unknown points, it is possible that the similarity in the results is simply a function of the inherent similarities within the models themselves. However, further research would be required in order to quantify this.

As fuzzy clustering is at the core of the TS fuzzy modeling approach, it can be surmised that the performance of the method may be limited by the quality of the fuzzy clusters or by their spherical shape, which agrees with findings from the literature [15]. Therefore, further investigation into other clustering methods is warranted. Additionally, the data used in this research have a negatively skewed distribution. A factor that was not considered was normalizing the distribution of the data prior to the analysis, which would impact the performance of the models. However, for testing the relative performance of the TS fuzzy model against IDW and OK for spatial prediction using incrementally less spatially distributed data, the results do show that the TS fuzzy model in this research is not a superior prediction method.

A useful confirmation from this research is that fuzzy set theory is highly flexible [14] and thus, has many other applications including the prediction of time dependent variables [26]. Fuzzy-based methods lend a flexibility to environmental modelling and assessment that crisp methods do not [34,35]. It can be incorporated into other types of modelling schemes such as those that consider connectivity in earth system processes, for example. Connectivity is a subject that describes the degree to which water flow and sediment are related or connected [36,37], and can be used to acquire enough accurate data of past and current conditions to assess the role of landscape processes. When contaminants are costly to analyze, as is the case for many emerging contaminants of concern in coastal regions, Keesstra et al. [36] showed that a connectivity modelling scheme can be used to help define a more cost-effective monitoring and measurement plan.

#### **5. Conclusions**

A TS fuzzy model using the fuzzy c-means clustering algorithm was used to predict lead in soil concentrations under increasingly sparse spatial data conditions. The results from the TS fuzzy model were compared with that of OK and IDW using the same training and validation datasets. The ability of the three methods to predict outlier points within the respective validation sets appeared to be very similar. As with all other performance metrics, the performance of the three methods was difficult to separate and the similarity in model performance may be related to the models themselves. Both IDW and OK generate weights for the surrounding points based on a similar model shape, while the TS fuzzy model utilizes Gaussian shaped membership functions to determine the weighting for each of the cluster's predictions. This research showed that the TS fuzzy model used here did not outperform OK and IDW for the spatial predictions of lead in soil under increasingly sparse data conditions. Future research should consider different types of data, normalizing data prior to fuzzy modeling, other clustering methods, and further optimization of the model parameters.

**Author Contributions:** Conceptualization, R.T., U.T.K. and C.V.; methodology, R.T., U.T.K. and C.V.; validation, R.T.; formal analysis, R.T.; investigation, R.T.; resources, C.V.; data curation, R.T.; writing—original draft preparation, R.T. and F.T.; writing—review and editing, R.T., U.T.K., C.V. and F.T.; supervision, C.V. and U.T.K.; project administration, C.V. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Acknowledgments:** The authors would like to thank the reviewers for their insightful comments in their reviews of this work. The authors would also like to thank the Department of Mechanical Engineering at the University of Victoria for access to software used in this work.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **An In Silico and In Vitro Study for Investigating Estrogenic Endocrine Effects of Emerging Persistent Pollutants Using Primary Hepatocytes from Grey Mullet (***Mugil cephalus***)**

**Paolo Cocci, Gilberto Mosconi and Francesco A. Palermo \***

School of Biosciences and Veterinary Medicine, University of Camerino, Via Gentile III Da Varano, I-62032 Camerino, MC, Italy; paolo.cocci@unicam.it (P.C.); gilberto.mosconi@unicam.it (G.M.) **\*** Correspondence: francesco.palermo@unicam.it

**Abstract:** There is growing concern about the environmentally relevant concentrations of new emerging persistent organic pollutants, such as perfluorinated compounds and pharmaceuticals, which are found to bioaccumulate in aquatic organisms at concentrations suspected to cause reproductive toxicity due to the activation of estrogen receptor (ER) α and β subtypes. Here, we use a combined in silico and in vitro approach to evaluate the impact of perfluorononanoic acid (PFNA) and Enalapril (ENA) on grey mullet (*Mugil cephalus*) hepatic estrogen signaling pathway. ENA had weak agonist activity on ERα while PFNA showed moderate to high agonist binding to both ERs. According to these effects, hepatocytes incubation for 48 h to PFNA resulted in a concentration-dependent upregulation of ER and vitellogenin gene expression profiles, whereas only a small increase was observed in ERα mRNA levels for the highest ENA concentration. These data suggest a structure–activity relationship between hepatic ERs and these emerging pollutants.

**Keywords:** endocrine disruptors; *Mugil cephalus*; PFNA

#### **1. Introduction**

Chemicals interfering with the endocrine system known as endocrine disrupting chemicals (EDCs) are pollutants that typically occur in aquatic environments as a result of municipal wastewater discharge, landfill leachates, and agricultural and urban runoff [1]. EDCs are considered a major cause of aquatic wildlife decline and loss of biodiversity [2]. Aquatic organisms such as fish may experience life-long exposures to EDCs and may bioaccumulate them developing a wide range of hormonal abnormalities [3]. Today, there is growing concern about the environmentally relevant concentrations of new emerging persistent pollutants, such as perfluorinated compounds (PFCs) and pharmaceuticals (e.g., contraceptives and anti-depressants), which are found to bioaccumulate in aquatic food webs at concentrations suspected to perturb neuro-endocrine processes in living organisms including humans [4,5].

PFCs are synthetic chemical compounds that due to high stabilities and low surface tensions are increasingly used in various industrial applications and common consumer products [6]. Among PFCs, the perfluoroalkylated substances (PFAS) have been discovered as global pollutants remaining most persistently in each environmental compartment [7,8]. Although PFAS are considered moderately to highly toxic, some of these (e.g., perfluorooctane sulfonic acid (PFOS) and perfluorooctanoic acid (PFOA)) are suspected endocrine disruptors and have been found to cause adverse health effects, especially reproductive toxicity, in different vertebrate models [9,10]. Similarly, pharmaceuticals that include any chemical product used by individuals or agribusiness for promoting personal and livestock health, have aroused great interest as environmental pollutants for their ecotoxicological potentials [11,12]. These compounds have been detected, unchanged or as metabolites, in wastewater, surface and drinking waters throughout the world [13–15]. Calamari et al. [16]

**Citation:** Cocci, P.; Mosconi, G.; Palermo, F.A. An In Silico and In Vitro Study for Investigating Estrogenic Endocrine Effects of Emerging Persistent Pollutants Using Primary Hepatocytes from Grey Mullet (*Mugil cephalus*). *Environments* **2021**, *8*, 58. https://doi.org/10.3390/ environments8060058

Academic Editor: Sílvia C. Gonçalves

Received: 14 May 2021 Accepted: 16 June 2021 Published: 18 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

have defined pharmaceuticals as pseudo-persistent pollutants due to their continuous introduction into the environment, the biotic and/or abiotic transformation and the ability to exert subtle effects in non-target organisms. In this regard, the estrogenic potential of some pharmaceuticals has attracted great concern especially in the aquatic environment [17,18].

Estrogen-like EDCs (xenoestrogens) have the capability to bind to the estrogen receptors (ERs), mimicking the female steroid hormone, 17β-estradiol (E2), and thus activating intracellular signaling pathways. Activation of the ER-mediated signaling pathway has been extensively studied in several models, particularly fish in which feminization has been considered a direct result of xenoestrogen contamination [19–23]. In this regard, the ER-induced hepatic vitellogenin (Vtg) production is typically used to confirm exposure to estrogenic compounds in male fish [24,25]. Of the different fish organ cells, liver cells are widely used in in vitro primary culture models due to their ability to retain native liver properties including estrogen responsiveness [26–28]. For that reason, in vitro methods using primary cultures of fish hepatocytes represent a fundamental and recommendable alternative to in vivo studies for investigating several toxicologically relevant mechanisms [29,30].

In the present work, we focused our attention on the ability of perfluorononanoic acid (PFNA) and enalapril (ENA) to interfere with estrogen receptor signaling using a combined in silico/in vitro approach. The selected compounds belong to the most frequently detected classes of emerging EDCs in the aquatic environment such as PFAS and pharmaceuticals [31–35]. The Endocrine Disruptome program package was used to predict their potential interference with nuclear ERs in silico. A bioassay that uses primary hepatocytes from the grey mullet (*Mugil cephalus*) was then employed for screening estrogenic potential by assessing classical biomarker responses such as VTG protein and ER isoform mRNA expression. In addition, cytotoxicity using Alamar Blue and 3-(4,5 dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) assays was evaluated after 24 and 48 h exposure.

#### **2. Materials and Methods**

#### *2.1. Endocrine Disruptome Screening Tool*

A molecular docking approach for predicting interactions between PFNA/ENA and estrogen receptor α (ERα) and β (ERβ) ligand binding domains has been performed with Endocrine Disruptome Simulation (EDS) Tool. This web service has already been successfully adopted as a software tool for predicting the endocrine disruption potential of compounds, using well-validated crystal structures of 14 different human nuclear receptors including ER subtypes [36]. The crystal structures of 1A52, 3OLS, 1SJ0, and 1QKN have been chosen as templates on the basis of their sequence identity with fish receptors (higher than 60%). The docking scores reported are a measure of how the contaminants fit within the receptor-binding pocket, taking into account continuum and discreet parameters. According to the threshold calculations sensitivity (SE), it is possible to obtain four broad groups indicating predicted affinity for ER isoforms as follows: "red" (SE < 0.25), high probability; "orange" (0.25 < SE < 0.50) and "yellow" (0.50 < SE < 0.75), medium probability; and "green" (SE > 0.75), low probability of binding [36].

#### *2.2. Hepatocyte Isolation and Primary Cell Culture*

Flathead grey mullet (*Mugil cephalus*) males (95.5 ± 10.9 g initial weight) were provided by professional fishermen during fishing activities. Fish were acclimated for 2 weeks in 2.00 m × 2.00 m × 0.60 m tanks with constant aeration and natural photoperiod at Unità di Ricerca e Didattica of San Benedetto del Tronto (URDIS), University of Camerino in San Benedetto del Tronto (AP, Italy). Water quality parameters were monitored daily showing the following values: pH 8.4 ± 0.2, O2 = 10.3 ± 0.5 mg L−1, and temperature = 20–22 ◦C, salinity 36 ± 2 psu; undetectable level of nitrites and ammonia. Following the acclimation, fish were randomly euthanized using MS-222 within 5 min after capture. Animal manipulation was executed following the procedures established by the Italian law (Leg-

islative Decree 116/1992), the European Communities Council Directive (86/609/EEC and 2010/63/EU) for animal welfare and under the supervision of the authorized investigators. The liver tissue was collected to obtain hepatocytes under a laminar flow hood, according to Cocci et al. [37] and Palermo et al. [38]. Purified hepatocytes were suspended in Leibovitz (L-15) phenol red-free medium, antibiotic-antimycotic solution (100 U/mL) and 10 mM HEPES. The cell density was measured in a Burker Chamber and the viability of hepatocytes was over 90%, as assessed with the Trypan blue exclusion assay. Cells were seeded on 24-well Falcon Primaria culture plates (1 × 106 cells per well) in L-15 phenol red-free medium, antibiotic-antimycotic solution (100 U/mL) and 10 mM HEPES. Cells were cultured for 24 h in an incubator at 23 ◦C before chemical exposure to allow attachment. Then, 50% of the L-15 phenol red-free medium culture was removed, and hepatocytes were exposed to medium containing the vehicle (ethanol, final concentration 0.01%) and 1.0, 0.01, or 0.0001 μM of E2, ENA or PFNA. Hepatocytes were incubated in an incubator at 23 ◦C for 96 h. Media and cells were harvested separately at 0, 24, 48, 72 and 96 h with medium changes every 24 h. Doses of ENA and PFNA were chosen on the basis of environmentally relevant concentrations [33,39–41] and six independent wells were setup for both the control and each concentration of compound. The entire experiment was repeated 3 times.

#### *2.3. MTT Cytotoxicity Assay*

The 3-(4,5-dimethylthiazol-2-yl)-2,5-diphenyl-tetrazoliumbromide (MTT) activity was measured according to Smeets et al. [42] with slight modifications, using the MTT Cell Proliferation and Cytotoxicity Assay Kit assay (Boster Biological Technology, Pleasanton, CA, USA, Catalog # AR1156). After 0, 24, 48, 72, 96 h of treatment described above, incubation medium was removed and replaced with fresh culture medium containing MTT reagent (5 mg/mL MTT diluted in phosphate buffered saline, PBS). After an incubation of 40 min at 23 ◦C, the formazan crystals produced were solubilized by adding 200 μL Formazan solubilization solution. After the complete solubilization, 200 μL of medium was transferred to a 96-well microplate and absorbance values were measured at 570 nm using a microplate reader (BioChrom, Cambridge, UK).

#### *2.4. Alamar Blue Assay*

Cell viability was also quantified using the Alamar Blue™ assay reagent (Thermo Scientific, Waltham, MA, USA) as described by Cocci et al. [23] and following manufacturer's specifications. The incubation medium was removed after 24, 48, 72, 96 h of treatment, replaced with a fresh culture medium containing AB reagent at a final concentration of 10%, and incubated for an additional 1 h. The absorbance was monitored at 570/600 nm using a microplate reader. The cell viability was normalized to that of hepatocytes cultured in the regular media without any of the tested compounds.

#### *2.5. Quantitative Realtime PCR (q-PCR)*

After exposure, the medium was carefully removed, and cells were lysed by adding the TRIzol reagent (Invitrogen Life Technologies, Milan, Italy). Total RNA was isolated according to the manufacturer's specifications. RNA quality and concentration were measured spectrophotometrically at 260/280 nm, and purity was confirmed by electrophoresis through 1% agarose gels stained with SafeView Classic (abm). The cDNA was synthesized from 1.5 μg of total RNA in 20 μL using the 5X All-In-One RT MasterMix (with AccuRT Genomic DNA Removal Kit) according to manufacturer's instructions (abm). SYBR greenbased real-time PCR was used to evaluate expression profiles of ERα, ERβ, VTG target genes. 18s rRNA was selected as appropriate reference gene [28,43,44]. All the primer sequences are reported in Table 1 and were provided from Ribecco et al. [45], Vieira et al. [46], Cabas et al. [47], and Perez-Sanchez et al. [48]. The reaction included 10 μL of 2X BlasTaqTM qPCR MasterMix (abm), 0.5 μL each of forward and reverse primers (10 μM), 2 μL of cDNA template, and nuclease-free H2O to a final volume of 20 μL. The expression of individual

gene targets was analyzed using the ABI 7300 Real-Time PCR software. Thermo-cycling for all reactions was for 3 min at 95 ◦C, followed by 40 cycle of 15 s at 95 ◦C, and 60 s at 60 ◦C. Dissociation curve analysis revealed that a single peak was generated during the reaction demonstrating the production of a single product. Each amplified fragment was then compared with that obtained from amplification of *Sparus aurata* cDNA and verified with agarose gel electrophoresis (for details see Figure S1 in Supplementary Material). The efficiency of qPCR primer sets was reported in Table 1. Results were calculated using the 2−ΔΔCt method and reported as fold change corrected for 18s rRNA and with respect to vehicle levels. Values are the mean ± SD of three independent experiments.


**Table 1.** List of primers used in this study.

#### *2.6. Enzyme-Linked Immunosorbent Assay (ELISA)*

VTG concentrations in the culture medium of *Mugil cephalus* hepatocytes were determined using an ELISA method previously published [25]. Cell culture media were diluted 1:8 as reported for routinely diluted media samples by Navas and Segner [49]. All samples were analyzed in triplicate. Absorbance was recorded at 492 nm using a microplate reader (Biochrom).

#### *2.7. Statistical Analysis*

Data were assessed with Graphpad prism v6.01 software (GraphPad Software Inc., San Diego, CA, USA) and expressed as mean ± standard error of the mean (SEM). Statistical analysis was performed using ANOVA (one-way analysis of variance) followed by Bonferroni's multiple comparison test. Differences with *p* < 0.05 were considered statistically significant.

#### **3. Results and Discussion**

The prediction results obtained with the EDS model for ENA and PFNA are given in Table 2.

ENA is a drug of the class of angiotensin-converting enzyme inhibitors (ACEI) that is mainly used in the treatment of arterial hypertension. Several studies indicate a beneficial interaction between ACEI and estrogens which in turn are involved in reducing ACE mRNA concentrations [50]. Zilberman et al. [51] showed that chronic exposure to ENA significantly up-regulated ERα and β protein expression in rats. To date, however, there is no experimental evidence that ACEI can bind directly to ERs. Thus, it is not surprising that, according to EDS simulation, ENA presents moderate binding affinity against the agonist-active conformation of the ERα and low affinity for both conformations of the ERβ, respectively.

PFNA is one of the three main long chain PFCs, primarily used as an emulsifier for producing fluoropolymers, that can be found at high concentrations in the environment [52]. PFNA has been detected in various waters and animal tissues worldwide suggesting high bioaccumulation potential in the food chain [53,54]. In mammalian studies, an interference with gonadal development in neonatal mice was observed after gestational exposure to PFNA [55]. In addition, an increase in estrogenic activity was reported for exposure to different PFCs, including PFNA, in in vivo studies using fish as models [56,57]. The

results of in silico analysis obtained for PFNA allow us to support these effects because the simulation tool predicts its agonist activity on the ERs, showing a high probability of binding on ERβ and a slightly below probability on ERα. Only the agonist activity on the ERα from different species was previously described in the literature [56]. This latter paper reported that PFNA docked into the LBD region of ERα working as in vitro weak binders and activators. In the present study, PFNA docking scores for ERβ were more favorable than those for the ERα. In addition, an antagonist activity on the ERβ was also predicted but any evidence of this potential activity has been found in the literature.


**Table 2.** Prediction affinities for 17β-Estradiol (E2), Enalapril (ENA) and perfluorononanoic acid (PFNA). +/− indicate the crystal structures of estrogen receptors (ER) isoforms in complex with their respective agonist (+) or antagonist ligands (−).

> In order to determine whether a structure–activity relationship between ERs and tested compounds was evident, primary cultures of grey mullet hepatocytes were used to investigate the impact of PFNA and ENA exposure on ERα/β and VTG expression. In the first stage of our study, we examined whether the tested chemical compounds affected cellular viability by using common in vitro cytotoxicity assays such as Alamar Blue and MTT. The cell viability was expressed as metabolic activity, displaying 100% viability in the media of both negative (EtOH) and positive control (E2) at each time point (Figure 1A).

> In contrast, a significant inhibition of metabolic activity was obtained for both ENA and PFNA. We found that cell viability decreased in hepatocytes after exposure to the highest doses of ENA for 72 h (84.8% and 81.1% of solvent control at 0.01 and 1 μM ENA, respectively) and 96 h (89.3% and 65.0% at 0.01 and 1 μM ENA, respectively) (Figure 1B). Similarly, exposure of hepatocytes to 0.01–1 μM PFNA induced a significant change in viability at both 72 h (cell viability 86.7% and 67.5%, respectively) and 96 h (cell viability 71.7% and 53.9%, respectively) after the start of treatment (Figure 1C). Complete cell death was confirmed by microscopic examination of cells exposed to the tested chemicals. The observed effects on cell viability were further investigated using the MTT assay (Figure 2). ENA and PFNA performed similarly in both assays, showing the most significant decrease in hepatocyte viability (from 71.7% to 53.9%) exposed to the highest concentrations (0.01 and 1 μM) for 96 h (Figure 2B,C). However, the MTT assay failed to detect cytotoxicity for PFNA at 0.01 μM after 72 h exposure, thus proving to be slightly less sensitive than the Alamar Blue assay.

> To our knowledge, only one study has shown that ENA elicits specific cytotoxic effects in primary cultures of hepatocytes through the involvement of a glutathione-dependent detoxification pathway [58]. This effect was observed at concentrations ranging from 0.5 to 2 mM in an in vitro rat cell model. The potential cytotoxic and antiproliferative effects of ENA were also found at concentration- and time-dependent manners in human HL60 acute promyelocytic leukaemia cells [59]. Interestingly, viability of ENA-treated HL60 cells was observed to drop by about 20% after exposure to 3 μM for 48 h. On the other hand, the

effect of PFNA on hepatocyte viability was previously demonstrated in humans, using the WST-1 assay.

**Figure 1.** Alamar Blue cell viability of *Mugil cephalus* hepatocytes following exposure to E2 (**A**), ENA (**B**), PFNA (**C**) for up to 96 h. The dot line represents cell viability measured in the solvent control (assigned a survival of 100%). Values are given as mean ± SEM of three independent experiments and expressed as % relative to the solvent control. "\*" indicates significant differences between control and treated groups (*p* < 0.05).

PFNA was found to be more cytotoxic than PFOA and PFOS, causing a larger decrease in cell viability upon exposure for 6–72 h to a concentration range of 200– 400 μM [60]. According to the obtained results, treatment incubation for 24 or 48 h was applied in the further studies. However, the exposure duration and sampling time of 24 h for all molecular endpoints were found to be inappropriate to obtain a clear concentration response of the model compounds (Figure 3A–C). As expected, 48 h exposure to positive control (E2) produced a dose-dependent induced expression of ERs and VTG compared with control hepatocytes (Figure 3D). In contrast, a partial concentration-response curve (0.01–1 μM) for expression of all molecular endpoints was obtained for the model compound PFNA after 48 h of exposure (Figure 3E), whereas only a small increase was observed in ERα mRNA levels for 1 μM ENA (Figure 3F).

**Figure 2.** MTT assay of *Mugil cephalus* hepatocytes following exposure to E2 (**A**), ENA (**B**), PFNA (**C**) for up to 96 h. The dot line represents cell viability measured in the solvent control (assigned a survival of 100%). Values are given as mean ± SEM of three independent experiments in % relative to the solvent control. "\*" indicates significant differences between control and treated groups (*p* < 0.05).

We further discuss the mRNA expression of VTG at the level of protein concentrations in the medium used for the primary cultures of grey mullet hepatocytes. At 48 h after treatment, E2 caused a significant dose-dependent increase in VTG synthesis at any of the tested concentrations relative to control cultures (Figure 4). Both PFNA and ENA also increased the VTG synthesis, but to a lesser extent. Indeed, a significant increase in VTG levels occurred following treatment with the highest doses (0.01 and 1 μM) of PFNA. In contrast, VTG up-regulation was only induced in response to exposure at 1 μM ENA.

**Figure 3.** ERα, ERβ and VTG mRNA expression profile (fold change from vehicle) in *Mugil cephalus* hepatocytes exposed to different concentrations (0.0001–1 μM) of E2 (**A**,**D**), PFNA (**B**,**E**), ENA (**C**,**F**) for 24 and 48 h. Values are mean ± SEM of three independent experiments. "\*" indicates significant differences between control and treated groups (*p* < 0.05).

ER-mediated production of VTG is likewise by far the most used biomarker of xenoestrogen exposure in oviparous species [61–65]. VTG protein and gene expression has been shown to be up-regulated by various environmental pollutants in a number of in vitro–in vivo studies using fish models [66–68]. Interestingly, most of these works have found that VTG induction is accompanied by a clear increase in hepatic ER expression, mainly ERα. This latter is indeed considered the ER subtype with a major role in mediating VTG gene induction. However, there is evidence that ERβ subtype may have a functional role in the up-regulation of ERα enhancing hepatic VTG induction in response to E2 or xenoestrogen stimulation [69,70]. Our results on PFNA are in agreement with those of Benninghoff et al. [56] who, in a recent study, found both in vitro and in vivo

weak estrogenic activities of this PFC using a similar range of concentrations. Indeed, according to the relative binding affinity (RBA) values obtained in vitro, dietary PFNA also induced a consistent in vivo VTG induction in trout [56]. Collectively, our findings provide clear confirmation of data reported so far regarding the ability of PFNA to act as weak environmental xenoestrogen. PFNA has frequently been detected in surface waters at concentrations in the order of ng/L showing particularly high levels (up to a max of 100 ng/L) in the recycling sites due to the recycling activities of electrical and electronic waste [40,71]. Thus, the prominence of PFNA in water from this typology of sites suggests the need of a careful monitoring of this potential xenoestrogen in order to reduce its ecological impact in aquatic ecosystems.

Similarly, the occurrence of ENA in the environment can be related to incomplete removal of this drug from wastewater treatment plants (WWTPs). ENA was detected in wastewater influent at concentrations ranging from 35 to 1400 ng/L, while in the wastewater effluent this range was reduced to 0.85–290 ng/L [39]. These data demonstrate that total or high removal of this drug can be achieved in all WWTPs, thus suggesting a substantially lower accumulation rate. Given also our results about the weak estrogenic potential, one might predict that ENA has a mild impact on reproductive functions of aquatic vertebrates.

#### **4. Conclusions**

In summary, an in vitro hepatocyte bioassay was used to characterize estrogenic responses of gray mullet to PFNA and ENA as representative compounds of two classes of emerging pollutants. According to the environmental occurrence of these chemicals, the results would indicate potential adverse impacts especially on reproductive health. The observed effects are likely to be mediated through direct actions of these compounds on hepatic ERs suggesting a structure–activity relationship between ER and these emerging pollutants. While the extent of PFNA estrogenic potential is substantially supported by literature data, further studies such as in vivo investigations need to obtain more information on the estrogenic activity of ENA, especially to check its effectiveness following long-term accumulation.

**Supplementary Materials:** The following are available online at https://www.mdpi.com/article/10 .3390/environments8060058/s1, Figure S1. Results of PCR amplification using template DNA from *Mugil cephalus* or *Sparus aurata*.

**Author Contributions:** Conceptualization, P.C. and F.A.P.; methodology, P.C., G.M., F.A.P.; formal analysis, F.A.P.; investigation, P.C., G.M.; resources, F.A.P.; writing—original draft preparation, P.C., F.A.P.; writing—review and editing, F.A.P.; supervision, F.A.P. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** The data presented in this study are available in the article/supplementary material, further inquiries can be directed to the corresponding author.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


## *Article* **Assessment of the Levels of Pollution and of Their Risks by Radioactivity and Trace Metals on Marine Edible Fish and Crustaceans at the Bay of Bengal (Chattogram, Bangladesh)**

**Krishna Prasad Biswas 1, Shahadat Hossain 2,3, Nipa Deb 2, A.K.M. Saiful Islam Bhuian 2,\*, Sílvia C. Gonçalves 4,\*, Shahadat Hossain <sup>2</sup> and Mohammad Belal Hossen <sup>1</sup>**


**Citation:** Biswas, K.P.; Hossain, S.; Deb, N.; Bhuian, A.K.M.S.I.; Gonçalves, S.C.; Hossain, S.; Hossen, M.B. Assessment of the Levels of Pollution and of Their Risks by Radioactivity and Trace Metals on Marine Edible Fish and Crustaceans at the Bay of Bengal (Chattogram, Bangladesh). *Environments* **2021**, *8*, 13. https://doi.org/10.3390/


Academic Editors: Claude Fortin and Vernon Hodge

Received: 10 November 2020 Accepted: 1 February 2021 Published: 11 February 2021

environments8020013

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

**Abstract:** Marine environmental pollution is a longstanding global problem and has a particular impact on the Bay of Bengal. Effluent from different sources directly enters rivers of the region and eventually flows into the Bay of Bengal. This effluent may contain radioactive materials and trace metals and pose a serious threat to the coastal environment, in addition to aquatic ecosystems. Using gamma spectrometry and atomic absorption spectrometry, a comprehensive study was carried out on the radioactivity (226Ra, 232Th, 40K, and 137Cs) and trace metal (Cd, Pb, Zn, Cu, Ni, Fe, Mn, and Cr) concentrations, respectively, in fish and crustacean species collected from the coastal belt of the Bay of Bengal (Chattogram, Bangladesh). The analysis showed a noticeable increment in the levels of different radioactive pollutants in the marine samples, although the consumption of the studied fish and crustacean species should be considered safe for human health. Anthropogenic radionuclide ( 137Cs) was not detected in any sample. Furthermore, the metal concentrations of a small number of trace elements (Pb, Cd, Cr) were found to be higher in most of the samples, which indicates aquatic fauna are subject to pollution. The estimated daily intake (EDI), target hazard quotient (THQ), hazard index (HI), and target cancer risk (TR) were calculated and compared with the permissible safety limits. It was found that consuming the seafood from the Bay of Bengal may cause adverse health impacts if consumption and/or means of pollution are not controlled.

**Keywords:** radioactive materials; trace metals; bioaccumulation; marine fish; crustaceans; marine environmental pollution; Bay of Bengal

#### **1. Introduction**

The marine environment is one of the most important sources of life on Earth and performs a number of key environmental functions for the lives and livelihood of humans and other organisms. Marine environmental pollution is a longstanding global problem. Marine pollution affects the Bay of Bengal because of the vulnerability of its aquatic habitats to such pollution following the industrial revolution in the 19th century. Most aquatic ecosystems can cope with a certain degree of pollution, but severe pollution is reflected in changes in the fauna and flora of their communities [1]. Developing countries, including Bangladesh, are most affected by this human-made problem. Bangladesh is considered to be one of the most suitable countries for the Blue Economy following the recent establishment of a vast maritime boundary in the Bay of Bengal. The result of two verdicts on maritime boundaries, also involving Myanmar and India, allows Bangladesh to exclusively

exercise its sovereign rights over 118,813 sq km of water. The affected area extends up to 12 nautical miles into territorial sea and includes a further Exclusive Economic Zone (EEZ) of 200 nautical miles [2]. However, agriculture runoff, untreated sewage, and industrial pollution are reportedly being discharged directly into rivers, eventually flowing into the Bay of Bengal. This effluent may contain radioactive materials and trace metals and pose a particularly serious threat to the coastal environment and the aquatic ecosystems [3,4].

The main causes of radioactivity in the marine environment are seabed movement originating from underwater volcanic activity and undersea earthquakes, natural processes of weathering, and mineral recycling of terrestrial rocks [5]. Different types of geological materials and minerals, such as ores and igneous rocks, which often contain large concentrations of natural radioactive elements, contribute to the transfer of radionuclides to water through leaching action [3]. The level of naturally occurring radioactive elements in the marine environment has gradually increased due to human activities, such as mining and processing of ores, production of natural oil and gas, and combustion of fossil fuels in coal-fired power plants. Moreover, anthropogenic contributions are made by underwater nuclear device tests, post-nuclear disposal of industrial and radioactive waste, recycling of spent nuclear fuel, and accidents, including leaks from nuclear power plants [4–7]. Among the radionuclides, uranium, radium, and radon are soluble in seawater, whereas thorium is almost completely insoluble. These can dissolve in seawater and then attach to sediment on the seabed and suspended plankton in the seawater. These dissolved radionuclides and plankton, in the long run, contaminate marine organisms, including fish, crustaceans, and several types of shellfish [8].

In addition, trace metals are released into the environment from different sources, such as transportation, industrial activities, fossil fuels, agriculture, urbanization, and other human activities [6–9]. The release of large quantities of trace metals into nature has resulted in several environmental problems due to their non-biodegradability and persistence. Marine life can have a considerable aptitude for bioaccumulation and biosorption of trace metals [10–12]. Under certain environmental conditions, such as consolidated effects of biotic and abiotic factors, such as plants, animals, and microbes; the amount of sunlight in the ecosystem; the amount of oxygen and nutrients dissolved in the water; proximity to land; depth; and temperature; these trace metals may accumulate to a toxic concentration and cause ecological damage by threatening the health of aquatic and terrestrial organisms, including humans [1,13,14]. The principal pathway leading to human exposure to radionuclides and trace metals is the consumption of seafood. A large variety of seafood is eaten. Fish and crustacean species, because of their important dual role as both an integral element of marine food webs and a commercial food product for humans, are frequently selected as indicator organisms in studies related to pollution or the bioaccumulation processes of various noxious substances [15]. The bioaccumulation of trace metals and isotopic contaminants by the tissues and organs of marine organisms has been studied globally and has led to the adoption of the bioindicator concept for environmental quality assessments [5–8,16,17].

Research on the bioaccumulation of radionuclides and trace metals in the most commonly consumed fish and crustacean species in the Bay of Bengal region is limited. Given the importance of such knowledge, this study was carried out to determine the radioactivity and trace metal levels in six marine fish and four crustacean species of the northern coastal belt of the Bay of Bengal, Bangladesh. The bioaccumulation levels of radioactive materials (226Ra, 232Th, 40K, and 137Cs) and trace metals (Cd, Pb, Zn, Cu, Ni, Fe, Mn, and Cr) were estimated by gamma spectrometry and atomic absorption spectrometry, respectively. Different radiological and health-hazard parameters were also calculated to assess the marine environmental quality in this region concerning these pollutants, in addition to the health risks resulting from the consumption of the studied seafood.

#### **2. Materials and Methods**

#### *2.1. Sample Collection and Preparation*

The studied field locations were openly accessible, and the organisms collected did not belong to any protected species, so permissions were not required at the time of collecting samples from the studied locations. The investigated fish and crustaceans were obtained from four fish/seafood markets of the northern coastal belt of the Bay of Bengal, Chattogram, Bangladesh (Figure 1), during the rainy (July) and autumn (September and October) seasons of 2017. Bangladesh is a country of six seasons (summer, rainy, fall, autumn, winter, and spring). The fish investigated in this study are mostly available in rainy and autumn seasons. Hence, we obtained our samples in these two seasons. The fish samples were obtained from fishing boats operating in the corresponding locations, while the crustacean samples were obtained from a local sea food market in Chattogram city. A total of 20 samples (10 samples from each season) were chosen for investigation. The studied organisms were six fish and four crustacean species, locally available and commercially important, i.e., *Tenualosa ilisha* (Ilish), *Harpodon neherues* (Lotia), *Sillaginopsis panijus* (Sundora Baila), *Sardinella longiceps* (Colombo), *Trichiurus lepturus* (Churi), and *Konosirus punctatus* (Shad) concerning the fish species, and *Scylla serrata* (mud crab), *Portunus sanguinolentus* (three-spot swimming crab), *Hemigrapsus takanoi* (Asian shore crab) and *Penaeus semisulcatus* (cat tiger shrimp) concerning the crustaceans.

**Figure 1.** Location of the study area on the (**a**) world map, (**b**) Bangladesh map, and (**c**) Bay of Bengal map.

The details of the collected samples are given in Table 1. All the samples were labelled, stored in ice, and, on the same day, transported to the laboratory and washed with clean water and distilled water. Afterwards, the samples were dried in filter paper, packed in polyethylene bags, and stored in a refrigerator of the Laboratory of Atomic Energy Centre, Chattogram, to avoid degradation, spoiling, contamination, or any other decomposition until further treatment and analysis were done. The collected fish and crustacean samples were measured for the total wet weight (WW) and recorded beforehand. For the analysis of radionuclides and trace metal concentrations, all the animals were used in their entirety, since the analysis was intended to quantify the radioactivity and the metal concentrations of each organism.


**Table 1.** Marine organisms collected and analyzed during rainy and autumn seasons of 2017 at the Bay of Bengal (Bangladesh).

Note: (L-1)—Fishery Ghat, Chattogram; (L-2)—Gohira Fishery Ghat, Anwara, Chattogram; (L-3)—Ganga Bari Fish Market, Chattogram; (L-4)—Hazari Goli Seafood Market, Chattogram; R—rainy season and A—autumn season.

> Each sample was sun-dried at least one week to remove the extra water and subsequently oven-dried at about 70 ◦C to obtain a constant weight. The samples were then re-weighed to determine the dry/wet ratio. Finally, the dried samples were ground, sieved for homogeneity, and stored in clean and dry uncontaminated empty cylindrical plastic containers of uniform size (2.8 cm × 8.0 cm), sealed with wide vinyl adhesive tapes around their screw necks to air tighten for succeeding uses.

> For the trace metal determination, a part of the powdered samples was kept aside in clean and dry sealable plastic bags and kept in an air-tight glass jar to save the samples from any types of chemical reactions. Five grams of the homogenized powder were taken from each specimen and placed in a 250 mL digestion beaker. A digestion mixture containing 6.0 mL of high purity nitric acid (70%) and 2.0 mL of hydrochloric acid (37%) were added and heated at 60 ◦C for half an hour. After cooling down for 15 min, 4.0 mL of hydrogen peroxide (30%) was added to each beaker and heated until a clear solution was obtained and the volume decreased into half of its original volume following the standard operating procedure described elsewhere [18]. The digested portions were filtered through Whatman filter paper (No. 42) and diluted to a final volume of 50 mL using de-ionized water.

#### *2.2. Methods for Determining Radioactivity and Trace Metals*

#### 2.2.1. Gamma Spectrometry Analysis

For the natural and artificial radioactivity measurement of each sample, about 100 g of dried fish/crustacean samples were weighed and placed in individual plastic containers. The samples were kept 4 weeks before the analysis in airtight conditions to allow secular equilibrium between thorium and radium and their short-lived progenies [14].

The activity concentrations of 226Ra, 232Th, 40K, and 137Cs in the samples were determined using an NATS GCD-40 190 p-type HPGe (Baltic Scientific Instruments, Riga, Latvia), gamma-ray spectrometer (63.3 mm crystal diameter, 61.3 mm thickness, +2200 V

operating bias voltage). The low background detector of relative efficiency of 43.1% and an energy resolution of 1.74 keV full width at half maximum (FWHM) at the 1332 keV peak of 60Co was enclosed in a cylindrical lead shield. The counting system was connected to an 8 k multi-channel analyzer (MCA 527, GBS Elektronik GmbH) with associated electronics for data acquisition of photo-peak areas. A Spectral Line GP© was used to analyze the gammaray counts received from the samples. Energy calibration and the absolute photo-peak efficiency evaluation were performed using standard reference material (Code: 8501-EG-SVE, Eckert and Ziegler Analytics), diluted with a multi-nuclide gamma-ray source (241Am, 109Cd, 57Co, 139Ce, 113Sn, 85Sr, 137Cs, 88Y, 60Co) having homogeneously distributed activity, and maintaining the same geometry and density as the plastic containers containing the samples. To minimize the statistical counting error, the samples were counted for a period of 30,000 s. An empty container was also counted under the same conditions to determine the background counts. To calculate the specific activities, the background counts for the same counting condition were deducted from the counts of each sample to obtain the net counts. For spectrum analyses, the single transition gamma-ray line 1460.822 keV was used to determine the activity concentrations of 40K. The gamma-ray photo-peaks of 295.221 keV and 351.922 keV from 214Pb, and 609.320 keV, 1120.310 keV and 1764.551 keV from 214Bi were used to determine the activity concentrations of 226Ra. The activity concentrations of 232Th were determined using the net counts under the 238.630 keV and 300.087 keV photopeaks from 212Pb, 911.205 keV and 968.970 keV photo peaks from 228Ac, and 583.190 keV and 2614.533 keV from 208Tl. For the evaluation of 226Ra and 232Th activity, a weighted mean approach was applied using the aforementioned gamma lines [19,20]. Much care was taken to prevent contamination during the investigation.

#### 2.2.2. Determination of Trace Metals

The fish and crustacean samples were analyzed for eight trace elements, Cd, Pb, Zn, Cu, Ni, Fe, Mn, and Cr, using an atomic absorption spectrophotometer (HITACHI Z-2000) [21]. Four standard solutions of different known concentrations were prepared, and the elemental concentrations in the unknown samples were determined by extrapolation from the calibration curve. The digested samples were directly aspirated into the flame (air–acetylene fuel mixture). The concentration corresponding to the absorption in the digest was determined by using the absorption mode. The minimum detection limits of the analyzer for the investigated trace metals in ppm are 0.002 (Cd), 0.010 (Cu), 0.005 (Zn), 0.050 (Pb), 0.010 (Mn), 0.020 (Fe), 0.020 (Ni), and 0.020 (Cr) [22]. The limits of detection (LODs) for all the elements analyzed in the samples were calculated as the blank signal plus three times its standard deviation, whereas the limit of quantification (LOQ) was calculated as ten times the standard deviation of the blank signal, following [23,24]. The recovery estimation was performed by the laboratory during method development and validation processes using standard reference materials (SRMs) and certified reference materials (CRMs). However, this time we used secondary reference material made from commercial standard solutions of each element. At each step of the digestion processes, acid blanks (laboratory blank) were prepared to ensure that the samples and chemicals used were not contaminated. Each set of digestion had its acid blank and was corrected by using its blank. The measurement of each sample was taken three times. All the experimental values were reported as mean value ± standard deviation (SD). Throughout the analysis, all the trace metal standards were prepared and run to check the precision of the instrument. The quality control and quality assurance protocol set by the U.S. Environmental Protection Agency for metal analysis were used. The quality assurance testing relied on the blank controls, and the yield of the chemical procedure was used.

#### *2.3. Radiation Dose and Risk Assessment*

2.3.1. Radiation Dose Assessment

Activity Concentration

The activity concentrations (Bq kg<sup>−</sup>1) of the natural and anthropogenic radionuclides in the measured samples were evaluated by the following equation [25]:

$$\mathbf{A} \left( \text{Bq } \text{kg}^{-1} \right) = \frac{\text{CPS} \times 1000}{\varepsilon\_{\text{Y}} \mathbf{I}\_{\text{Y}} \mathbf{M}} \tag{1}$$

where CPS = net count per second (i.e., CPS for sample—CPS for background), ε<sup>γ</sup> = detector efficiency of the specific γ-ray, I<sup>γ</sup> = intensity of the gamma-ray, and, M = mass of the sample in grams.

The errors of the measurements were expressed in terms of the standard deviation of the ±1σ level.

#### Total Effective Dose

Annual effective dose is a useful concept that enables the radiation doses from different radionuclides and from different types of sources of radioactivity to be added. Estimation of the radiation-induced health effects associated with the intake of radionuclides in the body is proportional to the total dose delivered by the radionuclides while resident in the various organs. Radiation doses ingested are obtained by measuring radionuclide activity in foodstuffs (Bqkg<sup>−</sup>1) and multiplying these by the masses of food consumed over a certain period (kgd−<sup>1</sup> or kgy<sup>−</sup>1). A dose conversion factor (SvBq−1) can then be applied to give an estimate of the ingestion dose. Thus, the ingested dose is given by [26,27]:

$$\text{annual effective dose} \left(\text{Svy}^{-1}\right) \text{ = concentration} \left(\text{Bq}, \text{kg}^{-1}\right) \times \text{annual intake} \left(\text{kg}^{-1}\right) \times \text{DCF} \left(\text{SvBq}^{-1}\right) \tag{2}$$

where DCF is the standard dose conversion factor, which is equal to 0.2800 μSvBq−<sup>1</sup> for 226Ra, 0.2300 μSvBq−<sup>1</sup> for 232Th, and 0.0062 μSvBq−<sup>1</sup> for 40K [25].

Therefore, the total effective dose (Svy−1) via ingestion is calculated by the following formula:

$$\text{total annual effective dose} = \left\{ \mathbf{C}\_{\text{R}} (^{226}\text{Ra}) \times \text{I}\_{\text{F}} \times \text{E}\_{\text{D}} \right\} + \left\{ \mathbf{C}\_{\text{R}} (^{222}\text{Th}) \times \text{I}\_{\text{F}} \times \text{E}\_{\text{D}} \right\} + \left\{ \mathbf{C}\_{\text{R}} (^{40}\text{K}) \times \text{I}\_{\text{F}} \times \text{E}\_{\text{D}} \right\} \tag{3}$$

where CR is the concentration of radionuclides in ingested fish/crustaceans (Bqkg<sup>−</sup>1), IF is the annual intake (kgy−1) of fish/crustaceans containing radionuclides, and ED is the ingestion dose conversion factor for radionuclides (SvBq−1). The intake rates for Bangladeshi consumers were taken from the "Year Book of Fisheries Statistics of Bangladesh 2018" [28].

2.3.2. Radiation Risk Assessment

Internal Hazard Index (Hin)

Radon and its short-lived descendants are hazardous to the respiratory organs. The internal hazard index (Hin) is used to quantify the internal exposure to radon and its progenies, which is given by the equation:

$$\mathrm{H}\_{\mathrm{in}} = \frac{\mathrm{C}\_{\mathrm{Ra}}}{185} + \frac{\mathrm{C}\_{\mathrm{Th}}}{259} + \frac{\mathrm{C}\_{\mathrm{K}}}{4810} \tag{4}$$

where CRa, CTh, and CK are the concentration (Bqkg<sup>−</sup>1) of Ra, Th, and K, respectively.

The values of the internal hazard index (Hin) must be less than unity for the radiation hazard to be negligible [29,30].

Excess Lifetime Cancer Risk (ELCR)

Nowadays, cancer is called a life-threatening disease, and the percentage of this disease increases all over the world, including in Bangladesh, due to various reasons. One of the reasons is the effect of radiation on the biological cells, which contributes to a greater extent to increasing cancer incidence. An effort was made to assess the excess lifetime cancer risk due to the ingestion of marine fish and crustaceans by the procedure proposed by the United States Environmental Protection Agency (USEPA) [31]. The following equation was used to calculate the excess lifetime cancer risk [12,32]:

$$\text{ELCR} = \text{A}\_{\text{ir}} \times \text{A}\_{\text{ls}} \times \text{R}\_{\text{c}} \tag{5}$$

where *ELCR*, *Air*, *Als*, and *Rc* are the excess lifetime cancer risk, the annual intake of radionuclide (Bq), the average lifespan 72 yr. (life expectancy of a male and a female are 70.6 years and 73.5 years, respectively in Bangladesh [33]), and the mortality cancer risk coefficient (Bq–1), respectively. The values of mortality cancer risk coefficients are 9.56 × <sup>10</sup>−<sup>9</sup> Bq–1 for 226Ra, 2.45 × <sup>10</sup>−<sup>9</sup> Bq–1 for 232Th, and 5.89 × <sup>10</sup>−<sup>10</sup> Bq–1 for 40K [12,31]. The acceptable ELCR limit is 10−<sup>3</sup> for radiological risk in general [32,34].

#### *2.4. Health Risk Assessment of Trace Metals*

#### 2.4.1. Estimated Daily Intake (EDI)

EDI is measured by the following equation in (mgkg−<sup>1</sup> body weight per day) [35]:

$$\text{EDI} = \frac{\text{E}\_{\text{F}} \times \text{E}\_{\text{D}} \times \text{F}\_{\text{IR}} \times \text{C}\_{\text{f}} \times \text{C}\_{\text{M}}}{\text{W}\_{\text{AB}} \times \text{AT}\_{\text{n}}} \times 10^{-3} \tag{6}$$

where EF is the exposure frequency (365 days per year), ED=72 years is the exposure duration, FIR is the ingestion rate, which is taken as 62.58 g per person per day [28], Cf is the conversion factor (Cf = 0.208) to convert fresh weight to dry weight considering 79% of moisture content of the fish fillet [35,36], CM is the metal concentration in fish fillet (mgkg−<sup>1</sup> dry weight basis), WAB is the average body weight (60 kg) of Bangladeshi adult people [37], and ATn (equal to EF × ED) is the average exposure time for non-carcinogens [35].

Several organizations such as WHO, FAO, etc., provided guidelines on the intake of metals by a human. The maximum tolerable daily intake (MTDI), provisional tolerable daily intake (PTDI), and the acceptable daily intake (ADI) are used to describe the safe levels of intake for several toxins, including toxic metals. The EDI of trace metals measures the amount of metal intake by an adult (60 kg) by ingestion of fish or crustaceans per day.

#### 2.4.2. Target Hazard Quotient (THQ)

The THQ is non-carcinogenic risk and is dimensionless. In this study, the noncarcinogenic health risks associated with the consumption of fish and crustacean species were assessed based on the target hazard quotients (THQs), and their calculations were done using the standard assumption for an integrated USEPA risk analysis as follows [38]:

$$\text{THQ} = \frac{\text{E}\_{\text{F}} \times \text{E}\_{\text{D}} \times \text{F}\_{\text{IR}} \times \text{C}\_{\text{f}} \times \text{C}\_{\text{M}}}{\text{W}\_{\text{AB}} \times \text{AT}\_{\text{n}} \times \text{RfD}} \times 10^{-3} \tag{7}$$

EF, ED, FIR, Cf, CM, WAB, and ATn are explained in the earlier section. RfD is the reference dose of individual metal (mg/kg/day) (Table 2). The RfD represents an estimate of the daily exposure to which the human population may be continually exposed over a lifetime without a significant risk of deleterious effects. If the THQ is less than 1, the exposed population is unlikely to experience obvious adverse effects. If the THQ is equal to or higher than 1, then there is a potential health risk [39], and related interventions and protective measurements should be taken.


**Table 2.** Reference dose (RfD) and carcinogenic slope factor (CSFo), oral.

2.4.3. Hazard Index (HI)

To estimate the overall potential health risk related with more than one metal, THQ of every metal is summed up and referred to as hazard index (HI). The HI can be calculated by the sum of the target hazard quotients (THQs) of each metal:

HI = THQMn + THQFe + THQCu + THQZn + THQPb + THQCd + THQCr + THQNi (8)

#### 2.4.4. Target Cancer Risk

For carcinogens, the target cancer risk (lifetime cancer risk) is estimated as the incremental probability of an individual to develop cancer over a lifetime exposure to that potential carcinogen (i.e., incremental or excess individual lifetime cancer risk) [38]. Acceptable risk levels for carcinogens range from 10−<sup>4</sup> (risk of developing cancer over a human lifetime is 1 in 10,000) to 10−<sup>6</sup> (risk of developing cancer over a human lifetime is 1 in 1,000,000). The equation used for estimating the target cancer risk, which is dimensionless, is as follows [38]:

$$\text{TR} = \frac{\text{E}\_{\text{F}} \times \text{E}\_{\text{D}} \times \text{F}\_{\text{IR}} \times \text{C}\_{\text{f}} \times \text{C}\_{\text{M}} \times \text{CSF}\_{\text{o}}}{W\_{\text{AB}} \times \text{AT}\_{\text{C}}} \times 10^{-3} \tag{9}$$

Here, CSFo is the oral carcinogenic slope factor (mg/kg BW/day)−1, and ATc (equal to EF × ED) is the average exposure time for carcinogens. Since Mn, Fe, Cu, and Zn do not cause any carcinogenic effects, their CSFo have yet not been established in USEPA 2012. Thus, TR values for the intake of Pb, Cd, Cr, and Ni were calculated to show the carcinogenic risk using oral carcinogenic slope factors (CSFo) of these toxic elements given in Table 2.

#### **3. Results and Discussion**

Chattogram is the second largest city in Bangladesh, where the main seaport of Bangladesh is located on the bank of the Karnaphuli river. The Chattogram seaport handles ninety per cent of Bangladesh's export–import trade and has been used by India, Nepal, and Bhutan to join for transshipment. There are thousands of industries including a textile mill, cement factory, tannery, oil refinery, dichlorodiphenyltrichloroethane (DDT) plant, chemical factory, fertilizer factory, paper mills, power plant, dry-dock, paint factory, rayon mills, etc., situated near to the bank of the Karnaphuli river, among which about 200 are identified as pollution causing units continuously discharging unlawfully a huge amount of pollutants that finally enter into the Bay of Bengal [42]. The wastes coming from the municipal sewage system through many canals of Chattogram city is another potential factor of the pollution problem of this bay. Moreover, illegal and/or accidental discharges of grease, fish oil, bilge, garbage, etc., from the merchant and fishing vessels are causing pollution to the Bay of Bengal [43].

#### *3.1. Radiation Dose and Risk Assessment by Radioactivity Analysis*

Marine fish and crustaceans are potential bio-indicators when they accumulate the target radionuclides from surrounding waters [15,44]. Monitoring radionuclide levels in fish and crustaceans is of great importance due to their significant contribution to the natural radiation dose received by human beings consuming them [45].

Among all the fish and crustacean samples, 226Ra was detected only in *Harpodon neherues*, and its value was 5 ± 2 Bqkg<sup>−</sup>1. There was no anthropogenic radionuclide (137Cs) in any sample. In all the studied fish samples, mean activity concentrations of 232Th and 40K ranged from 7 ± 1 to 190 ± 10 and 210 ± 50 to 360 ± 40 Bqkg<sup>−</sup>1, respectively. Similarly, in crustacean samples, the activity concentrations of 232Th and 40K ranged from 5.0 ± <sup>2</sup> to 53 ± 10 and 130 ± 40 to 240 ± 70 Bqkg−1, respectively. Further, note that the activity concentrations of 40K were greater than those of the other radionuclides for all samples (Table 3). Mean values of the activity concentrations of 226Ra, 232Th, and 40K in all the fish and crustacean samples were 5 ± 2, 67 ± 9, and 250 ± 50 Bqkg−<sup>1</sup> , respectively.

**Table 3.** Mean activity concentration (Bqkg−1) and excess lifetime cancer risk (ELCR) due to the consumption of natural radionuclide from marine fish and crustacean samples collected at the Bay of Bengal.


Note: BDL—below detection limit.

The comparative data of the present study with the previously reported studies of different radioactive elements in fish and crustacean samples of the Bay of Bengal region are reported in Table 4. It was reported in 2012 that the activity concentrations of 226Ra, 232Th, and 40K in sea-fish samples of the Black Sea region of Turkey were 0.06 ± 0.01 to 0.96 ± 0.36, 0.12 ± 0.04 to 1.03 ± 0.15, and 35.04 ± 0.24 to 127.4 ± 2.3 Bqkg<sup>−</sup>1, respectively [14].

In all these cases, the results were very much lower than those of the present observations. Conversely, the mean activity concentration of 226Ra, 232Th, and 40K in marine fish of the Straits of Malacca was reported as 4.05 ± 0.48 to 7.83 ± 0.78, 1.93 ± 0.24 to 6.21 ± 0.53, and <sup>288</sup> ± 15 to 399 ± 20 Bqkg<sup>−</sup>1, respectively, in 2015 [12], being considerably comparable to the present study except the concentration of 232Th. The obtained activity concentration of 232Th and 40K in all the fish and crustacean samples were 67 ± 9 and 250 ± 50 Bqkg<sup>−</sup>1, respectively, which are almost similar to the activity concentration obtained in crab samples from river Odeomi, Ogun State, at the Southwest of Nigeria [46]. It should be noted that if we liken the present study to previous studies for fish [47,48] at the same region, there is a noticeable increment in the levels of different radioactive elements in those marine organisms of the Bay of Bengal region (shown in Table 4). This increment might be due to the wide variety of activities (e.g., housing, tourism, power generation plants, petroleum, steel, shipbuilding, chemical, pharmaceutical, textile, vegetable oil refineries, glass manufacturing industries, etc.) around the Bay of Bengal region, which have been diversified over the years. Multifunctionality and diversification of industries, fisheries, and agriculture are probably linked with these changes in the marine environment. The absence of the anthropogenic radionuclide ( 137Cs) indicates that there was not any effect on the marine organisms analyzed due to the post-nuclear disposal of industrial and radioactive waste, underwater nuclear device tests, accidents including leaks from nuclear power plants and from recycling of spent nuclear fuel, etc. As Bangladesh is establishing its first nuclear power plant (NPP) at Rooppur, Pabna, this baseline data can help to further monitor these activities.

**Table 4.** Comparative data of the radioactivity (Bqkg<sup>−</sup>1) in marine fish and crustacean samples of the Bay of Bengal, Bangladesh.


Note: BDL—below detection limit.

The ELCR values of 226Ra, 232Th, and 40K in all the fish and crustaceans are shown in Table 3. It is to be mentioned here that ELCR for 226R was below detection level (BDL), except for in *Harpodon neherues*. However, all of the obtained values (shown in Table 3) were lower than the acceptable ELCR limit of 10−<sup>3</sup> for radiological risk in general, indicating no health hazard to the consumers [32,34].

The estimated annual effective doses due to intake of 226Ra, 232Th, and 40K from the consumption of seafood (fish and crabs) from the Bay of Bengal are presented in Table 5, and the comparison between the mean annual effective dose and its world average values are shown in Figure 2. The average annual effective dose for 226Ra and 40K were found to be within UNSCEAR acceptable limits. Though the average annual effective dose for 232Th was found to be three times greater than the UNSCEAR acceptable limits, the total values of annual effective dose were within the acceptable limits, indicating no threat to the consumers. Maximum and minimum values of the total internal hazard index were 0.79 ± 0.06 and 0.04 ± 0.01, respectively, and all the values were less than 1, showing that there is no health hazard to the consumers.


**Table 5.** Annual effective dose, total effective dose and total internal hazard index for different radionuclides in marine fish and crustacean samples collected at the Bay of Bengal (Bangladesh) during rainy and autumn seasons of 2017.

**Figure 2.** Comparison between the estimated annual effective doses obtained in the present study with the world average values [32].

#### *3.2. Pollution Level and Health Risk Assessment of Trace Metal Analysis*

#### 3.2.1. Trace Metal Concentration

The intake, bioassimilation, and subsequent bioaccumulation of trace metals in marine organisms are significantly affected by the dissimilar aquatic geochemistry of the trace metals as well as the diverse feeding habits and living methods of fish and crustaceans. Among marine animals, fish are possibly one of the most mobile and capable of travelling long distances. Compared to other types of organisms, fish are also on a high trophic level in the food chain; hence, their diet is probably the most diverse of the species studied here. Crustaceans are benthic organisms that usually live on or near the sea floor and are capable

of travelling some distance. They often feed on organic debris but also on small crustaceans and fish on or near the sea floor. Although the trace metal concentrations in the various species of marine fish and crustaceans analyzed in the present study were in a wide range of variations, these aquatic organisms also showed metal accumulation patterns that were noteworthy (see Table 6).

**Table 6.** Mean metal concentration (mgkg–1 DW) in different marine fish and crustacean samples collected at the Bay of Bengal (Bangladesh) during rainy and autumn seasons of 2017.


Note: BDL—below detection limit.

The trace metal concentrations found in individual fish muscles for both rainy and autumn seasons varied for Pb, 7.5–0.12, Cu: 4.3–1.3; Zn, 170–45; Mn, 24–5.5; Cd, 2.1–0.10; Ni, 13–0.28; Cr, 2.6–0.16; and Fe, 1300–11 mgkg−<sup>1</sup> dry weight. In addition, in individual crustacean samples of both rainy and autumn seasons, these values varied for Pb, 1.9–0.25; Cu, 100–39; Zn, 150–79; Mn, 640–15; Cd, 3.4–0.09; Ni, 2.3–0.31; Cr, 3.3–0.63; and Fe, 600–9.4 mgkg−<sup>1</sup> dry weight.

Most of the chemical elements present in fish and crustaceans are essential for biota at low concentrations; however, high concentrations of these elements can become toxic for them. Similarly, some metals (e.g., Fe, Mn, Ni, Cu, and Zn) are necessary for proper metabolic reactions in the human body, and other elements (e.g., Cd, Cr, and Pb) are of unknown benefits and may become toxic in the cases of chronic exposure to humans [50]. The obtained trace metal concentrations in the ten marine fish and crustacean samples analyzed in the present study in the rainy season followed the order of Fe > Mn > Zn >

Cu > Ni > Cr > Cd > Pb and in the autumn season followed the order of Fe > Mn > Zn > Cu >Cr > Pb > Cd > Ni. However, in both seasons, the average concentration followed the order of Fe > Zn > Mn > Cu > Ni > Pb > Cr > Cd in fish samples and of Fe > Mn > Zn > Cu > Cr > Cd > Pb > Ni in crustacean samples. This demonstrated that marine organisms in different groups had different accumulation mechanisms for trace metals. Additionally, the different concentrations of chemical elements in the marine environment and the age of the marine organisms analyzed may also be responsible for these trends.

The maximum allowed concentration of Pb in fish is 0.21 mgkg−<sup>1</sup> DW set by JECFA [51], whereas in our study we found that the fish *Sillaginopsis panijus*, *Sardinella longiceps*, *Trichiurus lepturus*, and all the crustacean species contained Pb concentrations exceeding the limit. According to FAO/WHO, the acceptable limit of Cd for human consumption is 0.20 mgkg−<sup>1</sup> [52]. Among the fish species, *Sardinella longiceps*, *Konosirus punctatus*, and among the crustacean species *Portunus sanguinolentus*, *Scylla serrata*, and *Penaeus semisulcatus* exceeded this limit. In addition, the permissible limit for Cr is 0.20 mgkg−<sup>1</sup> DW set by the National Academic Press, Washington DC, 1989 [51], whereas all the studied samples exceeded this limit. These three trace metals (Cd, Cr, and Pb) are considered to be toxic and all of them were in an excessive content in some of the studied marine organisms, demonstrating that the aquatic fauna was contaminated by these metals. Other essential trace elements (e.g., Cu, Zn, Mn, Ni, and Fe) were also observed in excessive contents in some of the studied samples. This indicates the overall pollution status of the studied area due to industrial activities, fossil fuels, agricultural run-off, ship-breaking activities, and other human activities around the coastline.

A comparison between available fish species metal concentration data and the present study is shown in Table 7. It is evident that the range of concentration of metal elements in other parts of the world are lower than that of the presently studied data [52–56]. On the other hand, the obtained data is comparable to similar studies in the neighboring regions [57,58]. If we liken the present study with similar kinds of studies conducted in the Bay of Bengal in the past, the range of metal element concentrations are also comparable [59–61], though it indicates a gradual increment in metal accumulation by the fish of those regions. It seems that day by day, the rate of environmental disturbances and the pollution levels increase due to different reasons, such as the establishment of different industries beside the rivers, which is gradually increasing and also affects the marine environment progressively, contaminating the marine organisms in the long run. Presently, it is an issue of concern for the local populations, but in the near future it might be the threat for the surrounding environment of all organisms, including humans.

**Table 7.** Trace metal concentrations (mgkg−<sup>1</sup> DW) in various species of fish from different areas of the world.



**Table 7.** *Cont.*

Note: BDL—below detection limit; ND—not detected; PS—present study.

#### 3.2.2. Estimated Daily Intake (EDI)

EDI and mean EDI values of respective trace elements (Cd, Pb, Zn, Cu, Ni, Fe, Mn, and Cr) for the fish and crustacean samples of both seasons are shown in Table 8. The corresponding Figure 3 shows the comparative results of estimated daily intake (EDI) values with the recommended daily intake (RDI) values. As the values of EDI < RDI for both fish and crustaceans, there is no risk from the consumption of the studied seafood.

**Table 8.** The values of estimated daily intake (EDI), target hazard quotients (THQ), and hazard index (HI) obtained for the marine fish and crustacean samples collected at the Bay of Bengal (Bangladesh) during rainy and autumn seasons of 2017.


**Figure 3.** Comparison between estimated daily intake (EDI) values and recommended daily intake (RDI) values obtained in the present study.

#### 3.2.3. Target Hazard Quotients (THQ)

The dimensionless target hazard quotient (THQ) is the indicator of non-carcinogenic health risks associated with the consumption of food (fish and crustaceans). The mean THQ values for the trace metals Pb, Cu, Zn, Mn, Cd, Ni, Cr and Fe are given in Table 8. After analysis of all the fish and crustacean samples, the obtained THQ values for different metals were less than 1, which indicates that the exposed population is unlikely to experience obvious adverse effects [39].

#### 3.2.4. Hazard Index (HI)

The value of the hazard index, which is obtained by the summing of the target hazard quotients of each metal, was used to assess the overall potential health risk posed by more than one metal. The obtained hazard index value was 8.82 × <sup>10</sup>−4, which is less than 1, indicating that there are no health hazards to the local consumers (Table 8).

3.2.5. Target Cancer Risk (TR)

The obtained values of target cancer risk for Pb, Cd, Ni, and Cr due to exposure from the consumption of the targeted six fish and four crustacean species analyzed in the present study are shown in Figure 4. Generally, the values of TR lower than 10−<sup>6</sup> are considered as negligible, while those above 10−<sup>4</sup> are considered to be unacceptable, and those in between 10−<sup>6</sup> and 10−<sup>4</sup> are considered as an acceptable range [37]. The present study reveals that only TR for Pb were below the benchmark and those for Cd, Ni, and Cr were above the benchmark, indicating that the fish and the crustaceans were becoming polluted. This also increased the risk of chronic cancer due to exposure of Cd, Ni, and Cr through fish and crustacean consumption.

**Figure 4.** Targeted cancer risk (TR) values for associated toxic element consumption from the studied seafood (fish and crustaceans) at the Bay of Bengal (Bangladesh) during rainy and autumn seasons of 2017.

#### **4. Conclusions**

The radioactivity and the trace metal concentrations were studied in certain commercially important marine biota (fish and crustaceans) from the Bay of Bengal, Bangladesh, in two different seasons (rainy and autumn). Through the analysis of different radiological parameters, the present study revealed an elevated activity concentration compared to the acceptable limits in some of the samples. Comparing the present results with the reported results of the different regions of the world, including similar studies at the Bay of Bengal, it can be concluded that the studied region possesses more radioactive elements than the other marine regions. It was also observed that the marine radioactivity of the studied region is gradually increasing, day by day. The results reflect the contribution of technologically enhanced naturally occurring radioactive material (TENORM) pollutants, largely expected to be a result of power generation plants, petroleum, steel, shipbuilding, chemical, pharmaceutical, textile, vegetable oil refineries, glass manufacturing industries, etc. However, the present study indicates that radionuclide intake from the consumption of Bay of Bengal fish and crustaceans still poses an insignificant threat to public health.

Conversely, the metal concentrations of a few trace elements are higher than the acceptable limits in most of the samples, although the studied fish and crustaceans are still safe for human consumption. However, the target cancer risk (TR> 10<sup>−</sup>4) due to exposure to cadmium, nickel, and chromium indicates that consumer risk persists.

As a whole, the present study indicates an increase in the pollution by radioactivity and trace metals in the marine environment of the Bay of Bengal, mostly derived from an increase in human activities in the region. Therefore, consuming the seafood from the studied region has the potential to cause adverse health impacts if not controlled, and at the same time, all the stakeholders should take proper initiatives to prevent environmental pollution of the Bay of Bengal. The results reported here can be used as baseline data, though the sample size is not large enough. Similar studies should be planned and carried

out periodically with higher sample sizes to monitor any changes in marine pollution of the Bay of Bengal in the future. Statistical analyses should also be employed in future studies for better understanding of the scenario.

**Author Contributions:** Conceptualization, K.P.B., A.K.M.S.I.B., S.H.2 and M.B.H.; methodology, S.H.1, N.D., A.K.M.S.I.B. and S.H.2; software, K.P.B., S.H.1 and N.D.; validation, A.K.M.S.I.B., S.C.G., S.H.2 and M.B.H.; formal analysis, K.P.B., S.H.1 and N.D.; investigation, K.P.B., S.H.1, N.D. and A.K.M.S.I.B.; resources, A.K.M.S.I.B. and S.H.2; data curation, K.P.B., S.H.1 and N.D.; writing original draft preparation, K.P.B. and A.K.M.S.I.B.; writing—review and editing, S.H.1, A.K.M.S.I.B. and S.C.G.; visualization, A.K.M.S.I.B. and S.H.2; supervision, A.K.M.S.I.B., S.H.2 and M.B.H.; and project administration, S.H.2 and M.B.H. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research received no external funding.

**Data Availability Statement:** The data presented in this study are contained within the article.

**Acknowledgments:** The authors are thankful to the authority of the Atomic Energy Centre, Chattogram, Bangladesh Atomic Energy Commission, for providing support and technical co-operation to analyze the specimens using the HPGe and AAS. The authors are also delighted to express their gratefulness to the Ministry of Shipping, Peoples' Republic of Bangladesh, for their utmost support. Thanks also goes to Masuda Begum Sampa, Kyushu Institute of Technology, Japan for the fruitful discussion.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


### *Article* **A Global Meta-Analysis for Estimating Local Ecosystem Service Value Functions**

**Luiz Magalhães Filho 1,2,\*, Peter Roebeling 1,3, Maria Isabel Bastos 1, Waldecy Rodrigues <sup>4</sup> and Giulia Ometto <sup>5</sup>**


**Abstract:** The Meta-analysis has increasingly been used to synthesize the ecosystem services literature, with some testing of the use of such analyses to transfer benefits. These are typically based on local primary studies. However, meta-analyses associated with ecosystem services are a potentially powerful tool for transferring benefits, especially for environmental assets for which no primary studies are available. In this study we use the Ecosystem Service Valuation Database (ESVD), which brings together 1350 value estimates from more than 320 studies around the world, to estimate meta-regression functions for Provisioning, Regulating and maintenance, and Cultural ecosystem services across 12 biomes. We tested the reliability of these meta-regression functions and found that even using variables with high explanatory power, transfer errors could still be large. We show that meta-analytic transfer performs better than simple value transfer and, in addition, that local meta-analytical transfer (i.e., based on local explanatory variable values) provides more reliable estimates than global meta-analytical transfer (i.e., based on mean global explanatory variable values). Thus, we conclude that when taking into account the characteristics of the study area under analysis, including explanatory variables such as income, population density, and protection status, we can determine the value of ecosystem services with greater accuracy.

**Keywords:** ecosystem services; benefit transfer; meta-analysis; meta-regression function

#### **1. Introduction**

Jean-Baptiste Say poses the idea of nature's services as costless, free gifts of nature as follows: "the wind which turns our mills, and even the heat of the sun, work for us; but happily no one has yet been able to say, the wind and the sun are mine, and the service which they render must be paid for" [1]. However, currently, it is possible to observe that the overuse or misuse of some natural resources poses direct impacts on society. In the face of this problem came the concept of ecosystem services (ES), defined as the benefits that humans obtain from the natural environment and from properly-functioning ecosystems. Hence, several authors [2–8] argue that the sustainable management of natural resources requires correct valuation of the ecosystem defining their services to the society.

Ecosystem services support human life every day and contribute to human well-being in many ways, which are hard to define in a single notion. Hence, the Millennium Ecosystem Assessment [9] and the Common International Classification of Ecosystem Services [10] differentiate between the following ecosystem services: (a) Provisioning services (such as the supply of food via fishery production, fuel, wood, energy resources, and natural products); (b) Regulating and maintenance services (such as shoreline protection, nutrient

**Citation:** Magalhães Filho, L.; Roebeling, P.; Bastos, M.I.; Rodrigues, W.; Ometto, G. A Global Meta-Analysis for Estimating Local Ecosystem Service Value Functions. *Environments* **2021**, *8*, 76. https://doi.org/10.3390/ environments8080076

Academic Editor: Sílvia C. Gonçalves

Received: 18 June 2021 Accepted: 4 August 2021 Published: 9 August 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

regulation, carbon sequestration, detoxification of polluted waters, and waste disposal); and (c) Cultural services (such as tourism and recreation).

Ecosystems have great importance across many dimensions (ecological, socio-cultural, and economic) [5]. Thus, expressing the value of ecosystem services in monetary units (i.e., ecosystem service values; ESV), can prove to be of utmost importance to help raise consciousness and convey the (relative) importance of ecosystems and biodiversity to decision-makers. Indeed, monetized valuation pushes for more efficient use of limited resources and helps to select where protection and regeneration are economically more important and can be delivered at least cost [11,12]. It can also assist in determining "a fair compensation" to be paid for a loss of ES in liability regimes [5].

Historically, in the late 1990s and early 2000s, the concept of ES slowly found its way into the policy arena, e.g., through the "Ecosystem Approach" and the Global Biodiversity Assessment. In 2005, the concept of ES gained wider interest after the publication of the Millennium Ecosystem Assessment by the United Nations for policymakers [4,9]. ES are also entering the consciousness of mainstream media and business, namely through the World Business Council for Sustainable Development that has actively supported and developed this concept [13]. Many projects and groups are currently working toward better understanding, modeling, valuing, and managing ES and natural capital [4].

An increasing number of papers seeking the valuation of ES have been published over the last decades. Assessments have been conducted at local [14–16], national [17–19], continental [20,21], and global [4–6] scales. In the same way, databases compiling data from these primary valuation studies were created to aggregate information and facilitate public debate and policy action. Some examples of such databases include the Economic Valuation Reference Inventory [22], and the Ecosystem Service Valuation Database (ESVD; [23]).

Since the early 1990s, several researchers have investigated the applicability and the precision of benefit transfer. However, these past studies were primarily concerned with traditional methods of benefit transfer (in particular value transfer), replacing values directly from the study site to the policy site without amendments [24]. However, in the late 1990s meta-analysis started to be used, with multivariate regression being investigated for use in benefit transfer [25].

The meta-analysis (MA) is a technique that uses statistical models (meta-regressions) to summarize and evaluate previous research results. In benefit transfer, meta-regression results may be used qualitatively, to corroborate new primary results, or to transfer estimated values [26]. Meta-regression in benefit transfer summarizes the weight of the evidence and characterizes the degree of uncertainty about quality-adjusted ecosystem values. In meta-regression, the value estimates from primary valuation studies are thereby treated as individual observations [27]. Meta-regression also extends the range of primary valuation studies by allowing the estimation of values for services and functions that are constant within each primary valuation study but vary across different valuation studies [28].

Meta-analyses have been performed for specific ecosystem services, biomes, and locations. For example, Van Houtven et al. [15] assessed the cultural value of surface water quality in the United States, using 131 willingness-to-pay (WTP) estimates from 18 studies. Similarly, Hjerpe et al. [29] synthesized 127 WTP estimates from 22 different studies that provided estimations for preservation, forest restoration, and freshwater restoration also in the United States. Ghermandi et al. [30] performed a meta-analysis to determine the values of goods and services provided by wetland ecosystems, using 418 value observations derived from 170 valuation studies and 186 wetland sites worldwide. Finally, Hynes et al. [31] performed a marine recreational meta-analysis estimation, using 311 distinct value observations from 96 primary valuation studies. Nevertheless, there are no studies with a broader analysis, that estimate global meta-regression functions for Provisioning, Regulating and maintenance and Cultural ecosystem services across biomes and continents. In addition, testing the reliability of estimated meta-regression functions is relatively rare, (e.g., [32,33]). One of the main challenges is developing equations for ES that capture the local/regional characteristics of the biome and provide reliable value estimates.

Hence, the objective of this study is to estimate meta-regression functions for 3 different types of ecosystem services able to determine the ecosystem service value for 12 different types of biomes, with the possibility of these estimates being applied at the global scale. In this study, we provide the results of a meta-analysis based on the primary value estimates from the Ecosystem Service Valuation Database [23] for 3 ecosystem services (Provisioning; Regulating and maintenance; Cultural), provided by 12 main land covers (Coastal systems; Coastal wetlands; Coral reefs; Cultivated areas; Desert; Fresh water; Grasslands; Inland wetlands; Open Ocean; Temperate/Boreal forests; Tropical forests; Woodlands). In addition, complementary explanatory variables from the World Bank Data [34] and FAOSTAT [35] were gathered. Based on this review and meta-analysis, we aim to provide recommendations for future research that may enhance the use of ecosystem service valuation for policy analysis.

The remainder of this paper is structured as follows. The "Materials and Methods" section details the MA application and use in ES studies, the theoretical specification and validation method and, finally, the ESVD database and other variables used to build the models. In turn, in the "Results" section, we expose and analyze the functional forms of the models for the three ecosystem services, present the application of the models, and discusses the results. Finally, in the "Conclusions and recommendations" section, concluding remarks and observations are presented.

#### **2. Materials and Methods**

#### *2.1. Literature on Meta-Analysis*

Benefit transfer (BT) is an economic valuation tool, with the goal to adapt value estimates from past research to assess the value of a similar, but separate, change in a different resource [36]. Technically, BT uses valuation estimates from other areas (study sites) and applies them to a similar location (policy site) [3]. It is a technique that relies on primary studies and, therefore, allows for the reduction of field research constraints, both in terms of time and infrastructure. However, it can lead to over/underestimated values while the accuracy of an ESV estimate is determined by the quality of the reference studies used. Thus, peer-reviewed empirical studies from similar biophysical and socioeconomic contexts are preferred over any other type of data source [32].

BT is useful when the estimation of the economic service value cannot be obtained due to time and/or budget constraints and to, therefore, make the best possible use of the existing literature in order to evaluate the economic importance of a natural area [19]. This is possible by adopting and applying estimates from existing studies that best suit the new context, using one or more of the following BT methods: (i) benefit estimate or value transfer, which is the extrapolation of estimates from one site to another (i.e., values are directly transposed from the study site to the policy site without amendments), (ii) benefit function transfer, which is the transfer of economic functions between the sites (i.e., coefficients are used to determine the policy site values), (iii) meta-analysis, which combines the findings of independent studies related to the research topic as to summarize the body of evidence relating to a particular issue, and (iv) preference calibration, which uses existing benefit estimates derived from different methodologies and combines them to develop a theoretically consistent estimate for policy site values [37].

The meta-analysis (MA) technique can help reduce deviations in value estimates [26]. This technique was first put forward as a research synthesis method and has since been developed and applied in many fields of research, other than the area of environmental economics [38,39]. It is widely recognized that the large and increasing literature on economic valuation of ES and environmental impacts has become difficult to interpret and that there is a need for research synthesis, especially in statistical MA, to aggregate information and insights [27,40,41].

MA is by definition a quantitative analysis of statistical summary indicators reported in a series of similar empirical studies. It is a commonly used method for compiling and analyzing the data from studies towards the creation of a value function. The method synthesizes the results of multiple studies that examine the same phenomenon, through the identification of a common effect, which is then "explained" using regression techniques in a meta-regression model [40]. In the realm of environmental resource valuation, MA is commonly used in benefit transfer endeavors due to its usefulness in incorporating a structural utility framework with less strictly economic information [27,42].

#### *2.2. Specification of the Meta-Regression*

Based on consumer rationality and reasonableness, the microeconomic consumer theory is explained by two different approaches: the indifference curve approach and the utility function approach [43]. Indifference curves represent all combinations of goods and services that provide the same level of satisfaction to an individual (i.e., the same level of global utility). Implicit in an indifference curve is the marginal rate of substitution, which expresses the maximum amount of a good that one is willing to give up in exchange for one additional unit of another good, at the same level of satisfaction [43]. Utility functions represent the degree of profitability or satisfaction that we get from using goods and services, related to a measure of satisfaction relative to an economic agent. The analysis of its variation allows for explaining the behavior that results from the decisions taken by each agent to increase his/her satisfaction.

Any meta-analytic benefit transfer (MA-BT) must be based on the ecosystem service valuation theory and the utility functions theory (see Equation (1)), specific to microeconomics [42]. The general form of an MA-BT underlying the utility function is given by [24]:

$$\mathcal{U}\_{i} = f(P\_{i\prime}, \mathcal{Y}\_{i\prime} \mid \mathcal{Q}\_{i\prime} \mathcal{Q} l\_{i\prime} \mathcal{S} sub\_{i\prime} \ H\_{i\prime} l\_{i}) \tag{1}$$

where *Ui* is the utility (satisfaction) obtained by individual *i*, *Pi* is the general price level faced by individual *i*, *Yi* is the individual revenue, *Qi* is the quantity of ES available to individual *i*, *Qli* is the global quality of ES available to individual *i*, *Subi* represents the substitutes for Q available to individual *i*, *Hi* refers to other non-income attributes of individual *i*, and *Ii* is the information available to individual *i*.

Resorting to this microeconomic theoretic, we organize the MA-BT utility theory into three axes: the "strong structural utility theoretic (SSUT) approach", the "weak structural utility theoretic (WSUT) approach" and the "non-structural utility theoretic (NSUT) approach" (of which they only endorse the first two) [42].

Following the microeconomics reasoning, we assume that MA-BT is based on the utility function (see Equation (2)) and opt for analyzing the WSUT. Under the WSUT, each individual may choose between two alternative environmental options—ceteris paribus, a damaged ecosystem (*Q*0) and a restored ecosystem (*Q*1), which will assure an equilibrium situation (the maximum utility) [7,42], represented by:

$$\mathcal{U}\_{i}\begin{pmatrix} P\_{i\prime} \ \mathbf{Y}\_{i\prime} \ \mathbf{Q}\_{0} \end{pmatrix} = \mathcal{U}\_{i}\begin{pmatrix} P\_{i\prime} \ \mathbf{Y}\_{i\prime} - ESV\_{\prime} \ \mathbf{Q}\_{1} \end{pmatrix} \tag{2}$$

where *Ui* is the utility obtained by individual *i*, *Pi* is general price level faced by individual *i*, *Yi* is individual revenue, *Q*<sup>0</sup> quality/quantity of ES available to individual *i* in the absence of any payment, ESV is ecosystem service value paid by individual *i*, and *Q*<sup>1</sup> is the quality/quantity of ecosystem available to individual *i* after having paid for these ES.

Microeconomics utility theory will hold if both sides of this parity are equal. That is, an individual will stay at the same indifference curve if he/she gets the same level of satisfaction by consuming *Q*<sup>0</sup> with no payment or by consuming *Q*<sup>1</sup> and paying ESV in exchange. That is, the ESV the individual is willing to give up must be counterbalanced by an increase in *Q*. Thus, *Q*<sup>1</sup> > *Q*<sup>0</sup> after the amount has been spent.

In this study, we adopt the WSUT approach, where variables are added in the bidfunction (assumed to be derived from some unidentified utility function), but keeping the

flexibility to incorporate other explanatory variables into the ESV model, such as studysite characteristics, local price levels or local individual income [7]. This is the approach used in most previous MA-BT studies [7,31,44]. Our general theoretical model will focus on estimating the ESV (see Equation (3)), as a function of various explanatory variables according to the general form of the underlying conditional indirect utility function:

$$ESV = f(B\_{l'}, SQ\_{l'}, \mathcal{C}\_r QQ\_{r'}, I\_r, P\_r) \tag{3}$$

where, *Bl* is the biome and *SQl* the quality status for the location under analysis (*l*), *C* is the continent where the study area is located, and *QQr* is the quality/quantity of protected areas, *Ir* is the income and *Pr* is the population density in the region (*r*) where the study area is located.

The meta-modeling approach has several advantages for BT as compared to other methods (such as value transfer or function transfer). Different from those, which are based on single studies, MA resorts to information from a collection of studies and, thus, provides more rigorous measures of central tendencies that are sensitive to the distributions of underlying study values [24].

#### *2.3. Validity and Reliability of a Meta-Analytic Benefit Transfer*

The validity and reliability of the MA-BT can be assessed by applying the concept of transfer error (*TE*), defined as [32]:

$$TE = \frac{|ESV\_P - ESV\_B|}{ESV\_B} \tag{4}$$

where *ESVP* is the predicted value from the study site (s) and *ESVB* is the base value ("benchmark") at the policy site. The *TE* is often used as a validity measure of the acceptability of meta-models. Traditionally, validity requires that the values, or the value functions generated from the study site, be statistically identical to those estimated at the policy site [8]. The main objective is to find a target value of TE = 0, confirming that the estimated values from the MA-BT values are similar to those arising from the database.

There is no agreement on maximum TE levels for BT being reliable for different policy applications. The TE analysis is not supposed to judge which levels should be considered acceptable, or even conduct traditional statistical tests of BT validity. Instead, it remains a measure of reliability, especially if TE estimates are compared across metamodel specifications and restrictions, and between alternative ways of conducting BT based on the same data [7].

Therefore, we perform the following comparisons between the estimates from the meta-model and the original observations from the database:


#### *2.4. Background and Data*

MA in environmental valuation is, generally, based on brief statistics and analytical conclusions taking a group of studies as data. Therefore, MA estimates can reduce the time spent to acquire data—both in the case of older studies and unpublished work (where data may not be available) and current studies (where authors may be slow to disclose data). However, even within the same methodology, combining primary data is not always possible due to conflicting data structures and different estimation procedures [42]. This might limit the MA studies' representativeness.

A solution to this problem is the use of specialized ESV databases, which offer a wide range of detailed information about the studies taken into account, beyond the results found in the assessment. These databases give information on other factors crucial for the delimitation of a MA model, such as: the year of the study, protection status, location, type of environment, and method. In this analysis, we use the Ecosystem Service Valuation Database (ESVD, [23]), one of the biggest databases containing real values for a range of ES and biomes where the value estimates are systematized in monetary units (€/ha/year).

The ESVD was built to process and analyze the monetary estimates of ES values from different biomes in a way that is easily used by various end-users, worldwide. Composed by 267 studies and 1310 value estimates, the ESVD links various types of information from different studies with the value estimates and case study sites. These value estimates are organized by biome, ecosystem service, and country. The main biomes are "Coastal System" (*CSys*), "Coastal Wetland" (*CWet*), "Coral Reef" (*CoRf*), "Cultivated Area" (*CuAr*), "Desert" (*Dser*), "Grasslands" (*Gras*), "Inland Wetland" (*InWt*), "Marine" (*Mari*), "Temperate or Boreal forests" (*TeFo*), "Tropical forests" (*TrFo*), "Fresh water" (*FrWa*) and "Woodland" (*Wood*). The ecosystem services are Provisioning; Regulating & maintenance and; Cultural services, divided into 14 types of services (see in Figure 1). Finally, a total of 80 countries are included, 217 values from Africa; 352 values from Asia, 208 values from Europe, 180 values from Latin America and the Caribbean; 122 values from North America, 116 from Oceania, and 114 from the whole world.

**Figure 1.** Ecosystem service values division [45].

Initial criteria for selecting studies from the general ESVD database were: (1) original nature of the case study data (i.e., not based on value transfer or total ecosystem value); (2) the provision of a complete set of information, including the study site location, surface area and the scale of the study (i.e., not based on a "world" scale location); (3) clear characterization of valuation methodologies used (i.e., not unknown valuation methods); (4) clear mentioning of the surface area for which the ecosystem service valuation study is applied (so that estimates of monetary values per hectare can be obtained); and (5) ES or sub-service monetary value directly linked to a specific biome/ ecosystem and unit (i.e., not per person or household). Besides information on the location of each case study, the ESVD includes information on protection status and the size of the research area, enabling for the verification of whether more estimates about the same case study location are available from other sources or publications. Together with supplementary variables, coming from complementary socio-economic databases that are added to

ESVD, these variables allow for further socio-economic interpretation of the monetary output values.

In order to relate an estimate of an ecosystem service to the socio-economic context of a case study site, two additional variables were included in the country table—namely the Gross national income (*GNI*) per capita (based on purchasing power parity in current international prices) and the average Population density (*PDen*; people per square kilometer). This information was obtained from the World Bank Data, which provides world development indicators by country [34]. Collected values were obtained for the years in which the studies were carried out.

Regarding protection status, many of the data points in the ESVD pertain to case studies in protected areas. This information allows the assessment of the influence of the protection status on ES value, testing whether protection excludes the user's access to the site and consequently to the services generated or, alternatively, whether it allows for ecosystem conservation and subsequent appreciation of the services. Protection status is classified into 3 categories: Fully protected (*FProt*), Partially protected (*PProt*) or Not protected (*NProt*). Other complementary variables collected from the World Bank Data, used to verify the study-site protection status, were: Terrestrial Protected Areas (*TProt*; the percentage of protected land by country) and Marine Protected Areas (*MProt*; the percentage of protected territorial waters). From the Food and Agriculture Organization (FAO) statistical database [35], information on the land use characteristics was collected. Namely, the percentage of forest area (*FPer*) and the percentage of agricultural land (*APer*), which helped to understand land use and occupation characteristics with emphasis on agricultural activities and state of preservation/conservation of nature.

For each biome in the ESVD, 14 ecosystem services were identified and classified into the 3 main classes: Provisioning, Regulating and maintenance, and Cultural services (see Figure 1). This classification constitutes an important step in the linkage between ES and human well-being and will be used as a basis to perform MA-BT for ecosystem valuation. Provisioning services (*ESVProv*) are mainly composed of food provision, water provision (including regulation of water flows and water purification), fuels and fibers provision, and genetic resources provision [45]. This is an ES highly valued by humans, because of the direct impact on our day-to-day life. Regulating and maintenance services (*ESVReg&Main*) help maintaining air, climate, and water quality, moderating extreme events, maintaining soil quality, and preventing erosion. Albeit usually invisible and taken for granted, they are important for human well-being and the conservation of plants and animals [45]. Finally, cultural services (*ESVCult*) entail non-material benefits that people obtain from an ecosystem, such as aesthetic inspiration, recreation, and tourism as well as spiritual experience related to a natural environment [45].

All monetary values in the ESVD values are converted into a common reference unit, specifically 2015 'International' €/ha/year, using the Purchasing Power Parity (PPP) units expressed in Euros [34,45].

#### **3. Results**

#### *3.1. Data Summary*

Based on the above-mentioned criteria, the total number of monetary value estimates included in our sample amount to 636 observations. In this study, ES value functions are estimated for Provisioning, Regulating and maintenance, and Cultural services (see Section 3.2). The estimation of each ES value function draws on a different number of observations (see Table 1): Provisioning services (302; 47.5%), Regulating and maintenance services (225; 35.4%), and Cultural services (109; 17.1%).

Table 2 lists and describes the main variables used in the MA. Table 3 provides summary statistics for each of these variables for every service, with the exception of the dummy variables.


**Table 1.** Number of valuation studies, by service and biome, from the ESVD included.

Note: <sup>1</sup> ESVProv = Provisioning Ecosystem Service Values; ESVReg&Main = Regulating and maintenance Ecosystem Service Values; ESVCult = Cultural Ecosystem Service Values; <sup>2</sup> *CSys* = Coastal systems; *CWet* = Coastal wetlands; *CoRf* = Coral reefs; *CuAr* = Cultivated areas; *Dser* = Desert; *FrWa* = Fresh water; *Gras* = Grasslands; *InWt* = Inland wetlands; *Mari* = Marine; *TeFo* = Temp./Bor. forests; *TrFo* = Tropical forests; *Wood* = Woodland.

> The common variables in all models (Provisioning; Regulating and maintenance; Cultural) are Population density (*PDen*) and Gross national income per capita (*GNI* see Table 3). These variables show the largest mean, minimum, and maximum dispersion, representing the large differences in population and wealth in countries around the world. Additional variables were created to describe potentially influential study site characteristics. In the case of Provisioning services, these were: the agricultural areas (*APer*) and the terrestrial protected areas (*TProt*). The former represents the food, fuels, and fibers provisioned, and the latter represents regulation of flows and purification provided. In the case of Regulating and maintenance services, these were: the forest areas (*FPer*) and the terrestrial (*TProt*) and marine (*MProt*) protected areas. These variables express the quality/quantity of natural resources that directly influence their prevention, moderation, and support. In the case of Cultural services, these were the marine protected areas (*MProt*), which represent quality, namely related to the sea.

> > **Table 2.** Meta-analysis variables description.


#### **Table 3.** Summary statistics for meta-regression variables in ecosystem services.



**Table 3.** *Cont.*

Note: <sup>1</sup> See Table 2 for variable descriptions.

#### *3.2. Meta-Regression Model Specification*

We adopt a semi-log functional form specification for the ES value functions, which implies that the marginal effect of a change in ESV depends on income and population density [15].

The Provisioning ES value function is determined by the type of biome (*DBiome*), location of the continent (*DContinent*), terrestrial protected area (*TProt*) [46], percentage of agricultural land (*APer*) [47], population density (*PDen*) [5] and income (*GNI*) [15], and given by:

$$\begin{array}{ccccc} \ln(ESV\_{\text{Prov}}) = a\_0 + a\_1 \ \ast \ D\_{\text{Riouz}} + a\_2 \ \ast \ D\_{\text{Contint}} + a\_3 \ \ast \ Tport + a\_4 \ \ast \ APer + a\_5 \ \ast \ln(PDen) + \end{array} \tag{5}$$

where *α*<sup>0</sup> is a constant, *α*<sup>1</sup> and *α*<sup>2</sup> are dummy regression estimates, and *α*<sup>3</sup> to *α*<sup>6</sup> are variable regression estimates.

The Regulating and maintenance ES value function is determined by the type of biome (*Dbiome*), location of the continent (*DContinent*), level of protection in study area (*FProt*) [15], the terrestrial (*TProt*) and marine (*MProt*) protected area [46], percentage of forest land (*FPer*) [47], population density (*PDen*) [5] and income (*GNI*) [15], and given by:

$$\begin{array}{c} \ln\left(ESV\_{\text{RgckMain}}\right) = \beta\_0 + \beta\_1 \; \* \; D\_{\text{Bione}} + \beta\_2 \; \* \; D\_{\text{Cuntimet}} + \beta\_3 \; \* \; FProt + \beta\_4 \; \* \; FPer + \beta\_5 \; \* \; MProt + \beta\_6\\ \quad \*TProt + \beta\_7 \; \* \; \ln(PDen) + \beta\_8 \; \* \; \ln(GNI) \end{array} \tag{6}$$

where *β*<sup>0</sup> is a constant, *β*<sup>1</sup> and *β*<sup>2</sup> are dummy regression estimates, and *β*<sup>3</sup> to *β*<sup>8</sup> are variable regression estimates.

Finally, the Cultural ES value function is determined by the type of biome (*Dbiome*), location of the continent (*DContinent*), level of protection in study area (*PProt*) [15], marine protected area (*MProt*) [46], population density (*PDen*) [5] and income (*GNI*) [15], and given by:

$$\begin{array}{ccccc} \ln(ESV\_{\text{Cul}}) = \gamma\_0 + \gamma\_1 \, \ast \, D\_{\text{Rione}} + \gamma\_2 \, \ast \, D\_{\text{Contimet}} + \gamma\_3 \, \ast \, \text{PProt} + \gamma\_4 \, \ast \, M \text{Prot} + \gamma\_5 \, \ast \, \ln(PDm) + \gamma\_6 \\ \quad \ast \, \ln(GNI) \end{array} \tag{7}$$

where *γ*<sup>0</sup> is a constant, *γ*<sup>1</sup> and *γ*<sup>2</sup> are dummy regression estimates, and *γ*<sup>3</sup> to *γ*<sup>6</sup> are variable regression estimates.

#### *3.3. Meta-Regression Model Results*

Table 4 reports regression results for two model specifications; the "Full" model in which all variables are included and the "Restricted" model in which non-significant explanatory variables were excluded in a stepwise procedure (applying a cut-off significance level of 20% for the *t*-test). The following base values for the dummies are considered: Grasslands (*Gras*) for biomes; Not protected (*NProt*) for protection status; and Europe (*Euro*) for continents.

The main explanatory variables presented in all "Restricted" models were Population density (*ln\_PDen*) and Gross national income (*ln\_GNI*), with positive coefficient values and high significance (*t*-test < 0.09), which implies that an increase in population or income results in an increase ESV. As we adopt the logarithmic form for these variables, the marginal increase in ESV is decreasing in population or income.

We adopted additional explanatory variables for environmental quality, being *MProt* and *TProt* the percentage of, respectively, terrestrial and marine protected areas. Specifically, for the Provisioning model the *APer* (percentage of agricultural land) and for the Regulating & maintenance model the *FPer* (percentage of forest land) were used.

The Provisioning ES model provides a reasonable fit to the data, although it is the model with the smallest R<sup>2</sup> (0.19) and with the statistics of 0.01 in ANOVA for the restricted model. The signs of the explanatory variables are, as expected, positive for *Dbiome*, *LaAm, ln\_PDen*, and *ln\_GNI*, and negative for *TProt*. This confirms that the other land covers analyzed tend to have a higher value than Grasslands (used as a base for the dummy biomes) and that areas located in Latin America generate larger provisioning ecosystem service values, while the ecosystem service value decreases with an increase in the percentage of protected terrestrial area. The variable *APer* is an exception (*Coef* = −0.04 and *t*-test < 0.01), presenting a negative coefficient, for which a positive sign was expected—which could be explained by the fact that countries with larger agricultural areas present a greater supply of provisioning services, though lower productivity levels. Significant explanatory variables present *t*-test < 0.19, the remaining variables were dropped. Evaluating the dummy variables for biomes, the one that presented the highest coefficient for the ESVProv was *CuAr* (*Coef* = 3.69 and *t*-test < 0.01), indicating that Cultivated areas is the key variable explaining provisioning service values.

The Regulating and maintenance ES model provides a good fit to the data, being the model with the highest R<sup>2</sup> (0.46) and with the statistics of 0.01 in ANOVA for the restricted model. The sign of the explanatory variables is as expected positive for *Dbiome*, *ln\_PDen*, and *ln\_GNI*, and negative for *AFric*. This confirms that, as mentioned before, the other land covers analyzed tend to have a higher value than Grasslands and that areas located in Africa tend to have a lower value for this type of service (due to the lower aggregate income). The variables related to nature protection: *FProt* (*Coef* = −1.73 and *t*-test < 0.01), *FPer* (*Coef* = −0.02 and *t*-test < 0.05), *MProt* (*Coef* = −0.02 and *t*-test < 0.19) and *TProt* (*Coef* = −0.05 and *t*-test < 0.05), present negative coefficients, for which a positive sign was expected, revealing the theory that protected areas, which generally have low population density or are even inaccessible to the population, represent a low monetary value (i.e., people do not fully perceive the value of this service being generated). Significant explanatory variables present *t*-test < 0.19, the remaining variables were dropped. In the ESVReg&Main the largest coefficient for biome was observed in *InWt* (*Coef* = 4.77 and *t*-test < 0.01), although many others such as *CoRf*, *CWet*, and *CSys*, (*Coef* = 4.68; 4.19; 3.98 and *t*-test < 0.01, respectively) also presented high values, these biomes hold a series of important services, such as climate moderation, erosion prevention, maintenance,

and support for different species.

The Cultural ES model also presents a good fit to the data, with an R<sup>2</sup> (0.38) and with the statistics of 0.01 in ANOVA for the restricted model. The sign of the explanatory variables is as expected positive for *PProt*, *ln\_PDen*, and *ln\_GNI*, *LaAm* and negative for *MProt*, *Asia*, *Ocea*. This explains that partially protected areas make it possible for people to access and benefit from the services generated. Moreover, Latin America is the area that presents the largest Cultural ES (primary studies mainly from the Caribbean coast). The *Dbiome* variables *Mari* (*Coef* = −2.47 and *t*-test < 0.12) and *TeFo* (*Coef* = −3.09 and *t*-test < 0.01), present negative coefficients, for which a positive sign was expected, due to the small number of studies related to cultural services involving these land covers in the ESVD. In the ESVCult the largest coefficient was *CoRf* (*Coef* = 2.48 and *t*-test < 0.01), explaining the

high value of services associated with the Coral reefs biome, which provides services such as ecotourism, recreation, and aesthetics, receiving thousands of tourists annually.

**Table 4.** Meta-regression results for Provisioning (ESVProv), Regulating and maintenance (ESVReg&Main) and Cultural (ESVCult) ecosystem service values.


Notes: Dependent variable is ln\_ESV*i*. <sup>1</sup> See Table <sup>2</sup> for variable descriptions; <sup>2</sup> Variable used as the basis for analysis of the dummies; <sup>3</sup> F-test of joint restriction that coefficients of excluded variables are equal to zero.

The model with the least good fit was the Provisioning ES model (R<sup>2</sup> = 0.19), followed by the Cultural ES model with a reasonable fit (R<sup>2</sup> = 0.38) and the Regulating and maintenance ES model" with a reasonably good fit (R<sup>2</sup> = 0.46) for the restricted models. Although these values are low as compared to other ESV meta-analysis studies (see Table 5), a great variability is observed in these studies, with R<sup>2</sup> between 0.25 and 0.87. The explanation for these values is related to the large number of observed studies that presented different characteristics like the location, valuation method, and different years in which the study was performed. For example, [24,30,31] presented large samples, with 682, 416, and 311 observations, respectively. In addition, these studies were applied in wide areas, covering several countries.


**Table 5.** Studies applying the meta-analysis for ESV.

Note: <sup>1</sup> Values presented for the final/best model presented.

As previously exposed, the cut-off for the significance level adopted in the t-student test for the model variables was 20%, which eventually diminished the reliability of the models (i.e., it is common to use "cut-off points" of 0.5%, 1%, 5% or even 10%). Nevertheless, authors such as [7,24,29,41] used *t*-values similar to those adopted in our research. It will be demonstrated, in the next section, that the transfer errors obtained using these value functions are smaller than those obtained using other benefit transfer techniques.

#### *3.4. Value Function Transfer Errors and Estimates*

The validity of environmental benefit transfer has been the subject of a number of studies [7,48,49]. In all of them, the validity has been tested by stating a null hypothesis of no difference between the original study result and the benefit transfer estimate [50]. As in those studies, in this study we seek to verify the differences between the estimated values from MA-BT with the values from the ESVD database, using the Transfer Error technique (see Section 2.3).

#### 3.4.1. Transfer Errors

To assess the accuracy of the estimated ES value meta-models, in order to justify their adoption in future research covering different locations with varied characteristics, we determined the transfer errors associated with Value transfer, Global meta function transfer, and Local meta function transfer (see Section 2.3). This is done for the Provisioning, Regulating and maintenance, and Cultural ES value functions (see Tables 6–8, respectively).

The ecosystem service values and transfer errors per biome related to the estimates for the Provisioning ES are presented in Table 6. Overall, it can be concluded that the transfer error is reduced when moving from Value transfer to Global meta function transfer and, in turn, that the transfer error is further reduced when moving to Local meta function transfer. Notable exception holds for *Wood*, which demonstrates the lowest transfer error when using Value transfer. This is explained by the fact that this variable was dropped from the restricted model (not significant according to the *t*-test). Also, in some cases, the transfer error increases slightly when moving from Global meta function transfer to Local meta function transfer (such as for *FrWa*, *Mari*, and *TeFo*), which is explained by the large variation of values in the ESVD database that contained studies from different countries, continents, and years, and in the case of those biomes, ranging from 1.5 to 3000.0 €/ha/year.



Note: <sup>1</sup> *CSys* = Coastal systems; *CWet* = Coastal wetlands; *CoRf* = Coral reefs; *CuAr* = Cultivated areas; *Dser* = Desert; *FrWa* = Fresh water; *Gras* = Grasslands; *InWt* = Inland wetlands; *Mari* = Marine; *TeFo* = Temp./Bor. forests; *TrFo* = Tropical forests; *Wood* = Woodland.

**Table 7.** Comparison of values and transfer errors (TE) per biome for Regulating and maintenance ES, based on Value transfer, Global meta function transfer, and Local meta function transfer (in 2015 €/ha/yr).


Note: <sup>1</sup> See Table 6 for variable descriptions.

Table 7 presents the ecosystem service values and transfer errors per biome associated with the estimates for the Regulating and maintenance ES. According to the analysis of the previous table, the TE is reduced when moving from Value transfer to Global meta function transfer and then moving to Local meta function transfer. In this case, the exceptions hold for *CSys* and *Wood*, which demonstrate the lowest transfer error when using Global meta function transfer. This is explained by the variation of the values presented in the ESVD database for these biomes. No transfer error is observed for *FrWa* when using value transfer, as only one observation for this biome is available in the ESVD. Finally, no value estimate and transfer error were calculated for *Dser* because there are no primary value estimates data for this biome in the ESVD.


**Table 8.** Comparison of values and transfer errors (TE) per biome for Cultural ES, based on Value transfer, Global meta function transfer, and Local meta function transfer (in 2015 €/ha/yr).

Note: <sup>1</sup> See Table 6 for variable descriptions.

Finally, Table 8 presents the ecosystem service values and transfer errors per biome associated with the estimates for the Cultural ES. Again, it can be observed that the transfer error is reduced when moving from Value transfer to Global meta function transfer and next, to Local meta function transfer. Although there are exceptions, such as for *FrWa*, *InWt*, and *TrFo*, which presented similar TE across Global and Local meta function transfer. One prominent exception holds for *Gras*, which demonstrates the lowest transfer error when using Value transfer. This is justified because it contained only two observations for this biome in the database. No transfer error is observed for *Wood* when using value transfer, as only one observation for this biome is available in the ESVD. Finally, no value estimates and transfer errors were calculated for *CuAr* and *Dser* because there are no primary value estimates for these biomes in the ESVD.

Hence, it can be concluded that transfer errors are reduced significantly when using Global meta function transfer and, in particular, Local meta function transfer as compared to Value transfer. This is justified because value function transfers allow the analyst greater control over differences across sites, they can yield lower transfer errors than simple mean value transfers [51]. In fact, by comparison, value functions offer a greater reflection of the variability of a sample, because the study is dealing with a database with great variability. For this reason, finding a model that, for the most part, has obtained a superior result than other benefit transfers techniques, is an advance that justifies its application given the heterogeneity of the data.

Value functions should thereby draw upon common drivers of preferences reflected in economic theory, including only those variables applicable to all sites [52]. Economic theory suggests that the benefits from environmental improvements should be determined by [53]: (i) change in provision, (ii) distance to the site, (iii) distance to substitute sites, and (iv) characteristics of the valuing individual (in particular income). That is why

Local meta function transfer presents the lowest TE, for addressing these preferences and reflecting the context of each country.

3.4.2. Local Value Function Transfer Estimates

Ecosystem service value estimates per biome for Provisioning (ESVProv), Regulating and maintenance (ESVReg&Main), and Cultural (ESVCult) ecosystem services, are presented in Table 9. Value estimates are thereby based on the restricted models presented in Table 4, using local value function transfer and mean values for the explanatory variables (from Table 3).

**Table 9.** Estimated ES values per biome for Provisioning (ESVProv), Regulating and maintenance (ESVReg&Main), and Cultural (ESVCult) ecosystem services, using Local meta function transfer and mean national values for the explanatory variables (in 2015 €/ha/yr).


Note: <sup>1</sup> ESVProv = Provisioning Ecosystem Service Values; ESVReg&Main = Regulating and maintenance Ecosystem Service Values; ESVCult = Cultural Ecosystem Service Values; ESVTotal = Total Ecosystem Service Values.

> The values found in Table 9 show great variability, with values ranging from ESVTotal = 3.0 €/ha/year for Desert areas to ESVTotal = 1913.5 €/ha/year for Coral reefs. The biomes that provide largest total economic value are Coral reefs (*CoRf* = 1913.6 €/ha/year), Inland wetlands (*InWt* = 1004.2 €/ha/year) and Coastal wetlands (*CWet* = 757.7 €/ha/year). These biomes, in addition to standing out for providing a great diversity of ecosystem services, are also the smallest biomes in terms of the area around the globe and, consequently, the scarcest and, thus, most valuable. In fact, in studies that analyzed ES globally [4,5,54], these biomes were also those with the highest value.

> Provisioning services represent the lowest values and are related to the supplies of products (such as food, materials, or water) with values close to their direct use values [5]. The largest provisioning ES values are provided by Cultivated areas (*CuAr* = 121.9 €/ha/year) and Coastal System (*CSys* = 44.5 €/ha/year). The lowest values were found for Coral reefs, Desert, Grasslands, and Temp./Bor. forests (*CoRf, Dser, Gras*, and *TeFo*, with a value of 3.0 €/ha/year each).

> Regulating and maintenance services are linked to more indirect benefits, which are related to quality, moderation, and prevention in environmental factors (Rao et al., 2015). The largest Regulating & maintenance ES values are provided by Inland wetlands (*InWt* = 425.8 €/ha/year), followed by Coral reefs (*CoRf* = 389.9 €/ha/year) and Coastal wetlands (*CWet* = 238.1 €/ha/year), demonstrating a high added value for areas in transition, notably coastal areas. The lowest values were found for the Marine and Woodland areas (*Mari* and *Wood*, with a value of 3.6 €/ha/year each).

> Cultural services represent the largest values, because they involve complex issues such as aesthetics, generated inspiration, and spirituality, which can be considered incommensurable values as the perception about the environment varies from person to person [5,31]. The largest cultural ES values are provided by Coral reefs (*CoRf* = 1520.7 €/ha/year), Inland wetlands (*InWt* = 555.3 €/ha/year) and Coastal wetlands (*CSys* = 491.6 €/ha/year). The lowest values were found for the Marine areas (*Mari* = 10.8 €/ha/year) and Temp./Bor. forests (*TeFo* = 5.8 €/ha/year).

> It is necessary to be cautious when valuing ecosystem services since, although the aim of pricing is to use values in monetary units, they serve as a tool to provide better insight into the economic benefits of ecosystem goods and services. We do not try to find the

shortcomings and limitations of monetary valuation, both in relation to ecosystem services and man-made goods and services [5,55].

When ESV's models are created and values for biomes are estimated, this does not mean the biomes in question should be treated as private commodities that can be traded in private markets. Most of those ecosystem services are public goods or the product of common assets that cannot, or should not, be sold. Although the flowers, fruits, wood, and leaves enter the market as private goods, the ecosystems that produce them, as for example forests and woodlands, are common assets. Their values are an estimation of the benefits to society expressed in a way that communicates with a broad audience. This can help to raise awareness of the importance of ecosystem services to society and serve as a powerful and essential communication tool to inform better, more balanced decisions regarding trade-offs with policies that enhance the gross domestic product but damage ecosystem services [4].

#### **4. Conclusions and Recommendations**

Ecosystem service value (ESV) meta-models were designed to provide access to values in monetary units for ecosystem services (ES), taking into account the local context of the country and area under analysis. Through their application, it is possible to estimate values for 3 different types of ecosystem services (Provisioning; Regulating and maintenance; Cultural) and 12 different types of land covers (Coastal systems; Coastal wetlands; Coral reefs; Cultivated areas; Desert; Fresh water; Grasslands; Inland wetlands; Open ocean; Temperate/Boreal forests; Tropical forests; Woodlands) in the world. To this end, we built on the review and meta-analysis of the Ecosystem Service Valuation Database (ESVD).

The highest ES values were those associated with Cultural services, followed by Regulating and maintenance and, finally, Provisioning services. Among the biomes with greater associated ecosystem service values are Coral reefs, Inland wetlands, and Coastal wetlands that, among other characteristics, are transitional, aquatic-terrestrial biomes that are scarce and provide a great diversity of services.

It was observed that local independent variables, such as income, population, agricultural and forest area, and those related to the level of environmental protection, are significant explanatory variables and, thus, comprise the ESV meta-models. The application of the meta-functions provides values with greater accuracy as compared to simple value transfer and, as shown by the transfer error analysis, the application of local variables (local meta function transfer) further increases this precision.

A meta-analysis, thus, reduces value transfer errors by taking into account local specifications to determine ESV's. There are several studies that have used meta-models for the valuation of specific ecosystem services and biomes (e.g., [15,29–31]), however, we have not found such a comprehensive study in the literature that has determined the value of 3 ecosystem services for 12 different biomes in the world. Even considering that there are certain transfer errors with the application of meta-models, as compared to other benefit transfer techniques (such as value transfer and value function transfer) the meta-analysis technique has shown to be the best way to estimate the value of ecosystem services.

Some caveats to this study remain. First, there are improvements that can be made to the results, such as updating the database, adopting other explanatory variables, or even different functional forms. Second, the adoption of the ESVD, which although very broad, has some limitations, such as the necessity for further studies for biomes such as Fresh water, which presented only one study for Regulating and maintenance ES, and Woodland, which presented only one study for Cultural ES. Moreover, biomes such as Cultivated areas, Desert, and Marine presented few valuation studies, which could directly influence their estimated ES values. Third, ecosystem services and values from marine biomes face particular challenges as these are scarcely studied and poorly understood (e.g., [56,57]). Finally, it was not possible to estimate values for urban areas, albeit they are important because they have a constant relationship with human well-being through services provided

by areas such as parks, squares, and green spaces, as there were no studies analyzing this land cover in the ESVD database.

We expect this study to be a step further in studies that involve valuing ecosystem services and provide a basis for future research. Not in the least because ecosystem services and values are increasingly considered in environmental planning and nature conservation. Using reliable ecosystem service value estimates from local value functions for 3 ecosystem service types across 12 biomes will facilitate this process—in particular in data-poor circumstances.

**Author Contributions:** Conceptualization, L.M.F., P.R. and M.I.B.; data curation, L.M.F., P.R. and W.R.; formal analysis, L.M.F., P.R. and G.O.; investigation, L.M.F., P.R. and G.O.; methodology, L.M.F. and P.R.; supervision, P.R.; validation, W.R. and P.R.; writing—original draft, L.M.F., P.R. and M.I.B.; writing—review and editing, L.M.F., P.R., M.I.B. and W.R. All authors have read and agreed to the published version of the manuscript.

**Funding:** This research was funded by Coordenação de Aperfeiçoamento de Pessoal de Nível Superior—Brazil (CAPES), who financed part of this work through—Finance Code 001. Thanks are also due, for the financial support, to CESAM (UIDB/50017/2020 and UIDP/50017/2020), to FCT/MCTES through national funds, and to the co-funding by European funds when applicable.

**Institutional Review Board Statement:** Not applicable.

**Informed Consent Statement:** Not applicable.

**Data Availability Statement:** ESVD. Ecosystem Service Valuation Database. Foundation for Sustainable Development. Available online: https://www.es-partnership.org/esvd (accessed on 18 June 2021). World Bank. The World Bank—Databank. World Bank Group. Available online: https://databank.worldbank.org/ (accessed on 18 June 2021). FAOSTAT. Food and Agriculture Data. Food and Agriculture Organization of the United Nations. Available online: http://www.fao.org/ faostat/en/#data (accessed on 18 June 2021).

**Acknowledgments:** We would like to thank the Centre for Environmental and Marine Studies (CESAM) and the Department of Environment and Planning (DAO) at the University of Aveiro for facilitating this research. The authors also thank the anonymous referees for their helpful comments on earlier versions of this paper.

**Conflicts of Interest:** The authors declare no conflict of interest.

#### **References**


MDPI St. Alban-Anlage 66 4052 Basel Switzerland Tel. +41 61 683 77 34 Fax +41 61 302 89 18 www.mdpi.com

*Environments* Editorial Office E-mail: environments@mdpi.com www.mdpi.com/journal/environments

MDPI St. Alban-Anlage 66 4052 Basel Switzerland

Tel: +41 61 683 77 34 Fax: +41 61 302 89 18

www.mdpi.com

ISBN 978-3-0365-2236-4