Enzyme activities in soils under heavy metal pollution: a case study from the surroundings of a non-ferrous metal plant in Bulgaria

Abstract

Soil enzyme activities are potentially valuable indicators of soil health in cases of heavy metal and metalloid (HM) pollution. However, the extent to which HM pollution affects specific enzyme activities remains unclear. In this study we assessed the level of HM pollution and its effects on soil enzyme activities to select the most reliable biochemical indicators of soil health under HM stress. The initial hypothesis was that enzyme responses would be site-specific, depending not only on the level of pollution but also on soil abiotic and biotic properties. The Nemerow pollution index (NPI) and the enzyme activity index (EAI) were calculated. EAIs were based on the activities of dehydrogenase (Dha), β-glucosidase (BGl), acid (AcP) and alkaline (AlP) phosphatases, arylsulphatase (Ars) and urease (Ur). NPI indicated slight (1.0) to serious (67.46) soil pollution. The EAI varied from 4.50 to 11.54. The overall functional dissimilarity between soils was around 26%. Cluster and SIMPER analyses both indicated that the activities of Dha, BGl and AlP grouped together with Cu, Zn and Pb, and accounted for 85% of the between-soil functional dissimilarity. The results indicated that there was enzyme functional redundancy between soils. Soil enzymes showed a relatively high capacity to tolerate long-term HM pollution. Dha, BGl and AlP were highly sensitive to environmental fluctuations, including HM concentrations. The long-term HM stress on soil enzyme activities calls for further studies of the soil properties and the time scale involved in the evolution of metal tolerant enzyme producers in HM-impacted soils.

Keywords:

• Heavy metals
• soil pollution
• Nemerow pollution index
• index of enzyme activities
• soil hydrolases
• soil oxidoreductases

Introduction

Soil pollution with heavy metals and metalloids (HMs) is widespread as a result of industrial and agricultural activities. This highlights HM pollution as a global environmental issue. HMs adversely affect various biological processes in soil [1–3]. Soil organic matter decomposition, mineralisation and nutrient cycling are regulated via the activities of enzymes mainly secreted by diverse soil microorganisms [4]. In this context, soil enzymes can serve as biological indicators for the assessment of soil health because they are sensitive to HM pollution and directly contribute to soil nutrient cycles [3,5]. Metal toxicity may cause: (1) metal complexation with the enzyme substrate making it unavailable for the respective enzyme [6,7]; (2) disruption of the structure and the function of extracellular enzymes by binding with thiol and other groups on the protein molecules, which may replace metals naturally occurring in enzyme prosthetic groups [8], (3) disruption of the deoxyribonucleic acid of the microbial cells and a decrease in enzyme synthesis [9]. Additionally, the bioaccumulation of some metals is an important aspect of their toxicity, which may result in mobilisation through food chains [8,10] and possible effects on higher organisms.

Soil enzymes vary in their responses to HMs, depending on differences in enzyme function, soil properties (pH, organic matter, soil texture), and metal species and concentrations [3,11]. Many studies reported that low concentrations of HMs may stimulate enzyme activities, whereas higher ones are toxic for enzymes [12]. Soil pH, organic matter and soil texture are among the main soil factors influencing the relationships between HMs and enzyme activities [11,13]. Since these are the primary factors that affect the binding of metals to soil colloids and their uptake by biological systems, any changes in these soil characteristics will affect the interactions between heavy metals and soil enzymes [14]. In this context, studies need to analyse enzyme activities under local environmental conditions.

The enzymes which indicate soil quality include oxidoreductases and hydrolases [3]. In the context of HM polluted soil monitoring, the most commonly studied enzyme activities are those of dehydrogenase (Dha), arylsulphatase (Ars), urease (Ur), acid phosphatases (AcP), alkaline phosphatases (AlP) and catalase [5]. For example, researchers used dehydrogenase as a bioindicator in Cd, Cr, Cu, Fe, Mn, Ni, Pb and Zn polluted soils [11], showing a great impact of HMs on its activity. Extensive studies have also focussed on other enzymes, such as β-glucosidase (BGl), invertase and protease, but there are limited reports on their responses to HM pollution [5,15,16].

The aims of the present study were: (i) to assess the levels of soil HM pollution in the surroundings of the non-ferrous metal plant ‘KCM 2000 Group’ (Bulgaria), (ii) to examine the impact of soil pollution on key soil enzymes and (iii) to select the most reliable biochemical indicators for evaluation of soil health under HM stress.

Materials and methods

Study site and sampling procedure

The study area is located in the surroundings of the non-ferrous metal plant ‘KCM 2000 Group’ (Plovdiv, Southern Bulgaria), which is the successor of a factory with a similar industrial activity established in 1961 (Figure 1) [17]. Considering the direction of the pollution spread, five sampling sites were selected as follows: KCM_1 (42°03’31.68"N, 24°49’19.2"E) located 0.5 km south of the smelter, KCM_2 (42°03’5.76"N, 24°49’19.6"E) 2 km south of the smelter, KCM_3 (42°02’6"N, 24°49’19.2"E) 3 km south of the smelter, KCM_4 (42°03’31.68"N, 24°49’45.12"E) 1 km south-east of the smelter and KCM_5 (42°04’23.52"N, 24°49’45.12"E) 1 km east of the smelter. According to the sampling depth, these sites were named KCM_1.1, KCM_2.1, KCM_3.1, KCM_4.1 and KCM_5.1 for the surface layer, and KCM_1.2, KCM_2.2, KCM_3.2, KCM_4.2 and KCM_5.2 for the subsurface soil layer. Soil samples were collected randomly and kept at 4 °C during the sampling, transportation and analytical procedures.

Figure 1. Map of the study area around the non-ferrous metal plant ‘KCM 2000 Group’, Southern Bulgaria and locations of the sampling sites.

Soil physicochemical analysis and HM concentrations

Soil texture was determined via the Kachinsky method [18]. Soil pH was measured in a 1:2.5 soil deionized water slurry using a pH metre (HANNA pH-metre). Total organic carbon (TOC), soil nitrate (NO₃-N), ammonium (NH₄-N) nitrogen and inorganic phosphates (P₂O₅) were determined according to Chen et al. [19], Keeney and Nelson [20] and Olsen [21], respectively. Soil moisture (SM) was calculated after oven drying (105 °C). The concentrations of Zn, Pb, Cd, Cu and As were measured using an inductively-coupled plasma mass spectrometre (ICP-MS) according to EN 16170:2016 (ELAN 5000 Perkin-Elmer, Shelton, CT, USA) after soil decomposition via aqua regia.

Soil enzyme activities

Soil samples in triplicates per sampling site (1 g each) were incubated with the respective enzyme substrate at 22 °C in the dark for 1 h, followed by spectrophotometry of soil supernatant/extract and calculation of enzyme activity. The activities of dehydrogenase (EC 1.1; Dha), β-glucosidase (EC 3.2.1.21; BGl), urease (EC 3.5.1.5; Ur), alkaline (EC 3.1.3.1; AlP) and acid (EC 3.1.3.2; AcP) phosphatases and arylsulphatase (EC 3.1.6.1; Ars) were determined following the methods of Friedel et al. [22], Eivazi and Tabatabai [23], Kandeler and Gerber [24], Tabatabai and Bremner [25], Tabatabai and Bremner [26], respectively. One unit of enzyme activity was defined as µg g⁻¹ h⁻¹.

Data analysis

Each data point in the paper represents the mean of five subsamples with standard deviation (±SD). The Nemerow pollution index (NPI) was calculated to evaluate the levels of soil HM pollution, classifying soils into a five-grade scale: NPI < 0.7: safe soil, 0.7 ≥ NPI < 1.0: warning limit of pollution, 1.0 ≥ NPI < 2.0: slightly polluted soil, 2.0 ≥ NPI < 3.0: moderately polluted soil and NPI ≥ 3: heavily polluted soil [27]. The enzyme activity index (EAI) was calculated as the geometric mean of the studied enzyme activities [28]. One-way analysis of variance (ANOVA) followed by Tukey’s pairwise post-hoc test was performed to examine the significance of the differences in the mean values of soil parameters (pH, NO₃-N, NH₄-N, P₂O₅ and HMs) and enzyme activities (Dha, BGl, AcP, AlP, Ur and Ars) among sampling sites. Pearson correlation analysis was conducted to evaluate the effect of sampling distance from the non-ferrous metal plant ‘KCM 2000 Group’ on the level of soil HM pollution.

The significance of the differences in enzyme activities based on the HM pollution levels and soil depths was assessed through permutational multivariate analysis of variance (two-way PERMANOVA) based on Bray-Curtis similarity index after 9999 permutations. SIMPER (similarity percentage) analysis based on Bray-Curtis similarity index was applied to assess which enzyme activities were primarily responsible for the observed dissimilarity among sampling sites.

Two-way cluster analysis (UPMGA method, Bray-Curtis similarity index, Box-Cox transformed data) was conducted to determine both the relationships between sampling sites, and between soil abiotic and biotic parameters. Statistical analyses were performed using the package PAST [29] at a level of significance p < 0.05.

Results

Soil physicochemical properties

The soil texture was classified as clay (KCM_1.1 and KCM_1.2), loam (KCM_3.1 − 5.1 and KCM_3.2 − 5.2) and sandy loam (KCM_2.1 and KCM_2.2). The soil pH was neutral, with the lowest value of 6.6 at KCM_5.2 and the highest one of 7.2 at KCM_3.1. The total organic carbon ranged from 7.04 g kg⁻¹ (KCM_5.1) to 14.07 g kg⁻¹ (KCM_2.1). The inorganic nitrogen (sum of NH₄-N and NO₃-N) varied from 6.27 mg kg⁻¹ (KCM_3.1) to 50 mg kg⁻¹ (KCM_1.1). The phosphates ranged from 4.18 mg kg⁻¹ (KCM_3.2) to 29.8 mg kg⁻¹ (KCM_5.1). The values (log-transformed) of soil properties are shown in Figure 2A.

Figure 2. Log-transformed values of soil physicochemical properties (A) and HM concentrations (B) near the non-ferrous metal plant.

The level of HM pollution varied among the sampling sites (Figure 2B). The HM concentrations varied for Cu: 60 mg kg⁻¹ (KCM_3.2) − 891 mg kg⁻¹ (KCM_1.1); Zn: 250 mg kg⁻¹ (KCM_3.2) − 4843 mg kg⁻¹ (KCM_4.2); Pb: 130 mg kg⁻¹ (KCM_3.2) − 11569 mg kg⁻¹ (KCM_1.1); Cd: 3.6 mg kg⁻¹ (KCM_3.2) − 185 mg kg⁻¹ (KCM_1.1); and As: 5 mg kg⁻¹ (KCM_2.2) − 191 mg kg⁻¹ (KCM_1.1).

The least variable soil metric among the sampling sites was soil pH (Coef. var.: 2.87), and the most variable were the HM concentrations (Coef. var: Cu − 99.92; Zn − 105.24; Pb − 127.62; Cd − 120.07; As − 125.80).

The NPI indicated that the levels of soil pollution were slight (KCM_3.1 and KCM_3.2) through moderate (KCM_5.1 and KCM_5.2) to heavy (KCM_1.1, KCM_1.2, KCM_2.1, KCM_2.2, KCM_4.1 and KCM_4.2). The highest HM pollution was detected in the surface and subsurface layers of KCM_1 and KCM_4 (Figure 3).

Figure 3. Values of Nemerow pollution index (NPI) calculated for surface layers (blue bars) and subsurface layers (yellow bars) of the soils in the vicinity of the non-ferrous metal plant. KCM_1–5, five sampling sites.

Conducted Pearson correlation analysis indicated significant and negative correlation between the sampling sites’ distance to the source of HM pollution and NPI (-0.71, p = 0.0004). Such relationships were calculated also between the distance of sampling sites and the soil concentrations of Zn (-0.64, p = 0.002), Pb (-0.60, p = 0.005), Cd (-0.60, p = 0.005), Cu (-0.63, p = 0.003) and As (-0.61, p = 0.004).

Soil enzyme activities

We hypothesised that soil enzyme activities would change along the HM gradient, and the rates of changes would be influenced by the local soil properties. The EAIs of soils were calculated, and their values are shown on Figure 4.

Figure 4. Soil enzyme activity indices (EAIs) of surface soil layers (blue bars) and subsurface soil layers (yellow bars) in the vicinity of the non-ferrous metal plant.

Note: Same letters indicate significant differences (p < 0.05) according to Tukey’s pairwise post-hoc test.

The lowest and the highest EAI values were calculated for KCM_1.2 (4.50) and KCM_2.1 (11.54), respectively. The enzyme activities varied between sampling locations and soil depths (Table 1).

Table 1. Enzyme activities^a (µg g⁻¹ h⁻¹) of dehydrogenase, β-glucosidase, alkaline and acid phosphatases, arylsulphatase and urease.

Soil	Sampling depth	Dehydrogenase	β-glucosidase	Alkaline phosphatase	Acid phosphatase	Arylsulphatase	Urease
KCM_1.1.	Surface soil	24.21 ± 0.4	20.82 ± 6.7	18.06 ± 2.4	8.27 ± 2.5	2.57 ± 0.2	1.68 ± 0.1
KCM_4.1.		45.27 ± 6.2	23.59 ± 4.3	23.30 ± 3.2	13.04 ± 2.1	2.50 ± 0.6	1.86 ± 0.6
KCM_2.1.		60.91 ± 7.4	61.71 ± 6.6	49.51 ± 7.6	10.02 ± 1.0	6.82 ± 0.9	2.26 ± 0.2
KCM_5.1.		31.96 ± 5.6	14.73 ± 6.8	39.25 ± 0.3	8.19 ± 0.6	2.24 ± 0.8	3.45 ± 1.0
KCM_3.1.		40.32 ± 6.1	19.68 ± 6.5	21.22 ± 7.6	8.72 ± 1.5	6.98 ± 1.7	2.14 ± 0.2
KCM_1.2.	Subsurface soil	7.10 ± 3.2	28.57 ± 3.4	14.53 ± 1.3	5.48 ± 0.3	1.61 ± 0.5	1.62 ± 0.3
KCM_4.2.		16.84 ± 3.3	16.85 ± 3.9	27.26 ± 1.3	9.49 ± 1.6	2.36 ± 1.1	2.40 ± 0.7
KCM_2.2.		20.32 ± 7.4	33.21 ± 4.3	36.29 ± 7.4	6.60 ± 0.8	5.08 ± 0.8	2.04 ± 0.4
KCM_5.2.		14.05 ± 3.4	8.25 ± 0.8	23.64 ± 7.5	5.17 ± 0.5	3.22 ± 1.6	4.92 ± 0.9
KCM_3.2.		11.95 ± 0.1	8.78 ± 0.6	18.71 ± 0.8	7.90 ± 0.9	5.12 ± 0.6	1.40 ± 0.5

^aData are mean values (±SD).

In general, the enzymes of the surface soil layer were more active than those of the subsurface layer, except Ur, which had depth-independent activity. In particular, the enzyme activities of the surface soil layer were higher by 172% (Dha), 54% (BGl), 39% (AcP), 32% (AlP) and 7% (Ars) than those of the subsurface one. The lowest values of enzyme activities were those of Dha, AlP and Ars at KCM_1.2, as well as BGl and AcP at KCM_5.2. Two-way PERMANOVA indicated that the enzyme activities varied significantly by soil depth (p = 0.004), and non-significantly by the value of NPI (p = 0.057), and the interaction between the two factors was non-significant (p = 0.572). SIMPER analysis per soil layer demonstrated that the functional dissimilarity among soils was around 26%, and the main contributors to this dissimilarity were the activities of BGl, Dha and AlP. The variability of BGl, Dha and AlP activities accounted for more than of 85% in the calculated total functional dissimilarity (Table 2).

Table 2. SIMPER analysis of soil enzyme activities contributing to the functional dissimilarity between soils in the Bray-Curtis dissimilarity matrix.

Surface soil layer (overall dissimilarity: 25.66%)				Subsurface soil layer (overall dissimilarity: 26.46%)
Enzyme	Average dissimilarity	Individual contribution (%)	Cumulative contribution (%)	Enzyme	Average dissimilarity	Individual contribution (%)	Cumulative contribution (%)
BGl	8.92	34.73	34.73	BGl	10.02	37.80	37.80
Dha	7.24	28.24	62.96	AlP	7.40	27.96	65.77
AlP	6.74	26.27	89.23	Dha	5.07	19.16	84.92
AcP	1.29	5.01	98.24	AcP	1.65	6.25	91.17
Ars	1.15	4.49	98.74	Ars	1.51	5.76	96.93
Ur	0.32	1.26	100	Ur	0.81	3.07	100

To understand the similarity between the sampling sites and the relationships between enzyme activities and local soil parameters, we conducted two-way cluster analysis. The cluster analysis grouped the soils according to the level of HM soil pollution, yielding two main clusters: the first cluster included heavily polluted soils (KCM_1, KCM_2 and KCM_4), whereas the second cluster included slightly (KCM_3) to moderately (KCM_5) polluted ones (Figure 5).

Figure 5. Dendrogram of two-way clustering (UPGMA) of soil sites, abiotic and biotic metrics based on the Bray-Curtis similarity matrix.

Additionally, the cluster analysis grouped the soil metrics into four well-defined clusters. The first cluster consisted of three enzyme activities (Dha, BGl and AlP) and HMs (Pb, Zn and Cu); the second cluster grouped together all soil metrics (except inorganic nitrogen), HMs (As and Cd) and AcP; the third cluster contained the different forms of inorganic nitrogen and Ars, and the fourth cluster was formed by Ur alone. This mode of clustering defined close relationships between Dha, BGl and AlP activities and soil concentrations of Pb, Zn and Cu. Additionally, the activity of AcP was related to soil texture and to soil concentrations of As and Cd. Ars clustered closely with the soil inorganic nitrogen, suggesting some dependence on it. Ur activity was not closely associated with any of the analysed soil metrics.

Discussion

The present study focussed on the impact of HM stress on soil enzyme activities in long-term polluted soils in the surroundings of the non-ferrous metal processing plant ‘KCM 2000 Group’. The HM levels varied across the sampling sites, and most of the tested soils were moderately to heavily polluted. Pearson correlation analysis indicated that the main source of soil HM pollution is the long-term activity (since 1961) of the non-ferrous plant ‘KCM 2000 Group’.

It is common understanding that an excess amount of HM in soil exerts toxicity on indigenous biological systems [30]. On the other hand, a more favourable soil nutrient status could mitigate the impact of HM pollution on soil organisms [31]. For example, some of the heavily (NPI: 7.99 − 67.43) and moderately (NPI: 2.88 − 2.93) polluted soils were well abundant in TOC (KCM_2 and KCM_4), inorganic nitrogen (KCM_1 and KCM_2) and phosphates (KCM_2 and KCM_5). This suggested local mitigation of HM toxicity. The values of soil pH varied within a narrow range (pH: 6.6 − 7.2) and were close to the neutral point, suggesting slight to negligible effects of soil acidity on the mobility and toxicity of HMs.

The analysis of different enzyme activities can provide insight into the status of soil processes and functions [3,11,32]. Therefore, we examined six soil enzymes and their relationships to HM pollution and different soil physicochemical properties.

Dehydrogenase catalyses the transfer of hydrogen from organic substrates to inorganic acceptors, playing a central role in the oxidation of soil organic matter [33]. Various studies have described Dha as the most sensitive enzyme to HMs [5,13,34–36]. Herein, we found that the activity of Dha in the long-term HM polluted soils was not only sensitive to HMs, but also dependent on soil depth. These results are in agreement with Wolińska and Stępniewska [37], who identified soil depth as one of the most important factors influencing the Dha activity. Since Dha is strictly connected with the activity of living cells, we suggested that its strong depth dependency might be related to microbial abundance, which declines from surface to deeper soil layers [38–40]. Many authors have reported the ecological role of soil microorganisms as the main enzyme producers [41,42].

β-Glucosidase is a crucial enzyme for the degradation of soil organic matter [43], and can therefore be used as a reliable indicator of changes in the soil organic compounds. Previous studies reported a strong inhibitory impact of HMs, especially Cd, Zn and Pb on the BGl activity [44,45]. In contrast, Narendrula-Kotha and Nkongolo [43] found stimulated BGl activity in HM-polluted soils. Gong et al. [46] also demonstrated the stimulatory effects of Zn and Cu on BGl in a laboratory experiment with microbial producers of ethanol. Our results were in line with Narendrula-Kotha and Nkongolo [43] and Gong et al. [46] demonstrating resistance of BGl to HMs, even stimulatory effects at slight to moderate levels of pollution. This might be related to the enzyme impotence of microbial energy and carbon supply. In cases of pollution, impacted cells need more energy to synthesize enzymes and macromolecules that take part in their recovery process [47].

Urease catalyses the hydrolysis of soil urea into carbon dioxide and ammonia. Ur is highly sensitive to changes induced by external factors and for this reason, it is a general biological indicator of the soil health [48]. Kandziora-Ciupa et al. [49] reported that Cd, Zn and Pb decrease the Ur activity. However, we observed that the Ur activity was not significantly influenced by any of the studied soil factors, including HMs. A previous study [50] reported the resistance of Ur to other soil pollutants like agrochemicals in a mesocosm experiment with increasing concentrations of the fungicide Quadris®.

Soil arylsulphatase is an important enzyme that controls the acquisition of organic sulphur and thus regulates the soil sulphur cycling. The arylsulphatase enzyme activity was related also to the soil concentrations and bioavailability of carbon and nitrogen [51]. In the study, SIMPER analysis demonstrated that the Ars activity, like that of Ur, was the least variable (low contribution to soil dissimilarity) throughout the sampling locations. This suggested enzyme resistance to environmental fluctuations, with the exception of inorganic nitrogen – a soil parameter which was clustered together with Ars by the UPMGA method. Kang and Lee [52] also reported a strong effect of nitrogen on Ars activity.

Phosphatases are responsible for the transformation of soil organic phosphorus into a form suitable for plant use [53]. We observed higher enzyme activity of AlP than AcP, and in the surface- than in the subsurface soil layer. The two enzymes were closely related to soil HM concentrations with enzyme-specific interaction according to the cluster analysis. It grouped the AlP activity with the concentrations of Zn, Cu and Pb, whereas the AcP activity, with As and Cd. Additionally, the cluster analysis showed a closer relationship of AcP with soil texture than with the soil concentrations of HMs. Doelman and Haanstra [54] and Renella et al. [55] also reported a dependency of AcP activity on soil properties, especially soil pH and texture. Our findings also highlighted that AlP had higher contribution to EAI dissimilarity among soil locations (SIMPER) when compared to AcP. We suggested that the low sensitivity of AcP to environmental factors (including HMs) was related to soil pH (range: 6.8 − 7.2), which was not within the pH optimum for enzyme activity. This suggestion is in line with the findings of Renella et al. [55], who reported a linear link between enzyme sensitivity to HMs and soil pH. Overall, our findings demonstrated that in the context of the high NPI variation (1.00 − 67.43): (1) EAI changed within a narrow range (4.50 − 11.54) and (2) the calculated overall enzyme dissimilarity among soil locations was relatively low (26%). These facts indicated high functional redundancy of soil microbial communities, probably owing to the mitigation effects of local soil metrics (pH, organic matter, soil texture, etc.) and/or the potential of soil microbial communities to overcome HM stress via phenotypic or genotypic adaptations throughout the long-term history of soil pollution. However, cluster and SIMPER analyses indicated that three (Dha, BGl and AlP) of the tested enzyme activities clustered more closely with HMs (Zn, Cu and Pb) and explained 85% of the overall functional dissimilarity among soil sites. In this context, we assumed that the biochemical processes catalysed by Dha, BGl and AlP were much more vulnerable to HMs than the others, and this reflects negatively the cell energy metabolism (Dha) and hydrolyses of carbohydrates (BGl) and organic phosphates (AlP). Our results support the findings of previous studies reporting negative effects of HMs on Dha [52], BGl [49] and AlP [56], and provide new insights into the microbial adaptation and functional redundancy under HM stress.

Conclusions

HMs are fundamental factors that influence soil biology. When the concentrations of these elements exceed the background levels, soil functions can be greatly altered. Most of the alterations were negative, especially on the activity of dehydrogenase and alkaline phosphatase, but there were also some stimulatory effects (β-glucosidase). The long-term history of soil pollution and the functional soil redundancy in this study suggest: (1) local mitigation effects on the HM toxicity and/or (2) microbial adaptation to the newly created environments.

Data availability statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Disclosure statement

No potential conflict of interest was reported by the authors.

SM, soil moisture; NO₃-N, nitrate nitrogen; NH₄-N, ammonium nitrogen; P₂O₅, inorganic phosphates; and TOC, total organic carbon. KCM_1–5, five sampling sites; KCM_1.1–5.1, surface soil; KCM_1.2–5.2, subsurface soil.

Additional information

Funding

This work was supported by the National Research Fund of the Bulgarian Ministry of Education and Science under Grant [DN11/4-Dec.2017].