Wildfire and timber harvesting drive long-term modifications in soil properties of Southern Patagonian Nothofagus forests

1 Introduction

Species composition of temperate natural forests can influence soil physicochemical properties, nutrient concentrations, and nutrient stocks (Toro Manríquez et al., 2019; Han et al., 2022; Xiang et al., 2025). Tree species can differ in their effects on soil by many mechanisms, including rates of chemical weathering of parental material through changes in pH and composition of the complex organic acids in soil solutions, mostly in topsoil or near the roots (Courchesne and Gobran, 1997; Augusto et al., 2000, 2002), e.g., litter input contribution varies with canopy tree species because they have different decomposition rates and nutrient contents in the tissues (Augusto et al., 2002; Vesterdal et al., 2008; Yang et al., 2022).

Soil exhibits complex structures that support and regulate many ecological functions in natural ecosystems (Luo et al., 2022; Telo da Gama, 2023; Wang et al., 2024). Sampling efforts in monitoring usually include few variables in the analyses, e.g., carbon stock, soil acidity, and nutrient concentrations (Alauzis et al., 2004; Han et al., 2022; Nosetto et al., 2024; Barbosa Pimentel et al., 2025), however, other components (e.g., texture, roots, woody debris) showed greater performance as indicators in ecological studies of temperate forests (Gatica-Saavedra et al., 2023; Aravena Acuña et al., 2023; Chaves et al., 2023; Zhang et al., 2025). Besides, understory plants greatly influence topsoil, modifying their properties and generating positive and negative synergies with trees (Romanya et al., 2005; Fontaine et al., 2023; Xiong et al., 2024; Haritika and Negi, 2025). For this, the inclusion of more diverse soil variables can increase the understanding of changes produced by human impacts, e.g., wildfires or harvesting (Longo et al., 2011; Toro Manríquez et al., 2019; Martínez Pastur et al., 2021, 2022, 2024b; Soler et al., 2022; Toledo et al., 2024, 2025; Rosas et al., 2024).

Nothofagus pumilio (Poepp. and Endl.) Krasser and N. antarctica (G. Forst.) Oerst. are broad-leaf deciduous tree species growing in the southernmost forests of the world (Gutiérrez, 1994; Frangi et al., 2015). They present regional commercial values for harvesting (e.g., timber and firewood) and ranching activities (e.g., silvopastoral systems; Gea-Izquierdo et al., 2004; Martínez Pastur et al., 2022; Peri et al., 2022). These trees are medium shade-tolerant species and can grow in monospecific or mixed stands (Alauzis et al., 2004; Toro Manríquez et al., 2019, 2022). Natural regeneration occurs abundantly after disturbances, creating even- or uneven-aged forest structures, e.g., seedlings originated from seeds are mostly present in N. pumilio, while root sprouts are more frequent in N. antarctica (Soler et al., 2022; Martínez Pastur et al., 2023a; Rodríguez Souilla et al., 2023). The succession dynamics are predictable and mostly depend on periodic fine-scale disturbances (e.g., windthrow of individuals or small groups of trees; Pollmann and Veblen, 2004; Soto et al., 2019; Martínez Pastur et al., 2024a). Besides, soil characteristics in Nothofagus forests reflect the influence of the dominant pedogenic processes that differ between species (Frangi et al., 2005a; Romanya et al., 2005; Toro Manríquez et al., 2019). Thus, soil properties may contribute to the local occurrence of Nothofagus species under analysis (Frangi et al., 2005a), e.g., N. pumilio prefers deeper soils in the uplands while N. antarctica can live in diverse habitat conditions.

Human associated impacts on natural ecosystems have been identified as the main drivers of change in temperate forests around the world (Bače et al., 2017; Bousfield et al., 2023; Cerioni et al., 2024), as well as in southern Patagonia (Alauzis et al., 2004; Klein et al., 2008; Soler et al., 2022). Among them, wildfire and harvesting are identified as the most important drivers. These impacts modify forest structure, forest persistence, species assemblage, ecosystem functions (Lindenmayer et al., 2016; Martínez Pastur et al., 2021; Chaves et al., 2023), and soil properties (Longo et al., 2011; Toledo et al., 2024, 2025). However, soils are usually not included in long-term studies due to the complexity and time-consuming sampling effort, both during field and laboratory activities (Longo et al., 2011; Gatica-Saavedra et al., 2023; Loria et al., 2024). For this, the studies that do not include the below-ground forest properties leave aside relevant variables in the analyses (Martínez Pastur et al., 2022), e.g., soils are identified as one of the most important components of natural ecosystems, influencing productivity and supporting vital ecological processes (Peri et al., 2019; Nosetto et al., 2024; Barbosa Pimentel et al., 2025). The impacts include alterations on: (i) biophysical characteristics, (ii) nutrient content and biochemical cycles, (iii) ecosystem function, (iv) successional dynamics, and (v) community composition and soil species richness (Klein et al., 2008; Lencinas et al., 2019; Sokol et al., 2022; Rosas et al., 2024; Xiang et al., 2025; Zhang et al., 2025; Ma et al., 2026). In Patagonian forests, most research of anthropogenic impact studies have focused on above-ground components (Böswald et al., 2002; Huertas Herrera et al., 2018, 2023; Martínez Pastur et al., 2021; Arroyo-Vargas et al., 2024), and few have investigated the effects on soil properties (Klein et al., 2008; Gargaglione et al., 2022; Rosas et al., 2024; Arguiano et al., 2025).

Ecosystems are dynamic systems with complex responses to environmental variation (Falk et al., 2022), where ecosystem succession culminates in stabilized ecosystems which are characterized by closed biogeochemical cycles and high capacity for nutrient retention (Odum, 1969), directly linked to the net ecosystem productivity (Vitousek and Reiners, 1975; Huygens et al., 2007). One consequence of the human-related impacts (e.g., wildfire and harvesting) is the increasing acute and chronic stress from these vegetation and environmental changes, where natural forests can be pushed outside of their historical ranges of variation into alternative ecological states, e.g., conversion of forest to open forests or openlands such as grasslands (Batllori et al., 2019; Landesmann et al., 2021; Falk et al., 2022). These new states can be transitory and represent successional stages that may ultimately lead back toward pre-disturbance conditions; however, these states can persist and potentially self-reinforce on time (Frangi et al., 2015), especially under the cumulative effects of different impacts and presence of exotic species (Rastogi et al., 2023; Peller and Altermatt, 2024; Zandalinas et al., 2024). In this context, understanding the ecological processes related to forest resilience is crucial for sustainable management or ecosystem restoration strategies (Rist and Moen, 2013; Ren and Coffman, 2023; Triviño et al., 2023). Resilience can be defined as the ability of forests to recover or adapt to natural and anthropogenic disturbances, resulting from the expression of multiple mechanisms operating at levels of organism, population, and community interactions (Falk et al., 2022). Patagonian forests have experienced an increasing rate of disturbances during the last decades (Peri et al., 2022), which has sparked interest in assessing post-disturbance forest recovery (Martínez Pastur et al., 2025a) and the identification of forest recovery states (Peri et al., 2017; Martínez Pastur et al., 2021, 2022, 2026). Recovery from forest disturbances is a pivotal metric of forest resilience (Mandl et al., 2024), where forest recovery can be measured in different dimensions (e.g., forest structure, composition, and function), including different ecosystem components (e.g., trees, soils, plants, animals; Falk et al., 2022; Smith et al., 2022; Mandl et al., 2024). The inclusion of different dimensions and components allowed to determine the positive and negative synergies that influence forest resilience (Nikinmaa et al., 2020; Steel et al., 2021; Liu et al., 2025). This is a major gap to define new management proposals in the long-term (Rist and Moen, 2013; Triviño et al., 2023; Cerioni et al., 2024), and particularly in Patagonian forests (Martínez Pastur et al., 2026). In this context, the objective was to determine resilience of forest soils (Nothofagus pumilio, N. antarctica, and mixed) after long-term natural recovery from human impacts (wildfire and harvesting). We wanted to answer: (i) Can original soil properties are maintained after long-term impacts? (ii) Can impact type generate negative synergies in soil properties after long-term impacts? (iii) It is possible to identify recovery fails for the studied impacts that can be used to improve soil conservation and management proposals?

2 Methods

2.1 Study area

Sampling was conducted at La Catalana ranch (−54.1946, −67.3224) located in Tierra del Fuego (Argentina; Figure 1A). The ranch manages 10,339 ha for cattle breeding (4,290 ha of natural grasslands, 5,926 ha of native forests, and 123 ha of degraded forests; Martínez Pastur et al., 2026). The climate is characterized by short and cool summers, and long snowy winters. Precipitation in the study area ranged from 370 to 400 mm yr⁻¹, and mean annual temperature is 5.2 ± 0.1 °C (1.6 °C in the coldest and 9.6 °C in the warmest months), and where frost may occur at any time of the year. Soils are classified as podzols according to IUSS Working Group WRB (2014). The landscape occupied by these ecosystems is predominantly of glacial origin, composed of loess and alluvial deposits in the foothill region. Soils have a silty loam texture and slightly acidic pH, and are characterized by high permeability, good drainage, and high-water retention capacity (Martínez Pastur et al., 2022; Toledo et al., 2025). Forests include pure Nothofagus antarctica (75%) and N. pumilio (20%) stands, and mixed ones of both species (5%; 70% of basal area of N. antarctica and 30% basal area of N. pumilio). We selected one paddock inside the ranch with a well-known history of impacts, covering 300 ha used for harvesting of timber and firewood during 1960–1966, when heavy wildfire affected the area. After these impacts, the paddock was mainly used for sheep grazing for 50 years, followed by land abandonment, and conversion to cattle breeding during the last decade (Martínez Pastur et al., 2026). In this paper we will focus on the impacts of harvesting operations with or without wildfire damage; besides, we have not considered the potential impact of livestock or invasive plants in the understory. Livestock move freely across the study area (one paddock) and occur at a low density (one animal every 8–10 ha), and many invasive plants modify the understory assemblage and quality of the forage (e.g., Hieracium pilosella) and can leading to changes on the soil’s properties in the long-term (Martínez Pastur et al., 2022).

2.2 Theoretical framework of the study

We described the theoretical framework in a previous study (Martínez Pastur et al., 2026) based on the scientific background of Frangi et al. (2015) and Peri et al. (2022) (Figure 2). The proposal follows a natural ecological succession model from openlands to closed canopy N. pumilio forests, including: (i) Grasslands (G) dominated by Festuca gracillima intermixed with small shrubs of Empetrum rubrum, which increase their plant richness in wet areas (e.g., Juncus scheuchzerioides and Carex curta). (ii) Open N. antarctica forests (ONA) with greater presence of dicot plants. (iii) Close N. antarctica forests (CNA) with greater plant richness in the understory that includes more ferns and mosses presence. (iv) Mixed forests of N. antarctica and N. pumilio (MIX) representing an intermediate state of both pure forest types. (v) On time N. pumilio trees grow faster than N. antarctica trees, leading to pure N. pumilio stands (NP). The described succession model can follow a natural inverse pathway when the environmental conditions change (Toro Manríquez et al., 2022). These natural vegetation types can be modified due to wildfire and harvesting (timber and firewood), where we identified the following treatments: (vi) Harvested N. antarctica forests (remnant basal area BA > 50%) in open and closed forests (ONA-H and CNA-H). (vii) Light (NP-LH, remnant BA > 50%) and heavy harvested (NP-HH, remnant BA < 50%) N. pumilio forests. (viii) Harvested forests affected by wildfire (ONA-H+F, NP-LH+F, NP-HH+F), including harvested mixed (N. antarctica+N. pumilio) forests (MIX-H+F). The non-impacted status and disturbance intensity were determined during the field work including the presence of: (i) stumps, (ii) logs or remnant coarse-woody-debris with wildfire or harvesting scars, and (iii) damage in standing trees due to harvesting or scars due to wildfire. This survey was possible due to the latitude where these forest growth and low decomposition rates that allowed the conservation of the woody material in the stands (Frangi et al., 1997, 2005a). The wildfire impact alone was not included due to the low presence in the landscape, and the analyses only consider the additive effect of wildfire in the harvested stands.

2.3 Sampling design and data taking

A total of 130 plots were measured in the study site, using a systematic grid of 150 m × 150 m across forestlands and openlands (Figure 1B). We use this sampling design to capture the variability of the wildfire event, where we sampled stands or openlands (replication unit) in all the affected and surrounded areas. With this design we obtained an unbalanced number of replicas (stands or openland areas) according to the abundance of each defined treatment in the landscape. A total of 41 plots corresponded to non-impacted environments: N. pumilio forests (NP, 6 plots), mixed N. pumilio and N. antarctica forests (MIX, 5 plots), closed N. antarctica forests (CNA, 9 plots), open N. antarctica forests (ONA, 16 plots), and grasslands (G, 5 plots). The other plots corresponded to different combinations of forest types and impacts: N. pumilio forests with light harvesting (NP-LH, 9 plots), N. pumilio forests with heavy harvesting (NP-HH, 5 plots), N. pumilio forests with light harvesting affected by wildfires (NP-LH+F, 7 plots), N. pumilio forests with heavy harvesting affected by wildfires (NP-HH+F, 26 plots), mixed forests with harvesting affected by wildfires (MIX-H+F, 6 plots), close N. antarctica forests with light harvesting (CNA-H, 10 plots), open N. antarctica forests with light harvesting (ONA-H, 18 plots), and open N. antarctica forests with light harvesting affected by wildfires (ONA-H+F, 8 plots).

At each plot center defined in the grid (latitude and longitude) we established a 25 m transect following sampling design proposed by Rosas et al. (2024). Forest structure was characterized at each plot, and outputs were previously presented elsewhere (Martínez Pastur et al., 2026). Naked soil cover (NS, %) and understory cover (UC, %) were obtained using the point-intercept method (Levy and Madden, 1933) with 50 intercept points (every 0.5 m) along the transects located at each plot. We also characterized the 30 cm soil layer at each plot using a composite of four samples using a field borer of known volume after removing the litter, including two locations along the transect (between 5–6 m and between 19–20 m) and two depths (0–15 cm and 15–30 cm layers). The subsamples’ identity was maintained for the analyses, but average values were used for the analyses. Samples were weighed before and after drying in laboratory oven (70 °C) until constant weight was reached, obtaining soil moisture (SM, %). We use 70 °C to avoid combustion or alteration of organic matter (Martinez et al., 2017) following recommendations of PEBANPA network and Field Manual for the Determination of Organic Carbon in Forest Soils (Ministerio de Ambiente y Desarrollo Sostenible (MayDS), 2021). The sieve soil through a 2 mm mesh separates the fine earth fraction from larger debris like stones, gravel, roots, and woody debris (e.g., dead roots, or dead wood that arrive in soil across the years). These soil weights allowed us to calculate effective soil bulk density (SD, gr cm³). Larger organic components were manually classified to estimate root soil content (ROOT, t ha⁻¹ in layer 0–30 cm) and dead wood debris soil content (WOOD, t ha⁻¹ in layer 0–30 cm).

For chemical analyses, we pooled the individual soil subsamples (locations into the transect) into one combined sample per plot but maintaining the layer depth identity (0–15 cm and 15–30 cm layers) for some variables. Soil acidity (pH) was estimated in a suspension of a soil/water ratio of 1:2.5 (Bao, 2000). Soil organic carbon was determined using dry combustion analyses using muffle furnace at 500 °C for 24 h and modeled following Martínez Pastur et al. (2022) using equations for each vegetation type based on organic matter and C determined with dry combustion (LECO CR-12, LECO Corporation, USA). With these analyses, we obtained: (i) upper soil organic carbon (C-UP, % in layer 0–15 cm), (ii) down soil organic carbon (C-DOWN, % in layer 15–30 cm), (iii) ratio of soil organic carbon between layers (RATIO, %) comparing C-UP over C-DOWN for each plot, and (iv) total soil organic carbon as the average of the four subsamples (C, % in layer 0–30 cm). Besides, soil nitrogen (N, % in layer 0–30 cm) was determined using a semi-micro Kjeldahl method (Carter and Gregorich, 2007), and soil available phosphorous (P, ppm in layer 0–30 cm) was determined using the Bray and Kurtz method (Bray and Kurtz, 1945) as the average of the four subsamples. Soil density and nutrient concentrations of each plot allowed us to calculate (i) soil organic carbon content (SOC, t ha⁻¹ in layer 0–30 cm), (ii) soil nitrogen content (SN, t ha⁻¹ in layer 0–30 cm), and (iii) soil available phosphorous content (SP, t ha⁻¹ in layer 0–30 cm).

2.4 Statistical analyses

To evaluate differences among vegetation types without impacts (NP, MIX, CNA, ONA, G) and different impact combinations (NP-LH, NP-HH, NP-LH+F, NP-HH+F, MIX-H+F, CNA-H, ONA-H, ONA-H+F), Generalized Linear Models (GLMs) were fitted using R-software version 4.3.3 (R Core Team, 2024). These models were applied to analyze: (i) soil physicochemical properties (SD, SM, pH, C-UP, C-DOWN, RATIO), (ii) nutrient concentrations and stocks (C, SOC, N, SN, P, SP), and (iii) roots, woody debris in soils (ROOT, WOOD), and understory covers (NS, UC). To properly accommodate strongly unbalanced sample sizes across treatments and handle non-normal, heteroscedastic data distributions, GLMs were specified with appropriate error distributions and link functions, e.g., Gaussian distributions with identity link functions were used for symmetric continuous variables (SD, SM, pH, SOC, SN, NS, UC), while Gamma distributions with log link functions were implemented for positively skewed continuous variables, ratios, and nutrient concentrations (C-UP, C-DOWN, RATIO, C, N, P, SP, ROOT, WOOD). Model assumptions, including residual uniformity/normality, variance dispersion (homoscedasticity), and outlier detection were validated using simulated residuals with DHARMa package (Hartig, 2022). Overall significance of environment factor was evaluated using Type-II Likelihood Ratio Chi-square test (χ², df = 12) at p < 0.05 via car package (Fox and Weisberg, 2019). Post-hoc pairwise comparisons among all environment variables were performed using Estimated Marginal Means (EMMs) with Tukey test (p < 0.05) adjusted using emmeans (Lenth, 2024), multcomp (Hothorn et al., 2008), and multcompView packages (Graves et al., 2019). Furthermore, planned pairwise contrasts were conducted to specifically compare each disturbed treatment against its corresponding unimpacted reference forest type, applying False Discovery Rate (FDR) corrections to account for multiple testing (Benjamini and Hochberg, 1995).

Principal Component Analyses (PCA) were conducted in R-software using prcomp function (stats package). To eliminate conceptual redundancy and multicollinearity (e.g., among overlapping carbon fractions and nutrient concentrations vs. stocks), variables were pre-screened using Pearson correlation coefficients (r > 0.70; Annex 1). A curated subset of 10 non-redundant variables was retained: soil physicochemical properties (SD, SM, pH), nutrient concentration and stocks (RATIO, SOC, SN, SP), and soil components and understory cover (ROOT, WOOD, UC). All variables were centered and standardized to unit variance (scale = TRUE) based on the ordination of correlation matrix. Separate exploratory PCAs were performed for: (i) unimpacted forests and grasslands (NP, CNA, ONA, MIX, G), (ii) N. pumilio stands, (iii) N. antarctica stands, and (iv) mixed stands. Convex hulls were computed using the chull function in R-software for visual group delimitation, while statistical inferences on disturbance effects were derived strictly from the GLM frameworks.

We developed one soil index (SI) to compare the non-impacted and impacted vegetation types based on soil properties and their associated components (see variables in Tables 1–3). Each variable was standardized (0–1) despite the studied treatments, where 1 is associated with climax vegetation type without impacts (NP). Three variables were converted (from 1 to 0) due to their increasing was inversely proportional to those variables found in the climax stage without impacts (NP), e.g. (i) an increase of SD, due to it was considered negative for recovery (see Fajardo and Gundale, 2015; Nosetto et al., 2024), (ii) RATIO was preferred as a minimum (no difference between layers) than maximum differences between topsoil and lower layer that can be associated with soil erosion (see Chaves et al., 2023; Aravena Acuña et al., 2023; Gatica-Saavedra et al., 2023), and (iii) UC was associated with natural regeneration failure (see Huertas Herrera et al., 2018; Martínez Pastur et al., 2021, 2026). After these conversions, all the variables were averaged without any weight. These outputs were standardized again to obtain the final indexes that varied between 0 and 1, where 0 showed greater dissimilarities and 1the most similar conditions with the climax stage without impacts (NP) found in the study area. The soil index (SI) were compared with forest structure index (FSI) described in Martínez Pastur et al. (2026), which included (i) dominant height (m), (ii) overstory crown cover (%), (iii) relative leaf area index (LAI), (iv) basal area (m² ha⁻¹), and (v) total volume of live trees (TV, m³ ha⁻¹). These relationships (average ± standard error of the indexes SI and FSI) were used to detect similarities and dispersion patterns among treatments and to characterize recovery values among treatments. Here, we characterize recovery according to its dispersion compared to non-impacted treatments.

3 Results

Most of the soil properties presented significant differences among non-impacted and impacted vegetation types, including soil density, soil acidity, and ratio between carbon concentration of upper (0–15 cm) and down (15–30 cm) layers (Table 1). Soil density greatly varied among non-impacted vegetation types, where mixed forests presented the lowest value (0.60 g cm³) and grasslands presented the highest value (0.95 g cm³). Impacted environments increased soil density (up to 40%) with different magnitudes and trends when they were compared to their corresponding non-impacted vegetation type. Soil acidity also changed among non-impacted vegetation types, presenting a natural gradient from open environments (e.g., 4.60 in grasslands) to close N. pumilio forests (e.g., 5.71). Soil acidity increased after the impacts in N. pumilio, mixed forests, and closed N. antarctica (−1% to −8%) and decreased in open N. antarctica impacted forests (2%–13%). Soil organic carbon contents at both studied layers did not change across non-impacted and impacted vegetation types; however, they tended to decrease (−2% to −44%). Finally, ratio between carbon concentration of upper (0–15 cm) and down layers (15–30 cm) changed among non-impacted vegetation types, presenting higher differences between N. pumilio and mixed forests (1.98–2.38) and other non-impacted treatments (1.33–1.46). This ratio was modified depending on impact treatments, e.g., harvesting in N. pumilio increased the values (5%–40%), while harvesting followed by wildfires decreased the values (−19% to −21%) after-60-years-impacts. However, in mixed and N. antarctica forests, this ratio tends to increase (5%–19%) compared to non-impacted treatments.

Table 1. Generalized linear model (GLM) analysis of soil properties values in non-impacted (NP = Nothofagus pumilio forests, MIX = mixed N. pumilio and N. antarctica forests, CNA = closed N. antarctica forests, ONA = open N. antarctica forests, G = grasslands) and impacted environments (H = harvesting, LH = light harvesting, HH = heavy harvesting, F = harvesting followed by wildfire) including soil bulk density (SD, gr cm³), soil moisture (SM, %), soil acidity (pH), upper soil organic carbon (C-UP, % in layer 0–15 cm), down soil organic carbon (C-down, % in layer 15–30 cm), and ratio between soil organic carbon layers (RATIO).

Values between brackets presented the changes in the values comparing impacted and non-impacted environments of the same type. See statistical details on Annexes 2, 3.

χ² = Likelihood Ratio Chi-square test (df = 12) from Generalized Linear Models (GLMs). Different letters within the same column show significant differences among environments using Estimated Marginal Means with Tukey adjustments (p <0.05).

Soil nutrients presented different behavior when non-impacted and impacted vegetation types were compared (Table 2). Soil organic carbon concentration (%) and content (t ha⁻¹) did not present differences 60-years-after-harvesting. Nevertheless, nitrogen and phosphorous, both in concentration (%) and content (t ha⁻¹), presented significant differences among treatments. Soil nitrogen concentration (%) in non-impacted vegetation types was greater in mixed and N. antarctica forests (0.38% to 0.43%) compared to N. pumilio forests and grasslands (0.28% and 0.29%, respectively). Soil nitrogen content (t ha⁻¹) presented a different trend in non-impacted vegetation types, where N. pumilio forests had the lowest values (5.55 t ha⁻¹) and N. antarctica had the greatest values (8.25–8.66 t ha⁻¹). These soil nitrogen values significantly changed in impacted forests. Harvested N. pumilio forests maintained or decreased soil nitrogen concentration (0% to −21%), while harvested with wildfire N. pumilio forests increased these values (4% to 11%). The other impacted forests (mixed and N. antarctica forests) decreased their soil nitrogen concentration (−10% to −21%). Soil nitrogen content (t ha⁻¹) slightly decreased in light harvested N. pumilio forests (−1%) while increased in other treatments (4% to 30% in heavy harvested N. pumilio forests, mixed and N. antarctica close forests with harvesting and harvesting followed by wildfire); however, in impacted open N. antarctica forests, the soil nitrogen content was reduced (−12% to −16%). Harvesting intensity can be the major response variable in N. pumilio forests, e.g. −1% in NP-LH compared to +25% in NP-HH, and +16% in NP-LH+F compared to +30% in NP-HH+F. Soil available phosphorous was maximum in N. pumilio forests (87.8 ppm and 0.186 t ha⁻¹) while mixed forests presented intermediate values (44.7 ppm and 0.081 t ha⁻¹) compared to N. antarctica forests and grasslands (10.8 to 37.9 ppm, and 0.029 to 0.074 t ha⁻¹). Soil available phosphorus concentration diminished in most of the impacted forests (−8% to −53% in N. pumilio, mixed and N. antarctica close forests) but increased in impacted N. antarctica open forests (5% to 30%). Finally, soil phosphorous content (t ha⁻¹) presented a dissimilar behavior among treatments, where some values presented a remarkable reduction, e.g., most impacted N. pumilio forests (−46% in harvested followed by wildfire stands) and impacted N. antarctica open forests (39% in harvested alone and harvested followed by wildfire).

Table 2. Generalized linear model (GLM) analysis of soil nutrients in non-impacted (NP = Nothofagus pumilio forests, MIX = mixed N. pumilio and N. antarctica forests, CNA = closed N. antarctica forests, ONA = open N. antarctica forests, G = grasslands) and impacted environments (H = harvesting, LH = light harvesting, HH = heavy harvesting, F = harvesting followed by wildfire) including soil organic carbon (C, % in layer 0–30 cm), soil organic carbon content (SOC, t ha⁻¹ in layer 0–30 cm), soil nitrogen (N, % in layer 0–30 cm), soil nitrogen content (SN, t ha⁻¹ in layer 0-30 cm), soil phosphorous (P, ppm in layer 0–30 cm), and soil phosphorous content (SP, t ha⁻¹ in layer 0–30 cm).

Values between brackets presented the changes in the values comparing impacted and non-impacted environments of the same type. See statistical details on Annexes 2, 3.

χ² = Likelihood Ratio Chi-square test (df = 12) from Generalized Linear Models (GLMs). Different letters within the same column show significant differences aong environments using Estimated Marginal Means with Tukey adjustments (p <0.05).

Root soil content (t ha⁻¹) did not vary among vegetation types (13.9–34.0 t ha⁻¹ in non-impacted and 18.4–30.3 t ha⁻¹ in impacted treatments; Table 3). Besides, wood debris soil content was higher in N. pumilio forests (22.4 t ha⁻¹) compared to mixed (16.3 t ha⁻¹), N. antarctica forests (6.3–9.0 t ha⁻¹), and grasslands (0.7 t ha⁻¹). Harvested N. pumilio forests decreased their wood debris soil content (−17% to −55%) while harvested followed by wildfire stands maintained or decreased these contents (1% to −61%). Mixed forests and harvested N. antarctica open forests also decreased the wood debris soil content (−25% to −28%), while other N. antarctica forests increased their values (4%–106%). Naked soil cover varied among non-impacted treatments, with minimum values in N. antarctica open forests (11.4%) and maximum values in grasslands (30.0%). Naked soil cover increased in harvested N. pumilio forests and stands with low harvesting followed by wildfire (8%–44%); however, these values decreased in heavy harvested followed by wildfire N. pumilio stands (−71%). This decrease in naked soil cover was also found in other impacted treatments (−31% to −58% in mixed and N. antarctica forests). Finally, understory cover followed an inverse behavior compared to naked soil cover; however, their magnitude, due to the number of plant layers, changed with vegetation type. Nothofagus antarctica forests presented the greatest understory covers (90.9%–97.0%) compared to other treatments (73.3%–74.4%). Understory cover decreased in harvested N. pumilio forests and stands with low harvesting followed by wildfire (−6% to −18%); however, these values increased in heavy harvested followed by wildfire N. pumilio stands (37%). Finally, this increase in understory cover was also found in other impacted treatments (9%–24% in mixed and N. antarctica forests).

Table 3. Generalized linear model (GLM) analysis of roots, woody debris and understory cover in non-impacted (NP = Nothofagus pumilio forests, MIX = mixed N. pumilio and N. antarctica forests, CNA = closed N. antarctica forests, ONA = open N. antarctica forests, G = grasslands) and impacted environments (H = harvesting, LH = light harvesting, HH = heavy harvesting, F = harvesting followed by wildfire) including root soil content (ROOT, t ha⁻¹ in layer 0–30 cm), wood debris soil content (WOOD, t ha⁻¹ in layer 0–30 cm), naked soil cover (NS, %), and understory cover (UC, %).

Values between brackets presented the changes in the values comparing impacted and non-impacted environments of the same type. See statistical details on Annexes 2, 3.

χ² = Likelihood Ratio Chi-square test (df = 12) from Generalized Linear Models (GLMs). Different letters within the same column show significant differences among environments using Estimated Marginal Means with Tukey adjustments (p <0.05).

Planned pairwise contrasts between impacted and their corresponding non-impacted environment showed distinct legacy patterns and vulnerability thresholds across 60-year post-disturbance treatments (Annexes 2, 3). In N. pumilio forests, heavy harvesting followed by wildfire (NP-HH+F) induced persistent stand modifications, characterized by a reduction of available soil phosphorus, an increase in understory cover, and reduction in naked soil cover compared to non-impacted stands (NP). In contrast, N. pumilio stands impacted with light harvesting followed by wildfire (NP-LH+F) not present differences across the studied soil properties compared to non-impacted forests, showed higher long-term soil resilience compared to higher harvesting intensity (NP-HH+F). Impacted mixed forests (MIX-H+F) showed persistent soil degradation, evidenced by a significantly higher soil density compared to non-impacted stands. Open N. antarctica forests affected by harvesting followed by wildfire (ONA-H+F) significantly differed from non-impacted stands, showing higher soil pH and greater soil root content. Conversely, soil carbon (C and SOC); its distribution across profile (C-UP, C-DOWN, RATIO); and soil moisture not showed differences between impacted and non-impacted stands after 60-years.

Principal Component Analysis (PCA) before impact in the studied environments (Figure 3A) reinforced the theoretical framework, revealing that environmental gradients were driven by soil properties. Non-impacted N. pumilio stands (NP) distinctly segregated on the positive PC2 versus negative PC1 quadrant, influenced by higher SP, pH, RATIO, and WOOD. Conversely, grasslands (G) and open N. antarctica stands (ONA) occupied the opposite ordination space, strongly influenced by SD, UC, and SN, while closed N. antarctica (CNA) and mixed stands (MIX) showed an intermediate position associated with SM and ROOT. The second ordination for N. pumilio stands (Figure 3B) showed a remarkable expansion of multivariate soil variables dispersion as disturbance intensity increased. Non-impacted and light harvested stands with or without impact of wildfire (NP-LH and NP-LH+F) occupied a more constrained space associated with a better-preserved carbon ratio, SP, and soil wood debris. In contrast, heavy harvested following by wildfire stands (NP-HH+F) expanded the multivariate polygon, positively associated to PC1 and negatively associated to PC2 toward increased UC, SM, and SN, or toward increasing soil compaction (SD), confirming a profound alteration of the original soil characteristics. The third ordination analysis for N. antarctica stands (Figure 3C) showed substantial polygon overlapping among non-impacted (CNA, ONA) and impacted stands (CNA-H, ONA-H, ONA-H+F), showing the resilient capacity of this tree species. The main differentiation occurred between closed canopies (associated with higher SOC, SN, and WOOD on the bottom-left) and open or harvested followed by wildfire stands with higher ROOT, UC, and SD. Finally, the ordination for mixed forests (Figure 3D) clearly separated non-impacted stands (MIX on the left, associated with SM, ROOT, and SN) from impacted stands (MIX-H+F on the right), which showed increased soil compaction (SD) and altered C distribution across the soil profile (RATIO).

Non-impacted and impacted vegetation types showed an increasing value of both developed indexes (FSI vs. SI, average ± standard error; Figure 4). Extreme values for both indexes (FSI and SI) were occupied for early (grasslands) and climax stages (non-impacted N. pumilio forests). The impacted forests decreased their FSI and SI values compared to non-impacted forests, but the magnitude of these changes depends on index and forest type, in example: (i) Light harvesting (NP-LH) presented lower changes compared to heavy harvesting (NP-HH) on forest structure index (FSI) in N. pumilio stands, (−5.4% and −16.2% compared to non-impacted forests), but was more affected in soil index (SI; −28.4% and −18.2% respectively). (ii) This trend changed when harvesting was combined with wildfire in N. pumilio forests (NP-LH+F and NP-HH+F), where values tend to increase in heavy harvesting treatments (−27.5% to −71.3% in FSI, and −9.7% to −36.6% in SI). (iii) Mixed forests (MIX) also changed their trend in both indexes (−41.9% in FSI and −33.6% in SI). Finally, (iv) N. antarctica stands presented a different pattern of change compared to those forests with N. pumilio trees, where FSI increased the values (−1.6% to −48.1%) but where SI decreased their values (−15.2% to −29.4%). These comparisons allowed us to determine the main recovery characteristics of the studied vegetation types: (i) Heavy impacted (wildfire+harvesting) N. pumilio forests (NP-HH+F) are separated from the other studied treatments for the same forest type. (ii) Open N. antarctica impacted stands (ONA-H and ONA-H+F) presented less recovery than close stands (CNA-H). And (iii) impacted mixed forests also presented a lower recovery after 60-years-after-impacts. These combined outputs of forest structure and soil variables confirmed the theoretical framework presented in methods.

4 Discussion

The study site presents a gradient of natural vegetation types, which is representative of the ecotone zone in Tierra del Fuego (Argentina), varying from openlands (grasslands) to closed N. pumilio mature forests (Rodríguez Souilla et al., 2026a). The origin of soils is glacio–fluvial (Martínez Pastur et al., 2023b), but the soil development differed during the last millennia and was related to vegetation that influenced over main characteristics (Gutiérrez, 1994; Romanya et al., 2005; Toro Manríquez et al., 2019). Our study is limited to the topsoil layer (0–30 cm), which is influenced by roots and microbial activity (Ross and Täte, 1993; Martínez Pastur et al., 2024b), as well as human impacts (Toledo et al., 2024, 2025). Forestlands and openlands present different soil properties (Tables 1–3), mainly due to leaf litter and woody-debris inputs, modifying soil characteristics on time (Jia et al., 2022; Rodríguez Souilla et al., 2023). Besides, microclimate also changes across the studied vegetation gradient, where extreme climate values occur in openlands (e.g., maximum and minimum temperatures, as well as the quantity of rainfall that reaches the soil) during short summers and long winters, e.g., at this latitude the growing season is about 4 months (Martínez Pastur et al., 2023a; Rodríguez Souilla et al., 2026b). The difference in leaf litter inputs and microclimate influences decomposition rates and microbial communities, e.g., soil acidity, nutrient contents, and dead wood soil content (Frangi et al., 1997, 2005b; Bahamonde et al., 2012; Soto, 2025). For this, it is necessary to compare treatment by treatment the impacts (harvesting with or without wildfire influence) and the recovery after-60-years-impacts.

The greatest impacts after disturbances occurred immediately after the impact events due to changes in the forest structure, and due to overstory dieback in the successive years, e.g., many medium- and long-term studies suggest 7–9 years for the greater differences compared to non-impacted forests, followed by a fast recovery in forest structure and functions of soils and ecosystem components (Chaves et al., 2024; Fachin et al., 2024; Toledo et al., 2024, 2025). In our study, we analyze the recovery after-60-years-impacts without considering the pathways followed by the different treatments in between. In impacted treatments, soil density generally increased after the events, as were cited before for wildfire and harvesting (Dymov, 2023; Peter-Contesse et al., 2024; Medeiros Naval et al., 2025; Chung et al., 2026). These increasing values of soil density can be related to losses of soil carbon (Martínez Pastur et al., 2021, 2022; Lin et al., 2025; Alves et al., 2026). In forest environments, higher soil organic carbon content was associated with lower soil density and pH (Teng et al., 2017; Kupka and Gruba, 2022; Peter-Contesse et al., 2024; Medeiros Naval et al., 2025). These changes can be explained by losses due to the impacts (harvested products or burned soil components), lower leaf litter or wood debris inputs, or increasing soil microbial activity (Chaves et al., 2024; Toledo et al., 2024, 2025). Besides, harvesting machinery traffic disrupts soil aggregates, exacerbates soil compaction, degrades soil structure, and significantly increases soil density (da Luz et al., 2022; Mileusnić et al., 2022; Dymov, 2023; Chung et al., 2026). Soil carbon depletion generates lower soil moisture capacity and increasing soil acidity (Pastor and Post, 1986; Berthrong et al., 2009; Vanderkelen et al., 2026), where this difference in acidity can be generated through several mechanisms, including increased production of organic acids or generation of carbonic acid from higher rates of autotrophic respiration (Richter and Markewitz, 1995). In addition, this carbon depletion affects nutrient cycles and nutrient availability for plants and microfauna, compromising the support for plant root growth (Jobbagy and Jackson, 2003; Tiunov and Scheu, 2004; Anitha et al., 2024; Medeiros Naval et al., 2025; Kaushal et al., 2026).

The magnitude of changes due to human disturbances depends on: (i) type of impacts (harvesting with or without wildfire), (ii) magnitude of impacts, (iii) negative synergies generated among ecosystem components, and (iv) vegetation types. Forests affected by harvesting presented higher soil recovery than stands affected by harvesting followed by wildfires, as well as higher forest structure values and coarse-woody-debris after the events (Alauzis et al., 2004; Fajardo and Gundale, 2015; Rosas et al., 2024; Martínez Pastur et al., 2026). Besides, soil impacts can be related to the magnitude of events, e.g., the quantity of harvested trees or the combined damage due to harvesting followed by wildfire. Many authors have related the magnitude of human impacts over stand characteristics in Southern Patagonia (Klein et al., 2008; Peri et al., 2019; Chaves et al., 2023; Aravena Acuña et al., 2023; Ruggirello et al., 2023, 2025; Martínez Pastur et al., 2024b; Nosetto et al., 2024; Toledo et al., 2024, 2025; Barberá et al., 2025); however, few papers have explored the negative synergies and recovery rates (Rosas et al., 2024; Martínez Pastur et al., 2021, 2026). In our study, when harvesting was affected by wildfires, the magnitude of soil impacts was magnified more than harvesting alone, generating negative synergies that affect recovery on time. Besides, the recovery was different at each forest type, where N. antarctica recovered faster than N. pumilio stands. These findings suggest that N. pumilio has lower resilience to face human related impacts than N. antarctica forests. This is consistent with previous studies that identified N. antarctica as a more plastic and resilient species, e.g., root sprout capacity allowed the species to quickly recover the lost tree crown cover after regenerating the impacted stands (Gyenge et al., 2011; Soto et al., 2019; Soliani et al., 2021; Soler et al., 2022; Martínez Pastur et al., 2025b, 2026). Besides, N. antarctica was identified as a seral species, while N. pumilio was pointed out as the climax species (Gutiérrez, 1994; Frangi et al., 2005a, 2015). In this context, multivariate and index analyses (Figures 3, 4) allowed us to determine different recovery characteristics that differed among forest types: (i) Soil properties after-60-years-impacts in light harvesting (NP-LH and NP-LH+F) and heavy harvesting (NP-HH) are closely related to natural N. pumilio forests (NP), while heavy impacts (NP-HH+F) presented major differences and soil value losses. (ii) In N. antarctica forests, soil properties after-60-years-impacts in harvested stands (CNA-H and ONA-H) are closely related to natural forests (CNA and ONA), while stands affected by harvesting followed by wildfire (ONA-H+F) presented greater losses in their soil values. (iii) Finally, mixed stands presented intermediate values, showing greater resilience than pure stands. Many studies analyzed the natural dynamics of these species (e.g., Gutiérrez, 1994; Frangi et al., 2005a, 2015), and their resilience and recovery values must be integrated into the future management proposals when it is necessary longer periods of time to revert the impacts. Our study suggests that treatments that get closer to grassland characteristics must be considered that have surpassed the ecological limits of sustainability. These findings contribute to developing better soil conservation practices and to improving management proposals (Peri et al., 2017, 2022), giving clues for restoration activities considering soil recovery in impacted stands (Fajardo and Gundale, 2015; Valenzuela et al., 2018; Zamorano-Elgueta and Moreno, 2021; Garcia-Braga et al., 2024).

As presented here, multiple high-severity impacts (e.g., wildfire and harvesting) within short intervals of time can produce novel conditions that exceed the adaptive forest resilience (Enright et al., 2015; Turner et al., 2019), permanently modifying soil physicochemical properties, nutrient concentrations, and nutrient stocks (Johnstone et al., 2016; Dove et al., 2020; Santiesteban-Serrano et al., 2025; Walia et al., 2025). In this context, one major challenge of sustainable forest management and conservation lies in defining the resilience of managed ecosystems (Peri et al., 2017, 2022; Martínez Pastur et al., 2026), e.g., hybrid and novel ecosystems can be created due to the lack of recovery of the managed or impacted stands (Radeloff et al., 2015; Frangi et al., 2015; Ramírez et al., 2025). In this context, adaptive forest management and long-term monitoring must be implemented in the managed landscapes, reducing the frequency and magnitude of impacts to maintain the natural resilience of the natural environments.

5 Conclusion

Natural landscapes generate a gradient of soil characteristics associated with forest structure and vegetation types, from closed Nothofagus pumilio to open N. antarctica forests and grasslands. The human-related impacts (e.g., harvesting with or without wildfire) modified soil physicochemical properties, nutrient concentrations, and nutrient stocks according to their magnitude and vegetation types. Harvesting followed by wildfire generates negative synergies on forest soil when occurring together, being more severe in N. pumilio than in N. antarctica forests (mixed or pure stands). Multivariate and indexes based on forest structure and soil characteristics allow us to characterize the recovery of the different forest types, and the potential generation of hybrid or novel ecosystems linked to soil and forest structure degradation. These findings allow us to identify soil indicators and design more efficient long-term monitoring to implement better adaptive management in temperate forests. In this study we cannot split the effect of the wildfire from previous harvesting impact and represent the main limitation to split the effect of both impact types. Besides, it was not possible to include many other variables that can be considered in future studies (e.g., livestock).

Generative AI statement

The author(s) declared that Generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffgc.2026.1945424/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Yanbo Yang, Zhejiang Agriculture and Forestry University, China

Reviewed by: Ignacio Nuñez-Hidalgo, Pontificia Universidad Católica de Chile, Chile

Gaia Mascetti, University of Milan, Italy