Spatial distribution and potential ecological risk assessment of heavy metals in agricultural soils of Northeastern Iran

ABSTRACT Illustrating the spatial distribution and potential sources of soil properties, and heavy metals, viz., Fe, Zn, Mn, and Cu, are the vital prerequisites for decreasing their pollution. The 68 composite agricultural soil samples in triplicates were collected by employing grid method to evaluate the concentration of pH, sand, silt, clay, organic carbon, P, K, and heavy metals. Multivariate techniques (Pearson’s correlation, heatmap, principal component analysis, and nonmetric multidimensional scaling), geostatistcal techniques, and contamination indices were employed. The contents of Fe, Zn, Mn, and Cu were lower than the limits for Iran EPA guidelines and Earth’s crust. The results of contamination factor and potential ecological risk index (RI) showed that agricultural soils have less contamination and low ecological risks. The enrichment factor (EF), geoaccumulation index (Igeo) and modified ecological risk index (MRI) indicated that 99%, 86.7% and 52.4% agricultural soil samples showed very high enrichment and ecological risks of heavy metals. Both anthropogenic activities and natural factors were responsible for heavy metal contents. The results of geostatistcal analysis revealed that Zn is accumulated more in Central regions, whereas Cu and Mn accumulated more in South and Northeastern regions of the studied area for EF, Igeo, and modified potential ecological RI.


Introduction
Soil is important nonrenewable resources that act as origin and pool of various contaminants. Soil contamination by heavy metals is an important issue, and various activities such as agricultural, urbanization, and industrialization are responsible for enhancing the heavy metal content in the soil (Hu et al., 2013;Keshavarzi & Kumar, 2018;Kumar et al., 2019;Zou, Dai, Gong, & Ma, 2015). The content of heavy metals in the agricultural soils is a matter of great apprehension due to their accumulative and nondegradable characteristics (Facchinelli, Sacchi, & Mallen, 2001). The pollution of heavy metals in agricultural soils was studied all over the world, including Iran, which is the main issue regarding the possibility of metal absorption by food Tian, Huang, Xing, & Hu, 2017;Zhang et al. 2018). Masoud, El-Horiny, Atwia, Gemail, and Koike (2018) while working on soils of Dakhla Oasis, Egypt applied multivariate and geostatistical techniques and reported that urban and agricultural activities were responsible for degrading the soil quality. Liu et al. (2017) used multivariate techniques, contamination indices, and geostatistical techniques in tobacco growing soils of Shandong, China and inferred that agricultural activities such as application of fertilizers, pesticides, irrigation water, etc., and industrial activities were responsible for heavy metal contamination in the soil. Further, the results of spatial analysis indicated that the distribution of heavy metals was affected by human activities and natural factors. Kelepertzis (2014) while working on heavy metals content in agricultural soils of Peloponnese, Greece found that anthropogenic activities were greatly responsible for heavy metals content. Further, the results of geostatistical analysis showed that high contents of Cu, Mn, and Zn were attributed to citrus soils cultivated for the production of oranges and mandarins.
The determination of heavy metal concentrations will not provide the degrading effect of heavy metals in the environment. In order to determine the pollution and ecological risks of heavy metals, various indices were applied. Scientists applied various factors, i.e., contamination factor (CF), enrichment factor (EF), geoaccumulation index (Igeo), ecological risk index (RI), and modified ecological risk index (MRI) for assessment of pollution and ecological risks (Ahmed et al., 2016;Kumar, Sharma, Minakshi, Bhardwaj, & Thukral, 2018;Tian et al., 2017). The EF and Igeo are based on relative assessment of heavy metals in polluted and unpolluted soil conditions Sakram, Machender, Dhakate, Saxena, & Prasad, 2015). Various researchers have done the spatial distribution of heavy metals by GISbased geostatistical techniques (Ander et al., 2013;Petrik, Thiombane, Albanese, Lima, & De Vivo, 2018;Tóth, Hermann, Szatmári, & Pásztor, 2016;Ungureanu, Iancu, Pintilei, & Chicoș, 2016;Vasu et al., 2017;Zou et al., 2015). The objective was to elucidate the causes for high heavy metal contents and to recognize the main regions for further evaluation.
However, from the last few years the pollution status of agricultural soils have been widely studied but less data were presented about source apportionment and ecological risks of heavy metals in agriculture soils of diverse land types in Mesopotamian countries and mainly in Northeastern Iran. Wheat, barley, corn, etc. are the major crops of Northeastern Iran (Esmaeilzadeh et al., 2018). Khorasan Steel Complex is the biggest incorporated steelmaking plant in Northeastern Iran. In addition to the production of iron and steel, various byproducts, i.e., sludge, oxide layers, slag, metal shell, and dust, are generated and are deposited in surrounding regions of this area which have great influence on the agricultural soil quality (Dogra et al., 2019). Keeping all these things in mind, the present study was done to assess the pH, sand, silt, clay, organic carbon (OC), phosphorus (P), potassium (K), and heavy metals (Fe, Mn, Cu, and Zn) in agricultural soils of Northeastern Iran. The multivariate techniques such as Pearson's correlation, heatmap, principal component analysis (PCA), and nonmetric multidimensional scaling (NMDS) in combination with contamination indices such as CF, EF, Igeo, and ecological RI and modified ecological risk index (MRI) were applied to determine the possible sources, level of contamination, and ecological risks posed by Cu, Zn, and Mn in agricultural soils of Neyshabur plain. Further, the spatial distribution maps of EF, Igeo, and MRI were also done to define the sites that contain high contamination and ecological risks that required further elucidated exploration. The outcomes of this study will present a baseline data about soil quality status in agricultural soils, which may be helpful in protecting the food crop quality and ultimately human health.

Study area
Neyshabur plain, Khorasan-e-Razavi Province, Northeast Iran was chosen for the present study ( Figure 1). It is situated between lat. 35°40′N to 36°4 0′N and long. 58°12′E to 59°31′E with a normal altitude of 1256 m above mean sea level. The climate of the studied area is semiarid with an average annual precipitation of 233.7 mm and an average annual temperature of 14.5°C. The irrigated farming is the major land-use approach in the studied area (Bagherzadeh, Ghadiri, Darban, & Gholizadeh, 2016). The general physiographic trend of the plain extends in the NW-SE direction. Figure 2(a-c) shows the geological map, land types, and land uses of study area. In the studied area, the major geological unit is Qft2, which describes the low level of Piedmont fan and valley terrace deposits. The description of all the geological units is given in Table 1. The major type of land is presented by code 4, which represents Piedmont plain. The irrigated farming and shortgrass rangeland are the dominated land uses in the studied area. As seen in Figure 2(d), Aridisols is the major soil type in the studied area. Based on earlier study (Bagherzadeh et al., 2016), the soils were classified as Aridisols and Entisols orders (Figure 2(d)).

Soil sampling and determination
Using ArcGIS 10.4.1 software, a fishnet sampling design including 300 grids was created ( Figure 1) for an appropriate identification of soil sampling areas to consider the spatial variation of the parameters affecting the agricultural soil quality. The grid interval was 500 × 500 m. During soil sampling, a portable Global Positioning System was applied to exactly find the sampling locations. Twelve locations of 300 grids were urbanized (mostly limited by fence) and were not sampled. A total of 68 sampling sites (0-30 cm) were selected and samples were collected in triplicates as composite samples from different sites. Approximate five subsamples of soil were pooled to make the composite sample. All the samples were stored in clean polythene bags and transported to the laboratory. After that, samples were air dried, grounded, and passed with 2 mm sieve for analysis. pH was measured with the help of digital pH-meter (Thomas, 1996). Soil textural properties of sand (0.05-2 mm), silt (0.002-0.05 mm), and clay (<0.002 mm) were measured by the hydrometer method (Gee & Bauder, 1986). Soil organic carbon (SOC) was determined by following the method of Walkley and Black method (Walkley & Black, 1934). Soil P was measured by using the method of Olsen et al. (1954). K was determined by following the method using extraction with 1 M ammonium acetate (NH 4 OAC) at pH7 (Thomas, 1996). Heavy metals such as Mn, Fe, Zn, and Cu were measured by atomic absorption spectrometer. The soil samples were digested with aqua-regia method (HNO 3 :HCl in a ratio of 1:3) by following the method of Thukral (2016, 2018). The digested samples were filtered and diluted with 20 ml double steam distilled water and used for study. The detection limits of instrument were Mn (1.0 mg/l), Fe (7.3 mg/l), Zn (1.6 mg/l), and Cu (1.2 mg/l). For quality assurance and quality control, the standards and blanks were run after every five samples to check the 95% accuracy of the machine (Arora et al. 2008). The 95-105% recovery rates for samples spiked with standards validated that the results are reasonable (Xiao et al., 2013).

Assessment of contamination level in agricultural soils
The various CFs and indices were used to different heavy metals to find the degree of pollution level and ecological risks posed by heavy metals such as CF, EF, Igeo, potential ecological risk (RI), modified potential ecological risk (MRI), etc. Tian et al., 2017). These pollution indices may illustrate a qualitative threshold or target on ecological risk measurement of individual heavy metals.

CF
It gives us indication about the anthropogenic inputs of heavy metals in the agricultural soils (Ahmed et al., 2016). It was calculated by following the equation given by Hakanson (1980): where HM s and HM b are the concentrations of heavy metals in sampling sites and background environment. The background values of heavy metals were taken from Taylor and McLennan (1995). CF is grouped into four grades: (CF < 1; low contamination), (1 < CF ≤ 3; moderate contamination), (3 < CF ≤ 6; high contamination), and (CF > 6; very high contamination) (Hakanson, 1980).

EF
It presents the enrichment of heavy metals against the content of background heavy metals (Sakram et al., 2015). The heavy metals geochemically differentiating with elevated content in the ecosystem and not competent of presenting antagonism or synergism toward the evaluated heavy metals are used as background heavy metals (Chandrasekaran et al., 2015). Fe was preferred to background heavy metal and EF was computed as where HM s and HM b are the heavy metal concentrations in the samples and background environment, respectively, whereas Fe s and Fe b are the iron concentrations in the samples and background environment, respectively. The EF is categorized into seven types such as (EF < 1; no enrichment), (1 ≤ EF < 3; less enrichment), (3 ≤ EF < 5; moderate enrichment), (5 ≤ EF < 10; moderately enrichment), (10 ≤ EF < 25; high enrichment), (25 ≤ EF < 50; very high enrichment), and (EF > 50; exceptionally high enrichment) (Marrugo-Negrete, Pinedo-Hernández, & Díez, 2017).

Ecological risk assessment (RI and MRI)
The potential ecological RI was computed to evaluate the ecological risk assessment of heavy metals in the agricultural soils. It is defined as multiplication of CF of each heavy metal and toxicological response factor (T r ) of individual heavy metals, viz., Cu (5), Zn, and Mn (1) . It was determined by the following equation: where CF n and T r are the CF and toxicological response factor of individual heavy metals, respectively. In order to determine the ecological risks of anthropogenic and lithogenic additions of heavy metals in agricultural soils, EF replaced for computation of potential ecological RI. The ecological RI computed from EF is called potential MPI . It is defined as the multiplication of EF and T r of individual heavy metals. It was determined by the following equation: where EF n and T r are the EF and toxicological response factor of individual heavy metals, respectively. The grades used for risk assessment are as follows: Er < 40 (low risk), 40-80 (moderate risk), 80-160 (considerable risk), 160-320 (high risk), and >320 (very high risk).

Geostatistical modeling
To illustrate the spatial distributions of contamination indices and factors such as EF, Igeo, and MRI for Mn, Cu, and Zn heavy metals, kriging was accepted as done by various researchers (Liu et al., 2017;Masoud, Koike, Mashaly, & Gergis, 2016). The kriging technique is one of the best linear unbiased techniques that makes satisfactory spatial maps in scanty data area and gives stochastic ambiguity of the maps (Burrough & McDonnell, 2015). The kriging interpolations of EF, Igeo, and MRI were computed by employing ArcGis software to show their spatial distribution maps (Chen et al., 2016).

Statistical analysis
The data were analyzed for various descriptive statistical analysis mean, standard deviation, standard error, coefficient of variation, skewness, and kurtosis in PAST software v. 3.15. The data were also analyzed for PCA and NMDS. PCA was implemented to assess the contamination sources in agricultural soils of Northeastern Iran. It was employed by using varimax rotation with Kaiser Normalization to assess the normalized data after evaluating the compatibility of datasets for PCA factors by employing SPSS v. 16 ([IBM, USA] software Zhang et al. 2018). In NMDS, grade alterations among the sampling sites in multidimensional space are sustained in 2D or 3D space using correlation as similarity measure. Heatmap was prepared by using R programming software v. 3.1.3.

Results and discussion
The descriptive statistical analysis of pH, sand, silt, clay, OC, P, K, Fe, Cu, Mn, and Zn is presented in Table 2. The pH was recorded in the range of 7.5-8.3 in different sampling sites. The slightly alkaline nature of pH is responsible for decreasing the mobility of heavy metals in the soils (Tian et al., 2017). The range of 0.17-0.73% for OC was recorded in the present study and affected the retention of heavy metal in the soils (Troeh & Thompson, 2005). The P and K were found in the range of 2.4-19.4 mg/kg and 73.08-261 mg/kg, respectively. Fe content varies from 20,000 to 550,000 mg/kg in worldwide soils (Bodek, Lyman, & Reehl, 1988) and changes extensively, even within same areas because of soil types.
In the present study, Fe content was found in the range of 2.31-1.24 mg/kg and their low content was attributed to sandy texture of agricultural soils. Normally, Fe content was found low in sandy soils and high in clayey soils (McGovern, 1987). Zn content varies from 10 to 300 mg/kg with a mean value of 50 mg/kg in worldwide soils (Alloway, 2008). In the present study, Zn concentration ranges from 0.28 to 2.84 mg/kg and low concentration of Zn may be attributed to sandy soil and low OC in agricultural soils of this area. The Zn content was also found low as compared to Indian permissible limits of soils, i.e., 300-600 mg/kg (Awashthi, 2000) and 30 mg/kg (European Union, 2009). The mean worldwide upper crustal concentration of Mn is 600 mg/kg and bulk continental crust content is 1400 mg/kg (Taylor and McLennan 1995). The Mn content for the present study was found in the range of 1.64-7.18 mg/kg which is lower than the average upper crustal and bulk continental crust concentrations. The Mn content found in the present study was found lesser than limit of 2000 mg/kg given by European Union (2009). The Mn content was also attributed to low OC and sandy texture of agricultural soils. The range of Cu recorded in the present study was found lower than the limits of Iran EPA Pale-red, polygenic conglomerate, and sandstone E1l Nummulitic limestone E1s Sandstone, conglomerate, marl, and sandy limestone E2c Conglomerate and sandstone E2f Sandstone, calcareous sandstone, and limestone E2m Pale red marl, gypsiferous marl, and limestone E2sht Tuffaceous shale and tuff E3c Conglomerate and sandstone E3m Marl   Mirsal, 2008). The skewness and kurtosis of pH, sand, silt, clay, Fe, and Cu were found less than one and are normally distributed (Beaver, Beaver, & Mendenhall, 2012). The skewness and kurtosis of OC, K, Mn, and Zn were found greater than one and indicate right-handed skewness and leptokurtic (Beaver et al., 2012). The coefficient of variation of heavy metals such as Zn, Mn, Fe, and Cu showed great variations representing anthropogenic activities have great influence on the heavy metals content (Cai et al. 2015). Some soil properties show high degree of standard deviation which may be due to the lack of uniformity of heavy metal distribution in agricultural soils of Northeastern Iran (Arenas-Lago, Pearson's correlation analysis was applied to soil properties to evaluate the relationship among different soil properties (Figure 3). pH showed negative correlation with Mn and Cu, and positive correlation with Zn. The sand content was positively associated with Zn. Clay content showed negative relationship with the Zn content. The P showed positive correlation with K and Mn. Mn and Cu are positively correlated with each other. The high correlations of heavy metals among each other indicated their similar source that may be due to natural as well as human activities. The parent rock material was the major natural sources because these metals are mainly components of Earth's crust (Taylor &    McLennan, 1995). The irrigated polluted water by discharge of sewage and industrial effluents was the major anthropogenic source (Nagajyoti, Lee, & Sreekanth, 2010;Tian et al., 2017). Heatmap was employed to both soil properties and sampling sites (Figure 4). The heatmap indicates that Zn, Fe, Cu, Mn, and OC are included in the same group. K formed the different group. Sand, silt, and clay formed the same group. Heatmap of sampling sites showed that sites 20, 24, 39 65, 2, 27, 46, 3, 30, 32, and 66 are included in the same cluster. The sites 9, 42, 15, 18, 40, 61, 11, 55, 12, 60, 53, 41, 43, 49, 45, 44, 56, 57, 52, and 51 are also included in the same group which may be attributed to similarities in soil properties of these sites.
PCA was also employed to soil properties and heavy metals (Table 3). The first four PC accounted for 73.4% of variation with Eigenvalues found greater than one. The loading of PCA in component matrix explained that PC1 was contributed by soil textural properties (sand, silt, and clay), K, and Zn and explained 26.5% of total variation. The high coefficient of variation of K and Zn indicated that anthropogenic activities have great affect on these parameters. The human activities such as smelting and industrial discharge affected the Zn content  (Singh & Kumar, 2017 (Zhang et al. 2018). The Co-efficient of variation (CV) of Fe and Cu was found low and belong to low spatial unevenness grade. Due to this, both Fe and Cu may be originated from natural sources. Ke et al. (2017) and Liang et al. (2017) also found such type of findings in their studies. The results of NDMS scatter plot showed four points such as point no. 14, point no. 36, point no. 38, and point no. 52 separated from the main group ( Figure 5(a)). Because the stress in NMDS Shepard curve ( Figure 5(b)) is 0.016 less than 0.05, the statistics demonstrated a good fit to NMDS data (Kaur et al., 2018).

Assessment of heavy metal pollution and ecological risks
The anthropogenic effect of heavy metals was assessed by using CF (Figure 6(a)). The CF of all Zn, Cu, and Mn was found less than one, indicating low contamination of these heavy metals in the agricultural soils. EF was also computed to determine the anthropogenic and lithogenic effects of Cu, Mn, and Zn in the agricultural soils ( Figure 6(b)). The results of EF showed that 99% agricultural samples were highly enriched with Cu, Zn, and Mn. The results further showed that each heavy metal showed very high enrichment in the studied area. The Igeo was also calculated for heavy metals (Figure 6 (c)). The results of Igeo showed that 86.7% agricultural samples were highly polluted with Cu, Zn, and Mn. The Igeo of Mn for all the sites showed extreme pollution, whereas Zn and Cu showed (92.6% and 66.1%, respectively) very high pollution in the studied area. The ecological risks of Cu, Mn, and Zn were evaluated by potential ecological RI (Figure 6(d)) and modified potential ecological risk index (MRI) ( Figure 6(e)). From the results of RI, it was found that Er values for heavy metals were found less than 40, indicating low ecological risks of these heavy metals. considerable ecological risks, whereas all the sampling sites for mEr of Cu showed very high ecological risks.
Spatial distribution of soil properties, heavy metals, and EF, Igeo, and MRI Geostatistical analysis presents an equitable determination of different parameters at sites without sampling. Plotting the spatial distribution of agricultural soils is important in estimating the threats of contamination and ecological risks at diverse sites. The spatial interpolation technique such as kriging assessment has been mainly employed in enlightening the spatial distribution of EF (Figure 7(a-c)), Igeo (Figure 7(d-f)), and MRI (Figure 7(g-i)) (Liu et al., 2017). Ten groups were used for each parameter on the basis of their original content. The geostatistical analysis was not employed to CF and RI because values of these indices showed low contamination and ecological risks. The spatial distribution maps of Zn for EF and Igeo showed that it was accumulated more in central regions as compared to Southern and Northeastern regions of the studied area. The spatial maps of Mn and Cu for EF showed that these heavy metals were predominately accumulated in Southern regions of the studied area. Similarly, spatial maps of Mn and Zn for Igeo observed that these are distributed more in

Conclusions
The average concentrations of Fe, Mn, Cu, and Zn were recorded lower than worldwide soils, Indian limits for soil, Iran EPA guidelines, and Earth's crust.
Heatmap and PCA showed that natural sources have great influence on the soil properties in addition to the anthropogenic sources. The results of CF and RI indicated that heavy metals posed low contamination and ecological risks in the agricultural soils. responsible for distribution of heavy metals. The data of this study may help the local government managers to ascertain the pollution reduction approaches and to protect food safety.

Disclosure statement
No potential conflict of interest was reported by the authors.