3.1 HM concentrations in soil-rice system
The descriptive statistics on the HM content in rice and soil were shown in Table 4. The mean concentrations of soil Cd, Cu, Pb, Ni, and Zn were 0.21, 28.65, 27.02, 38.5, 98.75 mg kg− 1, respectively, which were greater than the background values and lower than the soil risk screening value of hazardous elements in Chinese agricultural land (GB15618-2018). These mean values were lower than their counterparts reported by previous studies in rice-producing areas of Zhejiang province such as Wenling, Fuyang, and Nanxun (Zhao et al. 2010; Hu et al.,2017; Yu et al. 2019). The maximum values of soil Cd and Pb were six and three times greater than the soil background values, respectively.
Table 4
Descriptive statistics for heavy metals in soils and rice (mg/kg, n = 95)
Metals | Min | Max | Mean | SD | Skewness | Kurtosis | K-Sp | CV (%) | BVa | RSVb |
Cdsoil | 0.07 | 0.81 | 0.21 | 0.09 | 3.67 | 19.81 | 0.000(0.459) | 42.9 | 0.129 | 0.3 |
Cusoil | 11.8 | 87.03 | 28.65 | 13.35 | 1.96 | 5 | 0.001(0.097) | 46.6 | 30.54 | 50 |
Pbsoil | 5.2 | 113.99 | 27.02 | 22.03 | 1.97 | 3.69 | 0.000(0.589) | 81.5 | 30.46 | 90 |
Znsoil | 23.35 | 62.43 | 38.5 | 6.95 | 0.7 | 1.25 | 0.154 | 18.1 | 107.79 | 200 |
Nisoil | 36.02 | 257.35 | 98.75 | 32.05 | 1.72 | 5.71 | 0.0399(0.397) | 32.5 | 36.48 | 70 |
Cdrice | 0.01 | 0.64 | 0.09 | 0.1 | 10.7 | 2.73 | 0.000 (0.845) | 111.1 | n/a | 0.2 |
Curice | 0.94 | 5.34 | 2.98 | 1.08 | -0.69 | 0.14 | 0.004 (0.358) | 36.2 | n/a | 10 |
Nirice | 0.07 | 1.87 | 0.35 | 0.278 | 13.58 | 3.07 | 0.000 (0.998) | 79.4 | n/a | 10 |
Znrice | 14.11 | 34.9 | 22.42 | 3.54 | 1.03 | 0.54 | 0.164 | 15.8 | n/a | 50 |
Notes: SD: standard deviation; Max: maximum; Min: minimum; CV: coefficient of variation; K-Sp: Kolmogorov-Smirnov test for normality, the value in bracket were calculated after naturally-logarithmic transformation.“a”: The background values of heavy metals for agricultural soil in Zhejiang province (Environmental Monitoring Administration of China, 1990).“b”: Risk screening values for soil contamination of agricultural land based on GB15618-2018 Soil environmental quality, National Food Safety Standard Limits of Contaminants in Foods (GB 2762 − 2017); n/a not available. |
The proportions of exchangeable, residual, Fe-Mn oxide and carbonate bound forms of soil Cd were 49.6%, 35.5%, 4.9% and 10.0%, respectively, indicating Cd in paddy soil had high potential bio-availability (Table S3). And the main speciation of Cu, Ni, and Zn was residual, revealing low potential bio-availability.
The mean values of EI of Cd, Zn, Cu, and Ni in the soil-rice system were 0.521, 0.136, 0.031, and 0.245, respectively. The EI values of Cd and Ni were high (exceeding 0.15), in addition, the bio-availability of Cd and Ni were high (Fig.S1), indicating that the two HMs were easily transferred in the soil-rice system (soil to rice), and accumulated in rice plants. The EIs descended in the following order of Cd > Ni > Zn > Cu, which was consistent with some research on paddy soil (Chen et al., 2016; Mao et al., 2019).
The range of rice Cd concentration was 0.01–0.64 mg kg− 1, with a mean value of 0.11 mg kg− 1. The maximum value of rice Cd concentration was about three times greater than the standard value of National Food Safety Standard (GB 2062 − 2017). It indicated that Cd was enriched in rice, which was consistent with other studies (Chen et al., 2018; Zheng et al., 2020; Lü et al., 2021).
The coefficients of variation (CV) value less than 20% indicates weak variation, which means that the variation of target variable is mainly controlled by natural sources, while the CV value more than 50% means strong variation, which is influenced by external factors, such as anthropogenic activity, within this range, it is moderate variation (Fang et al., 2022). The CV of Cd, Cu, Zn, and Ni in rice were 111.1%, 36.2%, 15.8%, and 79.4%, respectively (Table 4). Among them, Cd showed strong variability. Cu, Zn, and Ni showed moderate variability. The CV of rice Cd was greater than that of soil Cd (42.9%), reflecting that soil properties (pH, soil organic matter (SOM), electrical conductivity (EC), etc.), which are external factors for rice heavy metals, also affected the HM concentration of rice, then increased the variability of rice Cd, which was in line with other researches (Wang et al., 2022)
3.2 Risk assessment of HM intake
The Ei values of soil HMs descended as Cd > Cu > Pb > Ni > Zn (Fig. 2). The mean Ei of all HMs less than 40. The RI was less than 100, which indicated moderate potential ecological risk in this area. Some soil samples Cd’s Ei values were > 40, suggesting potential ecological risks existed in some areas.
A risk quantitative calculation model recommended by USEPA was applied to assess the health risk of rice ingestion in the study. The method was widely applied to calculate the health risk from HM ingestion (Yang et al., 2020; Liu et al., 2021). The THQ values for adults and children caused by the HMs exposure route of food intake were in Table 5. The values of Cd for adult and child both were greater than 1 (Fig. 3), which indicating that local rice consuming posed a high non- carcinogenic health risk to the exposed people. The non-carcinogenic risk of children exposed to HMs was higher than that of adults, suggesting that children were more susceptible to HM exposure. This is because higher metabolic rate per unit body in children and their stronger response of the digestive system to poisons (Ma et al., 2018). This was consistent with other research results related to health risks assessment caused by HMs, which all reported that children had higher non-carcinogenic risk compared to adults (Peña-Fernández et al., 2014; Sun et al., 2017; Yadav et al., 2019). In term of carcinogenic risk (CR), due to the lack of reports on Cu and Zn, only the CRs of Cd and Ni were estimated. The results of CR caused by individual HM were listed in Table 5. Cd was the prominent contaminant with a higher CR than Ni, with CR of 2.32E-02 and 3.71E-02 for adults and children, respectively (Fig. 3). The CRs of Cd and Ni ingested from rice consumption were several orders of magnitude higher than acceptable level (1E-06 to 1E-04), 2.83E-02, 4.53E-02, respectively. The result indicated that local residents were under severely CRs. The results were consistent with previous research (Lei et al. 2015; Liu et al., 2021). Due to the high content of Cd in local rice, the CR risk of human body was high, so local residents should try to avoid eating local rice and import rice from safe rice production areas to meet their needs.
Table 5
Intake and health risk of heavy metals by consumed rice
| Non-carcinogenic hazard index (HQ) | | Carcinogenic hazard index (CR) |
Adults | Children | | Adults | Children |
Cd | 3.8 | 6.08 | | 0.0232 | 0.0371 |
Cu | 0.625 | 1 | | n/a | n/a |
Ni | 0.15 | 0.24 | | 5.10E-03 | 8.16E-03 |
Zn | 0.0125 | 0.2 | | n/a | n/a |
Total | 4.59 | 7.52 | | 0.0283 | 0.0453 |
Note: n/a not available. |
3.3 Spatial cluster and outlier of HMs
Significantly positive spatial autocorrelation of soil HMs (P < 0.05) were found, indicating the concentration of soil HMs was similar to that of surrounding locations with particular spatial cluster and outlier patterns. Based on LISA maps of soil HMs (Fig. 4), the High-High value areas of soil Cu, Ni, Zn, and Cd, identified by local Moran's I index, were mostly in the northern and central regions of Shengzhou, where long-term industrial activities existed. Similar spatial distribution patterns of soil Cd, Cu, Ni, and Zn (Fig. 4) indicated these HMs might come from identical pollution sources (Wang et al., 2019). The electrical-component manufactories are widely distributed in northern Shengzhou. The HMs discharged from the manufacturing and discharge processes of these factories were the primary sources of the accumulation of HMs in the local arable soil (Zhao et al., 2015). The Low-value spatial clusters of soil HMs concentrations were mostly in the southwest and southeast. These areas maintained traditional agricultural production with slight soil HM contamination. High values were found next to Low-Low clusters, named as High-Low outliers. According to our field investigation, there were no large-scale industrial activities in these areas, therefore, the High-Low outliers probably were related to the high-way transportation. The Low-High outliers tended to appear near high-high clusters (Lü et al., 2021).
3.4 Spatial structures of HMs
The best fitted models of each element of soil-rice system and corresponding parameters are given in Table S4 and S5, based on the largest coefficient of determination (R2). It was found that the closer R2 to 1, the more fitted model to the target variable (Hu et al., 2019). Soil Cd, Zn and Ni were best fitted using an exponential model. The spherical model was applied to fit soil Cu, and Pb, while the exponential model was applied to fit rice Cd, Cu, and Zn. When the spherical model was used to fit rice Ni it fit the best.
The nugget values of Cd, Cu, and Ni in the soil-rice system were 0.017–0.471, indicating that these HMs were affected by random and inherent factors lightly. The nugget-to-sill (N/S) ratios of soil HMs in the study area ranged from 3–49.2%. The (N/S) ratios of Cd, Ni, and Pb were all less than 25%, indicating a strong spatial autocorrelation existed, which suggested that they were primarily influenced by the variation of internal factors (parent material, soil formation process, and microbial activity, etc.). Generally, HMs in soil with strong spatial autocorrelation were mainly affected by structural factors (Fu et al., 2010). The Cu and Zn in soil had a moderate spatial autocorrelation, with (N/S) ratios of 28.4% and 49.2%, respectively. In the study, soil Cu and Zn may be significantly affected by random factors (fertilizer application and soil tillage) (Zhao et al., 2019). It significantly increased the pollution caused by similar structural factors, thus weakening the influence caused by random factors. It reflected that their variation was controlled by both structural factors and random factors (anthropogenic activities). Soil Cd and Pb had shorter ranges, showing that HMs were more easily affected by external factors (fertilizer application). The rice Cd had the longest range of 5.8 km, indicating that rice Cd was affected by the intrinsic factor. Ni and Zn in rice had ranges of about 2.5 km, which showed a weak spatial correlation. The range value of rice Cu, Ni and Zn were smaller than that of soil Cu, Ni and Zn (Table S4), meaning that the rice HMs were affected by both environmental and intrinsic factors, weakening its spatial dependency (Sadeghi et al., 2006).
3.5 Spatial distribution pattern of HMs
Figure 5 and Fig. 6 showed the spatial distributions of soil HM concentrations in soil-rice system. Their low concentration values were located in the east and west, and high concentration values in the north and center. The spatial characteristics of soil Pb were different from other soil HMs, showing soil in southwestern heavily polluted by Pb and exceeded the limited values regulated by GB15618-2018. The contents of Cu, Ni, and Zn in central areas exceeded the risk screening values. The kitchen utensils and machinery accessories industries are well developed in Sanjie and Xianyan towns. The pollutants discharged by these industries were considered to be the main sources of HMs contamination in local arable land. All HMs concentrations in most of the areas were still at safe levels. In addition, there is a local Chang-Tai expressway that goes straight through the middle of the city center, therefore, emissions from transportation may also be one of the causes of pollution. Overall, the concentrations of soil Cd, Cu, Pb, Ni, and Zn of the whole area all within the risk screening values (MEE, 2018).
A spatial similarity could be found between rice and soil HMs, identifying that the accumulation of HMs in rice was partially affected by the concentrations of soil HMs. The spatial patterns of Cd and Cu in rice were similar, with high values in the middle and east, and the concentration of these two HMs in north and west were low. The high concentration values of Zn were mainly in the north and center, and low concentration value were in the south. The spatial pattern of rice Ni was different from that of Ni in soil, and with its high concentration values were located in northern parts, indicating the HMs transfer from soil to rice was related to anthropogenic factors (soil tillage and fertilizer application) (Cai et al., 2015). Different from the spatial patterns of soil HMs, the spatial distribution patterns of rice HMs were quite a fragment, especially in high concentration value areas, showing island-like distribution, which was consistent with our result of the short range of rice HMs and weak spatial correlation. The concentration of rice Cd (Fig. 6) exceeded the national threshold value (MHRRC, 2017) in east and west, posing threats to crop quality in these fields, which needs an imminent need for attention and measurements to eliminate its ecological risk.
3.6 Spatial characteristics of HM EI and Transfer impact factors
High EI values of Cd were in the northeast, and low EI values of Cd were in the south (Fig. 7). The high EI values of Cu were in the eastern part, and low EI values of Cu were at the south and center. The EI value of soil Ni, Zn, and Cu had similar spatial distribution characteristics, which differed from the soil spatial distribution characteristics, which further identified that soil properties and rice varieties were also important factors affecting the soil HMs transfer and accumulation besides soil HMs concentration.
Some previous researches have reported that soil physicochemical properties, HM speciations, and rice genotypes significantly affect the bio-availability of HMs in soil (Chen et al. 2016; Pan et al. 2015; Papaioannou et al. 2018). In the study, the principal component analysis (PCA) was applied to find the correlation between HM and other soil physical and chemical properties to identify the prominent impact factor. The PCA results of hybrid rice indicated that the 3, 4, 4, 3 principal components of Cd, Cu, Ni, and Zn all cumulatively explain more than 70% of the variance (Table 6), so we determined the number of principal components of Cu, Ni, and Zn were 3, 4, 4, and 3, respectively. The first PCA components of Cd were EI, electrical conductivity, the exchangeable, Fe-Mn oxide and organic bound fraction of HMs, soil pH, OM, which indicated that the transfer of Cd of hybrid rice was mainly related to the Cd fractions, EC, and OM. The variance contribution rates of the three principal components of Cd were 31.5%, 22.8%, and 17.9%. Therefore, the soil properties and HM fractions were the main factor influencing the transfer of Cd in soil, especially, soil pH was the most influential factor. The result was consistent with previous studies (Chen et al., 2016; Liu et al., 2015; Yu et al., 2019). The variance contribution rates of the four principal components of Cu were 26.6%, 24.6%, 15.8%, and 14.5%, respectively. The contribution rates of the four principal components of Ni were 25.8%, 22.9%, 18.6%, and 13.8%, while those of Zn were 29.5%, 25.7%, and 16.8%. The Cu transfer was affected by electrical conductivity, organic matter and pH. The organic matter, exchangeable fraction of HM, and electrical conductivity were the main soil physical and chemical properties having effect on transfer of Ni. Zeng et al. (2011) reported that soil pH and OM had a significant effect on soil non-residual Cu an Zn. Overall, the transfer of HMs in hybrid rice were mainly influenced by non-residual forms of HMs, pH, and organic matter.
Table 6
The PCA results of soil environmental factors for Hybrid rice
Item | Cd | Cu | Ni | Zn |
PC1 | PC2 | PC3 | PC1 | PC2 | PC3 | PC4 | PC1 | PC2 | PC3 | PC4 | PC1 | PC2 | PC3 |
EI | -0.46 | 0.12 | -0.47 | 0.07 | -0.16 | -0.6 | 0.5 | 0.03 | 0.15 | 0.12 | -0.82 | -0.57 | 0.08 | 0.59 |
Exchangeable | -0.9 | 0.09 | 0.32 | / | / | / | / | 0.2 | -0.2 | 0.12 | -0.54 | 0.6 | -0.16 | 0.67 |
Fe-Mn Bound | 0.82 | -0.05 | -0.18 | -0.2 | 0.88 | 0.14 | -0.2 | 0.9 | 0.07 | -0.2 | 0.13 | 0.88 | 0.35 | -0.05 |
Organic Bound | 0.71 | -0.06 | 0.19 | 0.3 | 0.85 | 0.21 | 0.2 | 0.79 | 0.4 | 0.15 | 0.15 | 0.59 | 0.56 | 0.16 |
Residual | 0.14 | -0.12 | -0.82 | -0.2 | -0.93 | -0.1 | -0.2 | -0.32 | -0.2 | -0.5 | -0.14 | -0.92 | -0.17 | -0.27 |
pH | 0.88 | -0.23 | 0.08 | 0.34 | 0.07 | 0.11 | -0.7 | 0.23 | 0.11 | -0.9 | 0.16 | 0.11 | 0.33 | -0.85 |
OM | -0.45 | -0.32 | 0.5 | 0.32 | 0.23 | 0.63 | 0.4 | 0.16 | 0.47 | 0.52 | 0.37 | 0.33 | 0.28 | 0.16 |
EC | 0.64 | -0.08 | 0.32 | 0.15 | 0.21 | 0.76 | 0.1 | 0.32 | 0.16 | 0.06 | 0.78 | 0.69 | 0.05 | -0.07 |
Sand | 0.05 | -0.9 | 0.03 | 0.81 | -0.06 | 0.06 | 0.1 | 0.13 | 0.91 | -0.1 | -0.03 | -0.06 | 0.93 | 0.01 |
Silt | 0.03 | 0.87 | -0.18 | -0.9 | -0.08 | -0.2 | -0.1 | -0.15 | -0.9 | 0.04 | -0.04 | -0.21 | -0.86 | -0.02 |
Clay | -0.27 | 0.89 | 0.05 | -0.8 | -0.15 | 0.04 | 0.3 | -0.47 | -0.6 | 0.32 | 0.24 | -0.15 | -0.85 | 0.22 |
Variance(%) | 31.5 | 22.8 | 17.9 | 26.6 | 24.6 | 15.8 | 15 | 25.8 | 22.9 | 18.6 | 13.8 | 29.5 | 25.7 | 16.8 |
Accumulated variance(%) | 31.5 | 54.3 | 72.2 | 26.6 | 51.2 | 67 | 82 | 25.8 | 48.7 | 67.3 | 81.1 | 29.5 | 55.2 | 72 |
Note: “/” indicate that the contents of heavy metals were not detected. The values in bold mean the main contribution to the PCA components. |
3.7 The Best HMs Transfer Model
Deng et al. (2020) carried out a linear model for predicting As, Pb, Cd, Hg, and Cr concentrations in rice grain and soil using soil pH and heavy metal contents. Liang et al. (2015) also used the logarithmic model for predicting the mobility, bio-availability, as well as transfer of soil Cd. The stepwise multiple regression was applied to estimate the transfer of Pb, Zn and Cd in soil and their influencing factor with high coefficient of regression(R2) by (Elkribi-Boukhris et al., 2022) These two models have high accuracy in predicting soil heavy metal content, availability and other indicator.
The coefficients of determination (R2) differ depending on the fractions of HM and rice genotypes. The model-fitted R2 values of hybrid rice using Cd fractions and soil physicochemical properties were 0.68 and 0.76, respectively. For Japonica rice, the R2 values using Cd fractions and soil physicochemical properties were 0.47 and 0.66, respectively, both of them reached extremely significant level (P < 0.01). The model using soil physicochemical properties was greater than that using HM fractions, which indicated that the driving effect of soil physicochemical properties on the transfer of HMs in soil-rice systems was stronger than that of HM fractions.
The fitted model for Cu in hybrid rice using HM fractions and soil properties reached an extremely significant level (P < 0.01). The fitted model using soil properties for Cu in Japonica rice reached an extremely significant level (P < 0.01). The fitted model for HM fractions in soil reached a significant level (P < 0.05). The fitted model for Ni in hybrid rice using soil properties can reach an extremely significant level (P < 0.01). The fitted model for Ni in hybrid rice did not reach a significant level (R = 0.18, P > 0.05) while using HM fractions in soil. The fitted models for Ni in Japonica rice both reached an extremely significant level (P < 0.01). The transfer model of Zn reached an extremely significant level (P < 0.01), which indicated that the transfer of Zn in the two rice varieties was affected by the HM fractions and soil physicochemical properties. In general, the fitted model of soil properties was better than that of HM fractions.
Among the three model comparisons, the combined model of HM fractions and soil properties achieved the highest fitting degree for predicting EIs of soil-rice system at special rice production areas. However, even if the regression model reached an extremely significant level (P < 0.01), the prediction accuracy of models was still low compared similar models applied for small scale experiments (Liu et al., 2011; Liu et al., 2012). This is due to complexity of soil-water interface, practical management and spatial heterogeneity of soil HMs and nutrients. Therefore, more environmental factors should be included in the models to improve the prediction accuracy.