Coupling effect of key factors on ecosystem services in border areas: a study of the Pu’er region, Southwestern China

Location of the Pu’er region

The Pu’er region (99°09′–102°19′E, 22°02′–24°50′N) is in Yunnan Province, China. It covers an area of 45,000 km² and includes nine counties and one district. It is located in the southwestern part of Yunnan Province in Southwest China (Fig. 1) and shares the same border with Myanmar, Laos, and Vietnam. Its sub counties, Lancang, Ximeng, and Menglian have a border with Myanmar with a length of 486 km. Jiangcheng County borders Laos and Vietnam. The China-Laos border and the China-Vietnam border within Jiangcheng County are 116 km and 67 km in length, respectively. The Pu’er region has two national ports and one provincial port: Simao Port, Mengkang Port, and Menglian Port. In 2013, the total cargo volume of these three ports was 374,000 tons. In 2019 this reached 1,049,000 tons, marking a significant increase in the total cargo volume (China Association of Port-of-Entry, 2014; China Association of Port-of-Entry, 2020).

The estimated value of ESs in the Pu’er region administration

The Pu’er region has a subtropical climate influenced by monsoons. The warm, humid air from the Indian Ocean brings a significant amount of precipitation because of the southwest monsoon. Due to ample sunlight and humidity, the forest coverage rate in the Pu’er region is 74.59% (Ma et al., 2020). In 2015, the total value of the ESs in the district was 742.87 billion yuan; the value per unit area was estimated to be 0.163.7 billion yuan/km2 and the per capita value was 280.52 thousand yuan. Ren et al. (2021) reported that the value of the Pu’er region ranked first in Yunnan Province in terms of the total ES value. Provisioning services account for 3.71% of the total value at 27.53 billion yuan/year. Regulating services account for 95.13% at a value of 706.82 billion yuan/year. Cultural services account for 1.16% at a value of 8.64 billion yuan/year.

Transboundary factors and their CEs with local factors

By the end of 2018, the fall armyworm (Spodoptera frugiperda), originally found in Myanmar, had infested the Pu’er region on a large scale. This insect caused significant damage to farmland, resulting in reduced yields of major crops such as maize (Wu, Jiang & Wu, 2019). This reduction was further affected by a strengthened southwest monsoon in southwest Yunnan. The fall armyworm (S. frugiperda), yellow-spined bamboo locust (Ceracris kiangsu Tsai), and red imported fire ant (Solenopsis invicta Buren), which possess strong dispersal abilities, have harmed agriculture and human welfare in many counties in the Pu’er region through their direct impact as transboundary pests and carriers of disease (Chen, Zhang & Chen, 2022; Xiao, Zhou & Quan, 2009; Zhang et al., 2019).

In addition, as cross-border trade increases between Yunnan Province and the neighboring countries (ASEAN Secretariat, 2020), plants, animals, and microorganisms have been brought into the Pu’er region during trade and cultural exchanges. This has posed a greater threat of bio-invasion to the local ecosystem (Zong, Du & Huang, 2020). Road construction development has exacerbated the situation (Xie et al., 2021). This is an example of the CEs’ ability to impact biodiversity (Yin et al., 2020).

Selection of the main types of ESs in the Pu’er region

Based on the characteristics¹ of the Pu’er region (Pu’er Municipal Environmental Protection Bureau, 2020), four main types of ESs (Costanza et al., 1997; Xie et al., 2015) were selected to represent all of the ESs in the study area. The following items were considered:

• Food availability (FA)

Farmland accounts for 12.4% of the total area (Pu’er Municipal Environmental Protection Bureau, 2020); it is the main source for food production. Food availability is the main category within the provisioning services.

• Tourism & recreation (TR)

The Pu’er region has plentiful tourism resources. According to the Institute of Geographic Sciences and Natural Resources, Chinese Academy of Sciences (https://www.resdc.cn/), the Pu’er region contains 397 natural and cultural tourist attractions. Tourism & recreation represent the main source of cultural services.

• Soil conservation (SC)

This region has a complex geographical topography with high mountains and deep valleys. The Ailao Mountain, Wuliang Mountain, Weiyuan River, Xiaohei River, and Amo River experience different degrees of deterioration and soil and water loss. The strategic environmental assessment of the Pu’er region also showed that soil conservation is the main tool used to assess ecological areas. Soil conservation represents the main regulating service in the study area.

• Biodiversity (BD)

Forest covers 74.59% of the Pu’er region. The region contains two national, five provincial, and seven county level nature reserves. It represents one of the most abundant and biodiverse regions in China (Myers et al., 2000). Biodiversity is the main supporting service in the study area.

Factor sorting

After determining the main ESs in the Pu’er region, we chose 15 factors, as shown in Table 1 (all abbreviations for ESs and factors mentioned later are included). Transboundary factors are marked as T and local factors are marked as L (see Appendix A: Table A1 for a detailed description of these factors). Transboundary factors refer to the factors that cross geographical or political boundaries and originate from ecological, economical, societal, and cultural elements. These have direct and indirect impacts on ecosystem services (Liu et al., 2020; Mason et al., 2020). As can be seen from Table 1, cultural exchanges, economic and trade exchanges, land use, population, and other factors impact the various types of ESs, indicating that there is a complex logical network structure between these factors and the ESs.

It should be noted that the transboundary pests and diseases factor was distinguished from the species invasion factor. This was due to the fact that transboundary pests and diseases had a greater impact on farmland ecosystems, while other types of species invasion (such as invasive species like Mexican sunflower (Tithonia diversifolia) which is spread along roads and during cultural exchanges) mainly affect biodiversity.

Establishment of the logical network from correlations of the factors and the ESs

To characterize the logical network between the factors and the ESs with more accuracy, we invited five experts for the criteria and carried out three rounds of consultations with these experts. We identified the factors of each ES using the Delphi method (Vogel et al., 2019), and combined logical networks built by experts with the existing literatures (Dang et al., 2019; Feng et al., 2021) to establish the influencing pathways of each factor on the corresponding ES which reveals the relationship of network.

Valuation of ESs and factors

The data used in this study were selected (see Appendix A, Table A1), and we have selected 2020 as a reference point in our analysis. Some data of the ESs and factors, such as BD, SC, and INV were evaluated based on cumulative data over multiple years or were simulated from public datasets (see Appendix B).

To facilitate training at each node in the BN, we gridded and divided the territory of the Pu’er region into over 1,000 square units of equal size.

Building BN-GIS model for the quantification of the ESs and factors

BN is a graphical probabilistic model with two important components. The first component is the directed acyclic graph (DAG), which links the child nodes to the parent nodes using a set of arrows representing the causal relationships between connected nodes (Landuyt, Broekx & Goethals, 2016). The DAG is the part of BN used for qualitative analysis, which can reveal the logical relation networks of the study. The second component is the conditional probability tables (CPTs) of each node, which quantify the relationship among the variables and specify the state of the parent node in order to determine the degree of belief for a particular state (Forio et al., 2020; Pollino et al., 2007). The CPTs are the part of BN used for quantitative analysis, which is calculated by the baseline training. It supports the sensitivity analysis and the probability distribution adjustment of each node in the network (van der Gaag, 1996).

In this study, we established the BN-GIS model using Netica software, which is a specific program used widely in the construction and simulation of BN (Johnson, Low-Choy & Mengersen, 2012; Smith et al., 2007). We used GIS mainly for assisting in the acquisition of baseline training sample data and for the spatial display of the simulation results (Guo et al., 2020a; Guo et al., 2020b; Jensen & Nielsen, 2007).

Sensitivity analysis for identifying the key factors and their pathway

After baseline training, the conditional probabilities of the nodes were calculated. The sensitivity analysis was used to calculate the reduction in the variance of the real value of the target node, which was based on the conditional probability of each node in relation to other nodes (Pearl, 1988).

The main ESs of the Pu’er region were set as target nodes. The main ESs comprise continuous variable nodes. The degree to which other nodes influenced the target nodes was estimated according to the magnitude of the variance reduction (Wu et al., 2021). This was then used to select the key factors. We identified the factor nodes by withdrawing the ES nodes. The factor nodes with the top 25% variance reduction values were determined to be the key factors using the quartile method (Langford, 2006).

The influencing pathways of each key factor were determined by using the key factor nodes as the starting points and the corresponding ESs nodes as the ending points in the logical relation network.

Application of Netica

In this study, Netica was used to integrate the logical relationship networks of the four ecosystem services into a DAG and to preprocess prior data for the generation of CPTs. These CPTs, along with the DAG, collectively form the Bayesian network for the ESs. Netica was also used for sensitivity analysis of the four ESs to identify the key factors and influencing pathways. Furthermore, it allows for the adjustment of probability distributions of different factors to obtain posterior probabilities for the ESs.

Methodology for the determination of CEs

The formula for CEs is shown in Eq. (1), which was developed based on the physical chemistry of the multiple factors effect (Hestand & Spano, 2018). (1)${\displaystyle \Delta E{S}_{CT}-\sum _1^n\Delta E{S}_{T_f}\left\{\begin{array}{cc}>0\hfill & \hfill \left(Positivecoupling\right)\\ =0\hfill & \hfill \left(Nocoupling\right)\\ <0\hfill & \hfill \left(Negativecoupling\right)\end{array}\right.}$

In Equation (1): when a group of factors f₁, f₂…f_n (n ≥2) acts together on an ES, the CEs on the ES is expressed as $C{T}_{\left(f_1,f_2\dots f_n\right)}$; the individual action of each factor on this ES is T_f1, T_f2…T_fn.

The variation of CT_(f1,f2…fn) is ΔES_CT. The total ES variation, with T_f1, T_f2…T_fn , is ${\sum }_1^n\Delta E{S}_{T_f}.\Delta E{S}_{CT}$ and ${\sum }_1^n\Delta E{S}_{T_f}$ are related as shown in Eq. (1).

A positive result indicates that the CEs can increase the total effect of each factor on the corresponding ES, while a negative result indicates that the CEs can decrease the total effect of each factor on the corresponding ES.

Variations of the expected value of ESs

Using BN can establish a quantitative correlation between multiple factors and the corresponding ES. When we adjust the parametric probability of a factor in a certain state, that of its corresponding ES in each state will also change. This is considered the ES’s response to a single factor. When we simultaneously adjust the parametric probability of multiple factors in a certain state, the parametric change of the ES in each state is seen as a response to the CEs of multiple factors. When this response is applied to each unit in a geographic space, there are EVVs of the ESs.

Therefore, when the high state probability of the factors is adjusted to 100 simultaneously, ΔES_CT is the mean value of the EVV of the corresponding ES for each grid of the study area under the CE of multiple factors. ${\sum }_1^n\Delta E{S}_{T_f}$ is the sum of the mean values of the EVV of the ES for each grid of the study area affected by each factor in the same group above, when one adjusts the highest probability of each factor to 100.

Several permutation and combination scenarios of transboundary and local key factors were used to calculate the CE of multiple factors in different scenarios. Using this approach, we can identify which positive CEs should be promoted and which negative CEs should be closely monitored or addressed through the implementation of appropriate management measures. In this way, managers can develop effective strategies for the various scenarios in which ESs are influenced by CEs.

The expression of spatial variation of the EVVs with that of ESs in the main scenarios

Spatial variations of ESs in some scenarios can effectively illustrate the ecological risk and ESs management needs of each region. The EVVs of ESs were calculated in each grid in the Pu’er region using Netica software with multiple factors combination scenarios. These grids were expressed spatially using GIS to present the spatial heterogeneity of the corresponding ES. By depicting the spatial heterogeneity of the EVVs of ESs, we can pinpoint potential areas where negative CEs can occur. This can assist managers in adjusting management strategies for specific locations.

Baseline training

Once the DAGs of the nodes had been properly established, we integrated four logical relation networks into the initial network structure of BN. After evaluating and gridding each node in the network, the CPT of each node in the network (Fig. 4) was estimated using baseline training.

Figure 4 shows that the nodes of FA, TR, SC, and BD had high probability distributions (4.58%, 31.5%, 34.9%, and 21.7%, respectively), moderate probability distributions (77.4%, 42.0%, 48.4%, and 40.8%, respectively), and low probability distributions (18.0%, 26.5%, 16.6%, and 37.6%, respectively).

Key factors and the main influencing pathways

After the sensitivity analysis of each main ES, the variance reduction of each node to the corresponding ES was calculated (Table 2). After eliminating the ES node of FA, TR, SC, and BD, we used the quartile method (Langford, 2006) to identify the key factors of each ES. There are four key factors for each ES (see abbreviations in bold in Table 2). The influencing pathway of the key factors was then determined (Fig. 5). There are three key influencing pathways of FA and four main influencing pathways of TR, SC, and BD.

As shown in Fig. 5, the key factors for FA and BD include both local and transboundary factors, while SC and TR consist solely of local key factors.

CEs on FA

The key transboundary factor of FA was PD. The key local factors were AP, POP, and LU. Because POP and LU were on the same key pathway (pathway “b” in Fig. 5A) and POP was the parent node of LU, this study selected POP in FA as the main node in pathway “b” to calculate the CEs of the factors.

The transboundary factor PD was combined with the local factors AP and POP. The CE values of the factors on FA were calculated under each combination scenario. As shown in Fig. 6, the number of scenarios with a positive CE is significantly less than those with a negative CE. The highest positive CE is about twice as small as the negative CE. This indicates that FA is more likely to be affected by a large negative CE.

CEs on BD

The key transboundary factor of BD was INV. The key local factors were CON, PRO, and NDVI. Because PRO and NDVI were on the same key pathway (pathway “a” in Fig. 5D) and PRO was the parent node of NDVI, this paper selected PRO as the main node in pathway “a” of BD to be used in the calculation of the CEs of the factors.

The transboundary factor INV was arranged and combined with the local factors CON and PRO. The CE values of the factors on BD were calculated for each combination scenario. As shown in Fig. 7, the number of scenarios with a positive CE is higher than the number of scenarios with a negative CE. The highest positive CE was greater than the number of negative CEs. This indicated that BD was more likely to be affected by a large positive CE.

Spatial expression of the EVVs of the ESs in the main scenarios

In the transboundary and local factor combination scenarios for FA and BD, only the scenarios with the largest and smallest CEs were screened out for spatial expression, as shown in Fig. 8. PD (H) & AP (H) & POP (L) and PD (L) & AP (H) & POP (H) are the maximum and minimum CE scenarios for FA, respectively. INV (H) & CON (H) & PRO (M) and INV (M) & CON (L) & PRO (H) are the maximum and minimum CEs for BD, respectively.

As seen in Fig. 8, the expected value of FA and BD in each grid shows strong spatial heterogeneity in the Pu’er region. The CE on the ES is reflected in the expected values for all grids in the Pu’er region. The heterogeneity patterns are helpful for all levels of ecological management and can be used to adjust conservation strategies for different geographical features within the region.

For example, in the heterogeneity patterns in Fig. 8B, the strongest negative CE on FA occurs mainly in Lancang County and Mojiang County; this situation should be avoided and prevented in these counties. The heterogeneity patterns in Fig. 8C, the region where BD shows the strongest positive CE, appear mainly in Jingdong County, Simao District, and Lancang County. This is information these counties should take account into when considering conservation strategies.