Ecosystem services (ESs) are defined as all benefits human beings obtain from nature that are essential to their lives, providing various goods and services. The ESs are categorized into four services: 1) providing services such as timber production and food, 2) regulating services such as carbon storage and water purification, 3) supporting services such as biodiversity and nutrient cycle, and 4) cultural services such as recreation and aesthetic (Millennium Ecosystem Assessment 2005).

Human activities for enhancing certain ES such as clear-cut, were prevalent in the past. However, it turned out that those activities could unintentionally result in reductions in other ESs these days (Bennett et al. 2009; Egoh et al. 2008). Accordingly, it has gradually increased the need to protect and conserve ecosystems in the decision-making process for resource management. Simultaneously, there is a need to maintain a sustainable balance between benefits and costs in human activities for the short and long term (Carpenter et al. 2009; Liu et al. 2015). Despite this cognition, undesired results happen when silvicultural applications are conducted, which is attributed to a lack of knowledge or poor understanding of the ecosystem’s characteristics providing various and complex functions (Bennett et al. 2009; Chan et al. 2006; Peterson et al. 2003). Therefore, it is critical to deeply understand the relationships among ESs when forest applications are applied to increase desirable ESs. For instance, timber production increased by thinning, yet a decrease in aboveground carbon (AGC) storage (Zhou et al. 2018) and a reduction in biodiversity (Lafond et al. 2017) were found; in other words, trade-offs (TOs) occurred between the paired two ESs. Trade-offs indicate that when an ES increases, other ESs decrease between two or more ESs. On the contrary, synergy represents two or more ESs growing together. Trade-off analysis is a crucial method to study interaction, trends, and driving forces between ESs and helps understand ecosystem dynamics and mechanisms (Bennett et al. 2009; Maron and Cockfield 2008).

Various analysis methods such as descriptive, correlation, and multivariate analyses have been performed for TO analysis; however, those methods failed to quantify the TOs. Bradford and D’Amato (2012) suggested root mean squared error (RMSE) as a simple and effective method for quantifying TOs. Lu et al. (2014) reported that AGC storage and plant diversity were in a TO relationship, but soil total nitrogen and soil organic carbon in a synergy relationship using RMSE in the Loess Plateau in China. Liu et al. (2019) also revealed the TOs between net primary production, water yield, soil conservation, vegetation cover by elevation, and precipitation.

Although payments for and value assessments of ESs (Ahn and Rho 2016; Choi and Oh 2018) and model-based prediction for future ESs (Cho et al. 2023; Kang and Tenhunen 2010; Kim et al. 2021) have been mainly studied in Korea these days, the research on temporal change of TOs between ESs based on measured data is insufficient.

The study compared and analyzed the changes of relative benefits and TOs between AGC storage and species richness (SR) before and after thinning on Pinus koraiensis plantations in Chuncheon. The study’s objectives were 1) to measure the change of the TOs between two ESs and 2) to find the proper forest management method for increasing synergy with decreasing TOs.

Study site

The study was conducted in around 44 and 64-year-old Korean pine plantations at 330–500 m above sea level (N37°52´49.12, E127°52´30.11) on Mt. Gari, Chuncheon (Fig. 1). The annual mean air temperature was 11.1°C (± 0.12°C), ranging from –4.6°C in January to 24.6°C in August. Annual mean precipitation was 1,347.3 mm (± 70.8 mm), with 20.3 mm of the lowest in January and 383.8 mm of the highest in July in Chuncheon from 2000 to 2020 (Korea Meteorological Administration 2023). The site was generally west-facing with 24°–26° of slope, and its total area was about 118 ha. The site was divided into three regions, relying on thinning and stand age: 1) region A with P. koraiensis planted in early 1980 and thinned (20%–25% reduction) in 2002, 2) region B with P. koraiensis planted in early 1960 and thinned (20%–25% reduction) in 2002, and 3) region C with P. koraiensis planted in early 1960 without thinning application in 2002 (Table 1). The dominant woody species was P. koraiensis across three regions. The mean DBH and stand density regions A, B, and C in 2006 were 22.5 (± 0.97) cm, 28.9 (± 1.28) cm, 21.3 (± 2.15) cm, and 572.3 (± 22.4) trees/ha, 412.7 (± 28.5) trees/ha, and 1,050.2 (± 112.6) trees/ha, respectively. The experimental thinning was applied from 2007 to 2008 with three intensities: 1) light thinning (LT) with ~35% stand density reduction, 2) heavy thinning (HT) with ~65% stand density reduction, and 3) control (Con). Generally, the dominant species in 2006 across regions were Zanthoxylum schinifolium, Morus bombycis, and Aralia elata (Miq.) in shrubs layer and Oplismenus undulatifolius, Rubus crataegifolius, and Lindera obtusiloba in herbs layer, respectively.

Table 1 . Site description before thinning in 2006.

Region	Age class in 2006	Thinning history	Plot	Thinning intensity	Total area (ha)	SD before thinning in 2006 (No./ha)	DBH before thinning in 2006 (cm)
A	III	In 2002	A-1	LT	7.0	563	22.3 ± 0.6
			A-2	HT	4.9	600	21.9 ± 0.4
			A-3	Con	2.1	513	25.2 ± 0.5
B	V	In 2002	B-1	LT	18.5	363	30.8 ± 0.8
			B-2	HT	16.4	488	25.2 ± 0.6
			B-3	Con	3.7	425	29.7 ± 1.1
C	V	No thinning	C-1	LT	13.9	1,100	23.9 ± 0.7
			C-2	HT	14.9	726	24.7 ± 0.9
			C-3	Con	3.9	1,250	15.2 ± 0.2

Values are presented as mean ± standard error.

SD: stand density; DBH: diameter at breast height; LT: light thinning; HT: heavy thinning; Con: control.

Figure 1. The location of the study site in Mt. Gari, Chuncheon.

Forest inventory and vegetation survey

A 20 m × 40 m plot was set up on each treatment in 2005. Forest inventory in February and vegetation survey in July and August from 2006 to 2021 were measured. All the surveys were carried out in 2006, 2008, 2011, 2016, and 2021. In the forest inventory case, the diameter at breast height (DBH) and height of all the trees over 6 cm in DBH in each study plot were labeled and measured. With measured tree inventory data, aboveground biomass (AGB) was calculated using the allometric equation developed by Son et al. (2014) below. Then, AGC was calculated by multiplying the biomass by the carbon coefficient (0.5).

The AGB was the sum of three values obtained from the equations above.

In the 20 m × 40 m plots, we conducted an understory vegetation survey (shrub and herb layers), recording species presence and estimating cover scale (r, +, 1, 2, 3, 4, and 5) using the phytosociological method (Braun-Blanquet 1964) for SR. All the plant species in each plot were identified, following Korea National Arboretum (2003).

Relative benefits and trade-offs analyses

Standardization and RMSE were used to quantify the relative benefits and TOs suggested by Bradford and D’Amato (2012). A simple but effective method, RMSE quantifies the magnitude of TOs between two or more ecosystems. Data standardization is required to calculate TOs using the RMSE, removing the effect of different units (Bradford and D’Amato 2012). The equation is below.

Where ES_std indicates the standardized value of an ES, ranging from 0 to 1, ESosb, ESmax, and ESmin are the observed value, maximum value, and minimum value in the plot, respectively. The equation of RMSE is below.

Where ES_i is the standardized value of ith ES, and ES¯ is the expected value of the i number of ESs. Figure 2. is the diagram describing the TOs between two ESs. The degree of TOs was calculated using RMSE. The TO is the distance between a coordinate of a pair of ESs and the 1:1 line. As the distance is longer, RMSE and TOs are higher and greater, respectively. For instance, A, the coordinate point on the 1:1 line, represents zero TO. The level of TO of point B is greater than that of C. The degree of the TO of C is the same as D. The point’s relative position to the 1:1 line means which ES is more beneficial at the given condition. B is more beneficial for ES-2 (Bradford and D’Amato 2012; Lu et al. 2014).

Figure 2. Diagram of the trade-offs between two ecosystem services (ESs) (modified from the article of Bradford and D’Amato [Front Ecol Environ. 2012;10(4):210-6]). The blue lines indicate the magnitude of trade-offs between ESs.

Statistical analyses

One-way ANOVA (randomized block design) for AGC, SR, and relative benefits of AGC and SR were used to test the significance between treatments. Then, the post hoc test was performed with Tukey HSD (α = 0.05). All the statistical analyses were performed using R software.

The effects of thinning on DBH growth, AGC, and SR

There was a thinning effect on DBH growth after thinning (Table 2). In 2006, mean DBH among treatments did not differ from each other yet was significantly different (p < 0.05) from 2008 to 2021, representing LT = HT > Con. Also, the trend was observed in AGC. After thinning in 2007, the AGC across three treatments gradually increased until 2021 (Fig. 3A). AGC was not significantly different between treatments before thinning in 2006 yet showed a significant difference (p < 0.05) after thinning from 2008 to 2021, ranked following Con > LT > HT. Also, the block effect was observed from 2008 to 2021. Even though the block effect was statistically significant, it accounted for only a small proportion (around 0.20–0.27) of the variation, indicating that treatment was the major factor impacting AGC. In 2021, AGC (111.1 t C ha^-1) in the Con was about twice as high as that (60.4 t C ha^-1) in the HT, and AGC (82.0 t C ha^-1) in the LT in 2021 recovered to the level of AGC (78.6 t C ha^-1) before thinning in 2006.

Table 2 . The change of mean diameter at breast height.

Treatment	2006	2008	2011	2016	2021
LT	28.4 (± 2.3)ns	29.8 (± 2.3)a	31.4 (± 2.3)a	34.4 (± 2.8)a	36.1 (± 2.9)a
HT	26.5 (± 0.9)ns	29.0 (± 0.9)a	31.4 (± 0.9)a	34.6 (± 1.0)a	37.1 (± 0.6)a
Con	26.0 (± 2.1)ns	26.3 (± 1.2)b	27.1 (± 2.3)b	28.7 (± 2.6)b	29.9 (± 2.7)b

Values are presented as mean ± standard error.

LT: light thinning; HT: heavy thinning; Con: control.

Different letters by treatments indicate significantly different means (p < 0.05). ns indicates no significant difference.

Figure 3. The temporal changes of (A) aboveground carbon (AGC) storage and (B) species richness (SR) relying on thinning intensities. * and ∙ indicate significantly (p < 0.05) and marginally (p < 0.1) differences, respectively. LT: light thinning; HT: heavy thinning; Con: control.

SR in the LH in 2011 reached a peak of 94 species (± 10.4) and decreased to 41 species (± 8.4) in 2016 (Fig. 3). The dominant shrubs were Z. schinifolium and Clerodendrum trichotomum in 2011, and Z. schinifolium and Callicarpa japonica in 2016. In the herb layer O. undulatifolius and R. crataegifolius in 2011 and R. crataegifolius and Parthenocissus tricuspidate in 2016. In the HT, it was 88 species (± 4.3) in 2011 and 53 species (± 2.6) in 2016 (Fig. 3). The dominant shrubs were Z. schinifolium and A. elata (Miq.) in 2011, and Z. schinifolium and C. trichotomum in 2016. In the herb layer O. undulatifolius and Athyrium niponicum in 2011 and Actinidia arguta and R. crataegifolius in 2016. Generally, SR across treatments rapidly increased in 2011 and decreased in 2016. After thinning SR showed no differences between treatments before (2006) and after (2008) thinning; however, it was marginally higher in the LT (p < 0.1) than other treatments in 2011 (Fig. 3). The results of relative benefits had similar trends with obtained values (Fig. 4). AGC’s relative benefit significantly differed between treatments after thinning. However, the difference between obtained values and relative benefits was that the relative benefits of AGC in the Con and the LT were at peak in 2011 and 2008, respectively. Also, the relative benefit of SR in the HT increased over time, so it showed the highest value in 2021, unlike the obtained SR (Fig. 4).

Figure 4. The temporal changes of relative benefits of aboveground carbon (AGC) storage (A) and species richness (SR) (B) from 2006 to 2021. After thinning in 2006, AGC ranked following Con > LT > HT (p < 0.05), and SR on the LT was marginally higher in 2011 (p < 0.1). LT: light thinning; HT: heavy thinning; Con: control.

The effects of thinning on trade-offs

After thinning, the TOs across treatments decreased in 2008 and varied depending on thinning intensities (Figs. 5 and 6). Generally, the paired ESs in the LT relatively fell on the 1:1 line after thinning compared to other treatments, indicating that the TOs in the LT were relatively lower (overall mean: 0.21) than other treatments. On the other hand, the pairs of AGC and SR in the HT and Con deviated from the 1:1 line, representing that higher TOs occurred yet showed different patterns. HT’s TOs gradually increased along time courses. However, the TO on the Con reached a peak in 2011 and gradually decreased, which suggested that high TOs (overall mean: 0.45) and medium TOs (overall mean: 0.31) occurred in the HT and Con, respectively. Despite the TOs between the two ESs, LT, HT, and Con benefited SR, SR, and AGC, respectively (Fig. 6).

Figure 5. The pattern of trade-offs between aboveground carbon storage and species richness in the research period. RMSE: root mean squared error; LT: light thinning; HT: heavy thinning; Con: control.

Figure 6. Scatter-plot matrices of paired two ecosystem services from the thinning treatments. The distance between the circle and the 1:1 line indicates the degree of trade-offs (TOs). The green arrow represents the direction of the TO from 2006 to 2021. SR: species richness; AGC: aboveground carbon.

The study showed empirical evidence that TOs between two different ESs occurred, and the magnitude of TOs under thinning intensities varied. Also, we figured out that the temporal TO patterns among the two ESs changed depending on thinning intensities. The study demonstrated that the TOs between regulating and supporting services could either increase or decrease by thinning intensities compared to the Con in a relatively short time.

AGC is the highest in the Con among treatments over time after thinning due to the density of remaining trees, although the growth rate is lower in the Con than others. This result agrees with previous studies (Burton et al. 2013; Lee et al. 2005). Around thirteen years after thinning, AGC in the LT recovered from the disturbance to the level before thinning, suggesting the canopy gap in the LT also recovered or got closed to the level in 2006. The increase in SR after thinning may be due to an increase in light-preferred species related to the openness of the canopy in the overstory and increased early seral species on the disturbed forest floor (Ares et al. 2009; Fahey and Puettmann 2007; Kim et al. 2020). On the other hand, the decrease in SR across treatments in 2016 and 2021 may be considered a closed canopy and drought (Korea Meteorological Administration 2023), especially for LT. Another possible reason for the decrease in SR is probably that late seral species are competitively excluded by early seral species (Burton et al. 2013). Some light-preferred or early seral species invaded first after thinning, occupied certain areas, and became dominant (Kim et al. 2020). For example, Z. schinifolium at the shrub layer, R. crataegifolius, and A. arguta at the understory layer were dominant from 2011 to 2021 for the LT and HT. In addition, it is considered that the reduction in SR in the Con is attributed to limited growth areas and more closed canopy.

The thinning treatments had a great effect on the relative benefits and TOs. After thinning, among treatments, the HT and the LH were re-arranged on the upper left of the 1:1 line on the diagram (Fig. 6), where the relative benefit of SR was greater than that of AGC, although the degree of relative benefits between HT and LT was different; on the other hand, the Con was situated on the lower right of the diagonal line where the relative benefit of AGC was higher than that of SR. This explained that the lower thinning intensities were more beneficial to the accumulation of AGC storage, while the higher thinning intensities favored SR more. This is because of the stand density of remaining trees after thinning and the different gap sizes, microenvironmental changes such as light conditions, and soil moisture on thinned stands by thinning. The results agree with Zhu et al. (2018) regarding the relatively low thinning intensity favoring the AGC accumulation and previous studies on thinning intensities affecting SR (Li et al. 2020; Xu et al. 2020). This TO relationship between AGC and SR is also found in previous research (Mandal et al. 2013; Rana et al. 2017). The occurrence of TOs between ESs may be due to different rates of increase or decrease in each ES between paired ESs. This is because ESs have a tendency to change by time scales, yet the changes vary since ESs are the result of scale-dependent ecosystem processes (Deng et al. 2016; Lu et al. 2014; Rodríguez et al. 2006). Besides, supporting services such as biodiversity and nutrient cycle usually change more over time than other services (Millennium Ecosystem Assessment 2005). On the contrary, several papers discovered the synergy relationship between carbon stock and plant diversity in natural rainforests (Day et al. 2014) and subtropical eucalyptus forests (Zhou et al. 2017) probably due to site characteristic (high nutrient retention capacity) (Nilsson and Wardle 2005) and species characteristics (rapid growth and high carbon abortion capacity) (Zhou et al. 2017).

In the study, we analyzed the effect of thinning on ESs, relative benefits, and TOs between ESs. The response of each ES to thinning intensities and the degree of TO between the two ESs under different thinning intensities varied. We found a high relative benefit of AGC with a great TO on the Con but a high relative benefit of SR with a moderate TO on the HT. Moderate relative benefits of both AGC and SR with relatively lower TOs were observed on the LT. One objective of the study was to figure out an appropriate thinning method to enable an increase in synergy with a decrease in the TO between ESs. Unfortunately, there was no synergy among ESs after thinning in our study. However, we revealed that LT reduced the TOs between the two ESs with moderate increases in both AGC and SR. Thus, LT is a proper method in terms of the balance between enhancement and TO among ESs.

The authors are deeply grateful to the National Institute of Forest Science for funding during the research period.

LT: Light thinning

HT: Heavy thinning

Con: Control

ES: Ecosystem services

DBH: Diameter at breast height

AGB: Aboveground biomass

AGC: Aboveground carbon

SR: Species richness

TO: Trade-off

RMSE: Root mean squared error

ANOVA: Analysis of variation

KL did conceptualization, data curation, investigation, data analysis, and writing-original draft. SJY and MK did data curation and investigation. HMY did conceptualization, investigation, and writing-review and editing. ARK did investigation, supervision, and writing-review and editing. All the authors approved the manuscript.

This work was supported by a grant (Project no.: FE0100-2021-01) from the National Institute of Forest Science, Republic of Korea, Republic of Korea.

The datasets are available from the corresponding author upon reasonable request.

Not applicable.

Not applicable.

The authors declare that they have no competing interests.