Effect of landscape fragmentation on soil quality and ecosystem services in land use and landform types

This study was carried out to model the ecosystem services, i.e., habitat quality, carbon storage, nutrient export, and soil export using the InVEST software and soil quality through the integrated quality index across Shoor River basin, southwestern Iran, as well as their relationships to landscape metrics using the Pearson’s correlation and ordered least square regression. The interactions between ecosystem services efficiency, soil quality and landscape fragmentation were also evaluated using the ternary analysis. The results showed that the northern and eastern parts of the study area had more suitability to provide soil quality, habitat quality, carbon storage, and ecosystem service efficiency than the southern and western parts while this trend was the opposite for soil and nutrient exports. The highest soil quality, habitat quality, carbon storage, and ecosystem services efficiency were calculated for the forest's land use and the ridge/tops landform. The highest nutrient and soil exports were also seen in agriculture and rangeland land uses, respectively, as well as for plain and ridge/top landforms, respectively. Statistical analysis confirmed that there were positive and moderate relationships between soil quality to habitat quality and carbon storage and no certain relationships to soil and nutrient export. Moreover, patch density had a negative impact on habitat quality, carbon storage, soil quality, and ecosystem services efficiency and a positive effect on soil and nutrient exports. This trend was the opposite for the largest patch index, patch cohesion index, and effective mesh size. The findings demonstrated that increasing the landscape fragmentation decreased soil quality, habitat quality, carbon storage, and ecosystem services efficiency but increased nutrient and soil exports. Furthermore, the interactions between ecosystem services efficiency, soil quality, and landscape fragmentation were changed under the influence of the dominant land use and landform types across the sub-basins of the study area.


Introduction
Human is highly dependent directly and indirectly on nature to obtain plenty of goods and benefits, referred to as ecosystem services (De Groot et al. 2002;MEA 2005;Zhang et al. 2018). Meanwhile, ecosystem services are critical to human well-being and sustainable socio-economic development around the world (Costanza et al. 2017). They are affected by various ecological and biological characteristics of ecosystems, management conditions, and land use changes (Bryan et al. 2018;Chen et al. 2019). Ma et al. (2021a) evaluated the ecosystem services of soil conservation, habitat quality, carbon storage, and water yield using the integrated valuation of ecosystem services and trade-offs (InVEST) tool. De Leijster et al. (2021) investigated the trajectories of ecosystem services, such as carbon stock, erosion control, habitat provisioning, and pest control with the biotic and abiotic factors in the coffee system.
Soils are a fundamental and complex ecosystem consisted an important part of the natural environment (Delibas et al. 2021;Zhang et al. 2018). Soils establish many ecosystem services, such as purification of water, production of food and raw materials, regulation of climate, storage of carbon all essential for the survival and protection of human life (Costanza et al. 1997;Pereira et al. 2020;Wang et al. 2021). Hence, the essential role of soils is obvious in providing a wide range of goods and services for human well-being, creating sustainable socio-economic development, and determining ecosystem efficiency (Daily et al. 1997;Ellili-Bargaoui et al. 2021). Some basic processes in the soil, such as the nutrient cycles, contribute to providing ecosystem services. These processes reflect the soil support capacity for ecosystem services and can be assessed based on the chemical, physical, and biological properties of soils (Duran-Bautista et al. 2020;Kibblewhite et al. 2008).
Soil quality is an important and measurable characteristic that has a major impact on the function and sustainability of ecological systems and the provision of ecosystem services (Baveye et al. 2016;Soto et al. 2020). It depends mainly on soil fertility and reflects the capacity of soil to maintain and enhance climate quality, increase habitat diversity, produce plant and animal products, and support human health (dos . Assessment of soil quality should be based on a set of physical, chemical, and biological properties that interact with soil function. Various methods have been developed to evaluate soil quality. The most important of which are soil quality indicators (SQI) and fuzzy methods (Askari and Holden 2014;Jahany and Rezapour 2020;Tian et al. 2020). The SQI has been used extensively to evaluate different soil types due to its flexibility and performance (dos Nabiollahi et al. 2017;Raiesi 2017). Soil quality has been also studied concerning land use, topography and morphological features (Davari et al. 2020;Shao et al. 2020). Derakhshan-Babaei et al. (2021) evaluated the relationships of soil quality to erosion, land use, and landform characteristics in the Kan watershed, Iran. Karaca et al. (2021) assessed the changes of soil quality and vegetation intensity for pastures in northeastern Turkey based on the total and minimum data sets and using the multi-criteria decision analysis.
Land use is one of the most important criteria affecting ecosystem services and soil quality. It is closely related to human activities, such as urbanization and agricultural activities (Ma et al. 2021b;Pereira et al. 2020). Land use change alters the soil properties and the structure and function of ecosystems, resulting in a reduction and variation in ecosystem services, soil quality, and biodiversity (da Silva et al. 2018;Lawler et al. 2014). Yang et al. (2021) investigated the effect of land use change on water yield ecosystem services in China's Yellow River basin using the InVEST tool. Kusi et al. (2020) examined the impacts of land use change on carbon storage, soil conservation, and water yield through different scenarios. However, determining the effects of land use on ecosystem services and soil quality has received little attention based on the quantification of landscape characteristics. Landscape metrics can be applied to quantify the land use characteristics in terms of area, shape, composition, distribution, and geometry as quantitative values, which facilitates the analysis, evaluation, monitoring, and spatial planning of land use patterns (McGarigal et al. 2012;Peng et al. 2016).
Landforms are another important land features that affect ecosystem services and soil quality. They have different visual and physical properties that present a lot of information about environmental conditions and resources. Landforms shape the earth's surface in basic spatial units under the influence of various environmental processes over a long period (Derakhshan-Babaei et al. 2021;Du et al. 2019). They can be identified through land geomorphic characteristics, such as altitude, slope, aspect, and curvature using artificial intelligence models, and either supervised or unsupervised classification methods (Du et al. 2019;Li et al. 2020;Zhao et al. 2017).
Due to specific environmental conditions, the Shoor River basin is one of the most important basins located in the northeast of the Khuzestan Province, southwestern Iran. The basin has a semi-arid climate, extensive morphological and topographic changes, and diverse vegetation. These factors have led to diversity and change in characteristics and quality of soils following the provision of various ecosystem services. Human activities and agricultural development have severely affected environmental conditions, natural resources, and ecosystem services of the basin, especially in the southern and western parts. Identifying and evaluating the ecosystem services and soil conditions can be important in accurately understanding the differences and spatial potential of the basin. In addition, it leads to producing new information based on which local managers and decisionmakers can make more appropriate decisions to protect and manage the region. Accordingly, the present study was conducted to analyze the effects of landscape fragmentation on soil quality and ecosystem services according to the following objectives: (i) Modeling of four ecosystem services, i.e., habitat quality, carbon storage, nutrient (phosphorus and nitrogen) export, and soil export using the InVEST software and soil quality through the integrated quality index (IQI). (ii) Determination of the relationship between ecosystem services and soil quality at the sub-basins. (iii) Quantification of the effects of landscape fragmentation on the ecosystem services efficiency and soil quality.

Study area
The study region was the Shoor River basin with an area of 2865 km 2 located in the southwestern of Iran (Fig. 1). The region with a semi-arid climate has an average annual rainfall and temperature of 781 mm and 22 °C, respectively. The basin has severe topographic changes and high roughness having an average altitude and slope of 1136 m and 41%, respectively. Morphologically, the mountains are spread in the northern and eastern parts of the basin, whereas the southern and western parts contain plains. Moreover, forests predominantly cover the northern and eastern parts of the region but agriculture and rangeland are more widespread in the south and west. Forest, rangeland, and agricultural land uses together cover about 97% of the study area, while the rest of the land uses (constructed area and waterbody) include less than 3%.

Data preparation
The list of data used in this research and their descriptions are presented in Table 1. Random sampling was applied to prepare soil data during which 213 samples were taken from a depth of 0-20 cm throughout the study area. The supervised classification method (maximum likelihood) in the ENVI software was used to produce land use and landform maps of the study area (Fig. 1). Landsat 8 multispectral images, downloaded from the United States Geological Survey (USGS) site, were used for mapping the land use. The required data for creating the landform map were elevation, slope, aspect, topographic position index (TPI), topographic roughness index (TRI), valley depth extracted using digital elevation model (DEM) in the SAGA software. The statistics of climatological stations in the region created by the Meteorological Organization of Iran were considered for the preparation of climatic data. The sub-basins map was generated using the DEM in the Arc Hydro environment.

Ecosystem services
Due to simplicity and appropriate efficiency, the InVEST software was used to evaluate habitat quality, carbon storage, nutrient export, and soil export. Data collection and the conditions of the region have been influential to consider these ecosystem services (Ahmadi Mirghaed et al. 2020;Kusi et al. 2020;Li et al. 2021).

Habitat quality
Habitat quality indicates the ability of ecosystems to create suitable conditions for the survival of animal and plant species in ecosystems and is influenced by various factors such as human activities (Ma et al. 2021b;Sun et al. 2019). The InVEST uses land use and land cover (LULC) and threat resource information to estimate habitat quality according to Eq. 1. Threat sources are considered based on four factors: the relative impact of the threat, habitat distance from the threat source, habitat sensitivity, and conservation measures (Hou et al. 2021;Sharp et al. 2020). (1) where Q xj represents the habitat quality in the grid cell x with LULC j. H j and D xj show the suitability and habitat degradation in pixel x with LULC j, respectively. K is a semi-saturated constant. D xj is calculated as follows: ω r is the weight of the threat factor r, which varies from 0 to 1. r y represents the value of the total level of the threat factor in pixel y, which is 0 or 1. i rxy is the effect of threat r in pixel y for habitat in grid cell x. β x is the total level of access in pixel x, and S jx is the sensitivity of LULC j to the threat factor r. The information of threats and their sensitivity for this study is shown in the supplementary data (Tables S1  and S2), respectively, which were prepared based on scientific literature and the current condition of the study area (Ahmadi Mirghaed and Souri 2021; Feng et al. 2021;He et al. 2017).

Carbon storage
Four important sources including soil, surface biomass, underground biomass, and dead organic matter were applied Table 1 Description of soil properties used in this study alongside with fuzzy functions and control points BD bulk density, PD particle density, TP total porosity, OM organic matter, EC electrical conductivity, CEC cation exchange capacity, Ca exchangeable calcium, K exchangeable potassium, Na exchangeable sodium, Mg exchangeable magnesium, Fe total iron, P available phosphorus, N total nitrogen, SAR sodium absorption ratio, AAS atomic absorption spectrophotometer, LI linear increasing, LD linear decreasing, S symmetric  0.9, 9.22) to estimate the carbon storage based on the InVEST model calculated as follows Sharp et al. 2020).
where C t is the amount of carbon (Mg) stored in each pixel. C a , C u , C s , and C d are the carbon pools in the aboveground biomass, underground biomass, soil, and dead organic matter, respectively. Biophysical information required to estimate carbon storage in the study area was obtained based on scientific literature and field surveys (Eggleston et al. 2006;Feng et al. 2021;Ma et al. 2021b).

Nutrient export
The nutrient delivery ratio (NDR) model in the InVEST uses the mass balance approach to describe the travel of nutrient masses (especially nitrogen and phosphorus). The NDR model calculates the ratio of nutrient deliveries transmitted by surface and subsurface flows. Nutrient export for pixel i and its total amount for a basin was calculated as follows (Sharp et al. 2020;Wu et al. 2021): where x expi is the quantity of nutrients (nitrogen or phosphorus) in pixel i. load sur,i and load sub,i are the amounts of nutrient charge per pixel displaced by surface and subsurface flows, respectively. NDR suri and NDR subi are the nutrient delivery ratios in pixel i displaced by surface and subsurface flows, respectively (Sharp et al. 2020).

Soil export
The measurement of erosion in the InVEST software was carried out based on the revised universal soil loss equation where USLE i represents the amount of soil erosion (t ha −1 year −1 ). R shows rainfall erosivity (MJ mm (ha h) −1 ) and K is soil erodibility (t ha h (MJ ha mm) −1 ). LS is the slope length-gradient factor. C and P indicate cover management, and support practice factors, respectively. The SDR i ratio is derived from the conductivity index Sharp et al. 2020). The inputs of the RUSLE model in this study are explained in the supplementary data (Table S3). (3) Ecosystem services efficiency (ESE) was estimated per pixel relative to the whole basin for the ecosystem services of habitat quality, carbon storage, nutrient export, and soil export using the following equation (Asadolahi et al. 2017): where N is the number of ecosystem services types. ES i represents the value of ecosystem services type i in a given grid cell. ES i max and ES i min show the highest and lowest values of ecosystem services type i in the whole study area, respectively.

Soil quality evaluation
The soil quality of the study area was evaluated using the data listed in Table 1 and based on the IQI calculated as follows Mamehpour et al. (2021), Nabiollahi et al. (2017), Raiesi (2017): where n is the number of soil properties. W i and N i are the weight and standardized scores assigned to each soil property, respectively. Fuzzy functions were performed for standardization of indicators and the analytical hierarchy process (AHP) was used for criteria weighting (Derakhshan-Babaei et al. 2021;Hemati et al. 2020).

Standardization of indicators
Each soil property has a special scale. Therefore, a specific scale must be defined for all characteristics to integrate them for the soil quality assessment. In other words, all features must be standardized. Fuzzy membership functions are one of the most important standardization methods that can be used to score soil property values on a scale of zero to one. The values zero and one represent the minimum and maximum proportions, respectively. Standardization of soil properties was performed based on the effect of each feature on soil quality using fuzzy linear functions including increasing linear, decreasing linear, and symmetric in the TerrSet software (supplementary data, Fig. S1). As shown in Table 1, determining the effect of properties on soil quality was defined through control points (a, b, c, and d) adapted from the scientific literature (Derakhshan-Babaei et al. 2021;Eastman 2012).

Criteria weighting
The AHP was used to weigh soil properties. This method has been frequently used in various studies because of its simplicity, flexibility, and appropriate efficiency for determining the weights of criteria for evaluating a problem. The AHP was performed in three general steps in the Super Decision software: (i) Creating a hierarchical diagram based on purpose and physical-chemical properties of soil; (ii) Consisting of a matrix of pairwise comparisons and completing them based on the nine points scale (Supplementary Data, Table S4) and expert's opinions; (iii) Determining the global weights of soil properties and calculating the rate of the inconsistency of judgments. Inconsistency rates less than 0.1 for the AHP results were accepted (Saaty 1980).

Landscape fragmentation
Landscapes are fragmented and changed under the influence of human activities and the development of various land uses. It is possible to study these changes quantitatively using landscape metrics having the ability to quantify landscape in terms of pattern, composition, shape, diversity, continuity, and area at the patch, class, and landscape scales. Various metrics have been developed in this regard. In this study, some of the most important metrics including patch density (PD), patch cohesion index (PCI), largest patch index (LPI), and effective mesh size (EMS), explained in Table 2, were used to investigate the landscape fragmentation of the study area in the Fragstats software (Ahmadi Mirghaed et al. 2020;McGarigal et al. 2012).

Statistical analysis
The Pearson's correlation was used to assess the relationships of ecosystem services and soil quality to landscape metrics. In this regard, their averages in the sub-basins, landforms, and land uses were considered the input of the Pearson correlation in the SPSS software.
The spatial relationships between ecosystem services efficiency and soil quality to landscape metrics at the sub-basin level were measured using the ordinary least squares (OLS) in the Arc GIS 10.5 environment. OLS is a global form of linear regression used to generate a prediction or to model a dependent variable with one or more independent variables. The basic assumption in this method is that the coefficients are constant in the whole study area (ESRI 2016). The relation of the ecosystem services efficiency and soil quality to landscape metrics were also investigated based on the dominant type of land uses and landforms in the sub-basins using the ternary plot in the SigmaPlot software. For this purpose, their values were standardized on a scale of 0-100 and their mean values in the sub-basins were entered as the input of ternary analysis.

Ecosystem services
The produced ecosystem services maps of habitat quality, carbon storage, nutrient export, and soil export are shown A × 10000, PD > 0 n i is the number of patches in the landscape of class i A is the total landscape area (m 2 ) PD clarifies the density of the patches in a specific landscape LPI % LPI = (Maxa ij ∕A) × 100, 0 < LPI ≤ 100 a ij is the area of patch ij (m 2 ) A is the total landscape area (m 2 ) LPI is a simple measure of dominance that quantifies the percentage of the total landscape area contained by the largest patch 0 < PCI < 100 P ij is the perimeter of patch ij (pixel) a ij is the area of patch ij (pixel) Z is the total number of pixel in the landscape PCI estimates the physical connectedness of the corresponding patch type. It increases as the patch type becomes more aggregated or clumped in its distribution A is the total landscape area (m 2 ) EMS is based on the cumulative patch area distribution and is perceived as the size of the patches when the landscape is subdivided into the patches in Fig. 2. The results indicated that the suitability of habitat quality and carbon storage in the northern and eastern regions of the study area were higher than in the southern and western regions, whereas for nutrient export, and soil export, the opposite trends were observed. Carbon storage, nutrient export, and soil export were estimated at 193 Mt year −1 , 328 t year −1 , and 0.59 Mt year −1 , respectively, with averages of 67.4 Mg ha −1 , 1.4 kg ha −1 , and 2.6 t ha −1 , respectively. Figure 2 also illustrates that the efficiency of ecosystem services in the northern and eastern parts of the study area was more favorable than in the southern and western parts. For the average of ecosystem services in the landforms and land uses of the study area as presented in Table 3, the highest and lowest habitat quality is observed in the forests and built-up areas, respectively. The same is true for carbon storage, which averages 119 Mg ha −1 in the forests but there is almost 6 Mg ha −1 in the built-up areas. The maximum (3.2 kg ha −1 ) and minimum (0.81 kg ha −1 ) nutrient export were estimated for the agriculture and forest land uses, respectively. Also, rangeland and forest land uses have the highest and lowest amount of soil export with 41.5 and 8.8 t ha −1 , respectively.
In terms of landform type, the highest habitat quality was seen in the ridge/tops while the lowest value of it was spotted in the plains. The maximum amounts of carbon storage (88.6 Mg ha −1 ) and soil export (4.3 t ha −1 ) were found at the ridge/top landforms, but their lowest values were seen across the plains (28.8 Mg ha −1 and 0.5 t ha −1 , respectively). The highest (2.18 kg ha −1 ) and lowest (1.07 kg ha −1 ) nutrient export was also estimated in the plain and ridge/top landforms, respectively.
The efficiency of ecosystem services based on different landforms and land uses is presented in Table 3. The lowest and highest efficiencies are observed in the plains and ridge/tops landforms, respectively. The same trend was clarified to the built-up area and forest land use, respectively.

Soil quality evaluation
The weight of soil properties calculated according to the AHP method is brought as Table 4. Organic matter (OM) and particle density (PD) with weights of 0.208 and 0.004, respectively, are the most and least important criteria for evaluating the soil quality of the study area. The soil quality map is illustrated in Fig. 2. The lowest quality is observed in most of the soils of the northern section and parts of the south and west of the region, while the eastern and central parts had a better soil quality. The highest and lowest average of soil quality were observed in the forest and built-up areas (Table 3). In terms of landform type, the lowest and highest average of soil quality were considered for the plains and ridge/tops, respectively (Table 3).

Statistical analysis
The results of the Pearson's correlation confirmed that the relationships between habitat quality and carbon storage were positively significant at P value < 0.01 (R = 0.8), while their relationships to nutrient export were negatively significant at P value < 0.01 (R = − 0.53, R = − 0.50). Soil export showed a positive and relatively strong correlation (P-value < 0.01, R = 0.57) with nutrient export, while no correlation to habitat quality and carbon storage was found. The results indicated that soil quality had a positive and significant relationship to habitat quality and carbon storage, but no correlation was found to nutrient export and soil export (Table 5). Figure 3 illustrates the relationships between the ecosystem services and soil quality to landscape metrics at the sub-basin level based on Pearson's correlation. It was found that the patch density (PD) had a direct non-significant correlation to soil export and nutrient export (P value > 0.05, 0.23 < R < 0.34), while its relationship was an inverse significant to habitat quality, carbon storage, and soil quality (P value < 0.01, − 0.53 < R < − 0.4). The relationships between largest patch index (LPI), patch cohesion index (PCI), and effective mesh size (EMS) to habitat quality  Habitat quality (HQ) 0.88** 1 Soil export (SE) 0.12 0.13 1 Nutrient export (NE) − 0.5** − 0.53** 0.57** 1 Soil quality (SQ) 0.50** 0.43* 0.10 − 0.10 1 and carbon storage were positive and relatively strong (P value < 0.01, 0.51 < R < 0.73), while they are inversely and moderately related to nutrient export (P value < 0.01, − 0.52 < R < − 0.49). Soil export had an indirect and nonsignificant correlation to the LPI, PCI, and EMS metrics. Soil quality also had a positive correlation to the PCI and LPI metrics (P value < 0.01, 0.40 < R < 0.53), while it had a direct and insignificant relationship to the EMS metric (P value > 0.05, R = 0.38). The OLS regression results were presented in Fig. 4. It was clarified that the ecosystem services efficiency had a direct significant relationship to EMS, LPI, and PCI metrics with R values of 0.66, 0.63, and 0.53 (P value < 0.01), respectively, while it had a direct significant relationship to the PD metric (R = − 0.49, P value < 0.01). The relationships of soil quality to the LPI and PCI metrics were positive and relatively strong having R values 0.40 (P value < 0.05) and 0.53 (P value < 0.01), respectively, while it had inversely and moderately related to the PD metric (P value < 0.01, R = − 0.56). Soil quality also had a non-significant relationship to the EMS metric. The interactions between the ecosystem services efficiency, soil quality, and landscape metrics based on the dominant landform and land use groups in the sub-basins are shown in Fig. 5, respectively. The results confirmed that the interaction between ecosystem services and soil quality to the PDI and EMS than to the LPI and PCI had more variability based on the dominant land use and landform types at the sub-basins.

Discussion
The results clarified that the provision of ecosystem services and soil quality in the Shoor River basin were different according to the variety of the land use patterns and landforms. In addition, the ecosystem services variation and soil quality changes occurred under the impact of anthropogenic and agricultural activities, especially in the southern and western parts of the region. The northern and eastern parts of the study area had more suitable conditions than the southern and western parts in terms of habitat quality and

Fig. 4
Results of OLS regression between ecosystem services efficiency (ESE) and soil quality (SQ) with landscape metrics (PD patch density, LPI largest patch index, PCI patch cohesion index, EMS effective mesh size) at sub-basin scale. Jarque-Bera (JB) statistic shows the model bias. If JB is significant (p value < 0.01), it indicates that the residuals are not normally distributed carbon storage. The rate of nutrient export and soil export in the southern and western parts was higher than in the northern and eastern parts. It was also found that the efficiency of ecosystem services in the northern and eastern parts was higher than in the south and west. These conditions are created because of the denser vegetation in the northern and eastern parts and the intensive agricultural and human activities in the southern and western parts. Order of the land uses' value/importance was as follows: Forest > Rangeland > Agriculture > Built-up areas for carbon storage and habitat quality; Agriculture > Rangeland > Builtup areas > Forest for nutrient export; and Rangeland > Agriculture > Built-up areas > Forest for soil export. This trend confirms that forests play an important role in increasing habitat quality and carbon storage. Conversely, built-up areas and agricultural activities reduce the ecosystem services of habitat quality and carbon storage and increase nutrient export. The highest amount of soil export occurred in rangelands because the relatively low density of vegetation and relatively steep slopes were among the important factors for increasing soil export in the rangelands. Other studies pointed to the effect of land use patterns on ecosystem services changes (Aghsaei et al. 2020;Ahmadi Mirghaed et al. 2018;Yee et al. 2021).
In terms of the landform type, it was found that the highest habitat quality, carbon storage, and soil export were in the ridge/tops and the lowest in the plains. Instead, the highest amount of nutrient export is estimated in the plains and the lowest for the ridge/tops.
It was also clarified that organic matter, phosphorus, total nitrogen, clay content, and cation exchange capacity (CEC) were the most important indicators of soil quality in the study area. Hemati et al. (2020) confirmed that soil organic carbon and total nitrogen were the most important indicators of soil quality. Ma et al. (2020) also pointed out that organic matter, nitrogen, CEC, and clay percentage were among the most important factors in soil quality. It was found that the soils of the central and eastern parts of the basin had higher quality than other parts, which can be attributed to more diverse and relatively denser vegetation to preserve soil nutrients and increase soil fertility (Hemati et al. 2020). In addition, the intensity of anthropogenic and agricultural activities in these areas was relatively low, which was another factor to increase soil quality.
The highest soil quality was observed in the forest land use but it was the lowest among agricultural and built-up areas. Gong et al. (2015) noted that soil quality varied according to land use pattern and it was higher in natural forests than other land uses. In terms of landform type, the highest soil quality was observed in the ridge/tops, while the lowest quality was found in the plains where human activities including agricultural practices were relatively intense. Derakhshan-Babaei et al. (2021) pointed out that soil quality in rangelands was higher than in agricultural and man-made areas. In addition, they confirmed that landform type and topographic features may change soil quality.
The relation of landscape fragmentation to ecosystem services and soil quality was evaluated using the Pearson's correlation in the SPSS software and based on some landscape metrics including patch density (PD), largest patch index (LPI), patch cohesion index (PCI), and effective mesh size (EMS). The PD metric exhibits the land use density in a certain area, whereas the higher the PD, the more fragmented the landscapes. The LPI, PCI, and EMS metrics refer Results of the interaction between ecosystem services efficiency (ESE), soil quality (SQ) and landscape metrics (PD patch density, LPI largest patch index, PCI patch cohesion index, EMS effective mesh size) based on the dominant landform type (P plain, V valley, GH gentle hillslope, SH steep hillslope) and the dominant land use type (F forest, R rangeland, A agriculture) in each sub-basin to the dominance, continuity, and cumulative distribution of patches in a landscape, respectively, and the higher they are the less landscape fragmentation is expected (Ahmadi Mirghaed et al. 2020;McGarigal et al. 2012). The results illustrated that there was a direct relationship of the PD metric to nutrient export and soil export, while its relationship to habitat quality, carbon storage, and soil quality was inverse. These cases represented the opposite result about the EMS, PCI, and LPI metrics. Habitat quality, carbon storage, and nutrient export had the highest correlations to the EMS metric. Soil quality also presented the highest correlation to the PD metric. Accordingly, it can be concluded that further landscape fragmentation reduces ecosystem services, habitat quality, carbon storage, and soil quality but it increases nutrient export and soil export. Other studies also assessed the relationship between landscape metrics and ecosystem services. Zhang et al. (2018) examined the effect of land use patterns and landscape changes on soil erosion using 12 landscape metrics. They concluded that the largest patch index is the most important metric influencing soil erosion. Ahmadi Mirghaed and Souri (2021) considered the relationship between habitat quality and land use patterns using some landscape metrics. They confirmed that habitat quality increases in the forest and rangeland while it is decreased for agricultural and built-up areas.
The OLS results confirmed that the ecosystem services efficiency and soil quality can be affected by landscape fragmentation. It was found that the EMS metric had the greatest effect on ecosystem service efficiency compared to other indices, while soil quality is more affected by the PD metric. The interaction between the ecosystem services efficiency, soil quality, and landscape fragmentation was investigated using a ternary plot at the sub-basin level. It was revealed that this interaction can be different in the sub-basins with the similar dominant of land use and landform types. The PD and PCI metrics showed the highest and lowest differences for the interaction between ecosystem service efficiency and soil quality in the sub-basins, respectively. These results confirmed that the interaction between ecosystem services and soil quality can be different under the influence of the dominant land use and landform types in the sub-basins across the study area.

Conclusion
This study clarified that the northern and eastern parts of the study area were prioritized for providing habitat quality, carbon storage, and soil quality but the southern and western parts for nutrient and soil export. It was acknowledged that land use and landform patterns affect the change of ecosystem services and soil quality so that the highest habitat quality, carbon storage, soil quality and ecosystem services efficiency were obtained for the forest's land use and the landform of ridge/tops. The highest nutrient export was seen in agricultural land use and plains while the highest soil export was obtained for rangeland and ridge/top landform. Statistical analysis confirmed the existence of a positive and moderate relationship of soil quality to habitat quality and carbon storage while it had no certain relationship to soil and nutrient export. It was also revealed that patch density had a positive effect on soil and nutrient export while it had a negative impact on habitat quality, carbon storage, soil quality, and ecosystem services efficiency. This trend was obtained inverse for the largest patch index, patch cohesion index, and effective mesh size. It was concluded that the higher the landscape fragmentation, the less the habitat quality, carbon storage, soil quality, and ecosystem services efficiency but the more the nutrient and soil export. In addition, it was found that the interaction between ecosystem services efficiency, soil quality, and landscape fragmentation was different among sub-basins of the study area because these interactions may change under the influence of the dominant land use and landform types.