Landscape Patterns Modified by Harvesting and Ranching Allowed the Invasion of Understory Species, Affecting Sustainability of Management Practices in Southern Patagonia

1. Introduction

Forest landscape patterns reflect the combined influence of relatively constant factors (e.g. climate and topography), species distributions (e.g. grasses, herbs, and trees), and its ecological inter-relationships [1,2,3]. Besides, forest landscapes reflect the effects of natural and human disturbances [4,5]. Despite that most of the forest ecosystem are able to recover from major perturbations within decades to half-centuries [6,7], some legacies of forest disturbance can persist and alter the original forest understory [8], modifying forest stand structure [9], and increasing the susceptibility to subsequent disturbances [10]. Further, this susceptibility is closely related to the magnitude of these impacts [11].

Opportunistic invasive plants take advantage of the negative effects of disturbances and the lack of resilience of the impacted forests, which negatively affect the native vegetation communities [12,13]. These invasive plants can cause significant ecological and economic harms in natural and managed forests [14,15], e.g. altering ecological functions, including primary productivity, nutrient cycling, carbon sequestration, and tree regeneration [14,16]. The relationship between invasion success and resource availability is well-documented [13,17]. The likelihood of establishment by invasive species is higher in disturbed areas with more resources (e.g. nutrient-rich sites) and depends on the biotic interactions among native and exotic plants [18,19]. Many patterns and plant strategies of invasion associated have been described including environmental and landscape factors, such as soil, climate, land use or anthropogenic disturbances [13,18,20].

One major challenge of forest management and conservation lies in defining the threshold of resilience of natural ecosystems [11], e.g. where hybrid or de-novo ecosystems can be created due to the lack of recovery capacity of the managed or impacted stands [15]. Human land uses increase landscape fragmentation, and in consequence, the susceptibility to biological invasions [16,21,22]. In order to assess the impact of invasive species and to conserve biodiversity in human-dominated landscapes, we need to consider the role of the landscape context and how it modulates the modifications of natural species assemblages [23].

Over the last century in Patagonia, human activities have become the main driver of change for native forests [11], where the main economic activities (e.g. extensive grazing and harvesting) alter the natural vegetation cover at landscape level [24]. Furthermore, human disturbances (e.g. clear-cuts and fires) significantly change the forest structure, modify soil properties, main ecological processes, and vegetation composition [9,25]. In fact, forest degradation is the main consequence of these extreme impacts that totally change the provision of different ecosystem services, and reduce their resilience capacity [26]. Vegetation changes are frequently associated to the dominance of non-native over native plant species [27], however, the increase of native species associated to open environments (e.g. grasslands) could also be a consequence of changing environmental conditions in natural forested landscapes [9].

Since European colonization (1850-1950), the Patagonian forests in Argentina suffer different human impacts related to the improvement of provisioning ecosystem services which lead to changes in understory species to increase growth and palatability of forage for livestock [7,9,11]. One of the most affected forest landscapes was the Nothofagus antarctica (commonly named ñire) forests in ecotone areas, due to the relevant characteristics for grazing and timber for rural construction purposes [28]. In this context, one of the most extreme impacts was generated through intentional fires to decrease or remove the forest cover in order to maximize forage for livestock. More recently, silvopastoral systems was proposed as more sustainable alternative (e.g. thinning) to obtain poles for fences and lumber, and to open the forest canopy to stimulate the understory growth but maintaining shelter for cattle during winter [9,28]. However, grazing is the impact that prevails in all landscapes (open-lands and forest-lands), generating positive and negative trade-offs with the other described impacts [11].

Data collected through monitoring, especially over long periods, becomes indispensable for evaluating the consequences of ecosystem changes and supporting subsequent decision-making processes [29]. Monitoring provides key insights in ecology, environmental change, natural resource management and biodiversity conservation [30,31], and specifically, on forest ecosystems, the long-term studies have been widely used to monitor changes in forest structure, composition and services [29,32]. Over time, the need for measurable, simple, financially feasible and reliable indicators for biodiversity has increased, but their development has resulted in an inhomogeneous landscape of quantitative and qualitative biodiversity indicators [33,34]. These indicators can be positively or negatively impacted by harvesting or management actions [11,35]. Understanding the link between biodiversity indicators and management measures has been the focus of many comparative, descriptive, and experimental research studies in previous decades, but the relationship to management measures is still elusive [36]. Although the establishment and spread of invasive species is affected by a combination of biotic, abiotic, and landscape factors [16,37], studies that integrate multiple types of environmental data in predicting invasive species distributions are also scarce [16,38,39].

Bio-indication of abiotic site conditions from environmental relationships of plant species has a long tradition [40,41]. Vascular plants were used for environmental indicators due to their ecological behavior relating to main environmental factors modified by natural and human related impacts, including Patagonia in Argentina [24,42]. Data describing environmental indicators can be spatially and temporally referenced to understand changes to the environment over space and time, and are thus an important tool for decision-making [43]. The indicator plant species in silvopastoral systems in Patagonia were chosen according to its correlation to: (i) losses in forage productivity due to the introduction of exotic species (e.g. quantity and quality of forage palatability), and (ii) losses of environmental quality (e.g. soil erosion) that allowed the establishment of undesirable plant species [44]. Exotic species have invaded these austral forests since European colonization, and became one of the main drivers of change, modifying the species assemblage of the native forests. For example, Rumex acetosella has become one of the most frequent introduced species in disturbed environments since 1580 in the Magellan Strait [45,46], while Achillea millefolium was first reported in 1906 for Argentina, and was identified as a weed invading native grasslands [47]. Hieracium pilosella was more recently informed during the 1990s [48], especially after disturbances (e.g. fires and overgrazing) [49,50,51], and was pointed as the biggest threat for profitable livestock farming [52,53]. Further, some studies reported the increase of native plant species associated with open habitats in the impacted forests due to degradation processes [53], e.g. the encroachment of unpalatable plants in forested areas for cushion shrubs as Azorella caespitosa, A. trifurcata and Bolax gummifera, which decline the productivity of grasses associated to over-grazing, clear-cuttings or intense human-related fires [42,54].

Most of these studies try to understand the invasive ecology of the species [15,48,49,50,51,53], while other research analyzed vegetation cover changes related to specific impacts [24,32,42,44,47,52,55]. However, to understand how specific species change in a management forest landscape, it is necessary to analyze plant species considering all the environmental and the associated impacts. The objectives of this study is to define the landscape and environmental characteristics that allowed the invasion of understory species, and define the thresholds of the human related impacts to propose better sustainable management practices. We want to answer the following questions: (i) do the environmental characteristics of the well-conserved natural ecosystems (topography, soil, forest structure) reduce the vascular plant invasion?, (ii) do the changes of human related impacts (livestock, harvesting, fires) on the environmental characteristics facilitate the vascular plant invasion?, (iii) are there detectable thresholds for the changes generated by management practices (harvesting, livestock)?, and (iv) Are there indicator species for the different kinds of impacts (livestock, harvesting, fires) and environments (open-lands and forest-lands)?

2. Materials and Methods

2.1. Study Area

Sampling was conducted at El Roble ranch, located in the central area of Tierra del Fuego (Argentina) at 54º 01′ to 54º 02′ SL and 67º41′ to 67º43′ WL (Figure 1A,B). We selected a forested landscape with a long history of impacts and human uses, which covers 432 ha of ñire forests and associated open-lands (dry and wet grasslands). This forested landscape was mainly used for cattle grazing during the last century, and different management practices were implemented within this area: (i) thinning to promote silvopastoral systems (SPS), increasing pastures under tree cover, and obtain firewood and poles for fence construction, (ii) clear-cuts and pastures for livestock, and (iii) fires to remove trees and woody debris, and accommodate pastures for livestock [9,11,55]. A total of 165 plots were measured, using a systematic grid of 150 x 150 m across forest-lands and open-lands (Figure 1C).

2.2. Sampling Design, Data Taking and Calculations

Each plot was sampled using a 25 m transect from east to west of the center geographic point defined by the grid. We characterized each plot by their forest structure, environmental variables, animal uses, and understory plants. Forest structure was measured by one plot in the center of the transect using the point sampling method (BAF = 4) [56], obtaining the basal area (BA, m² ha^-1) with a Criterion RD-1000 (Laser Technology, USA) and the dominant height (DH, m) as the two tallest trees per plot using a TruPulse 200 (Laser Technology, USA). BA was also relativized considering the maximum value of the total sampling (BAmax, %) to define one tree occupancy index (values from 0% to 100%). These data were complemented with hemispherical photographs, which were taken at the center of transects at 1.5 m above ground level with an 8-mm fisheye lens (Sigma, Japan) mounted on a 35 mm full-frame digital camera (Nikon, Japan) with a tripod leveling head to ensure horizontal lens position. Each photograph was orientated with the upper edge towards the magnetic north, avoiding direct sunshine under evenly overcast skies or cloudless days. Gap Light Analyzer software v.2.0 [57] was used to define cover (overstorey crown and debris) (CC, %), as a percentage of open sky, and total direct radiation (DR, W m²) transmitted through the canopy and debris. The user-supplied input variables details were presented in Martínez Pastur et al. [58]. Slope (SLO, %) was also measured using an inclinometer in the center of transects. To characterize the soil layer, we collect two samples along each transect using a field borer (0-30 cm) of known volume after previously removing the litter layer. Samples were weighed before and after air-drying in laboratory conditions (24°C) until constant weight, obtaining soil moisture (SM, %) and soil bulk density (SBD, gr cm³) after coarse root debris and stones >2 mm were removed by sieving. For chemical analyses, we pooled the individual soil samples into one combined sample per transect. Soil samples were finely ground to below 2 mm using a tungsten-carbide mill, and determined: (i) soil acidity (pH) in a suspension (air-dried samples and deionized water) of a soil:water ratio of 1:2.5 [59], (ii) soil organic matter (SOM, %) and soil organic carbon (SOC, %) based on a dry combustion analysis (muffle furnace) and modeling (see [63]), (iii) soil nitrogen (SN, %) was determined with a semi-micro Kjeldahl method [60], and (iv) soil extractable phosphorus (SP, ppm) was determined by Bray and Kurtz [61] method.

We also recorded animal feces from native (Lama guanicoe, guanaco) and domestic animals (sheep, cows and horses) along each transect (25 × 4 m, 100 m²). We used the number of feces as a proxy of animal use (animals ha-1) by assuming: (i) Feces was maintained in the forest floor for a period not longer than a calendar year (decomposition rates was low during summer due to low temperatures and high during winter due to mechanical effect of snow accumulation) [62]; (ii) Average values of defecation per day were defined as 6.0 times per day for guanacos, and 12.3 times per day for cows and horses; (iii) The requirements of dry matter forage (palatable plants) varied according to the animals (650 kg DM yr^-1 for guanacos and 3250 kg DM yr^-1 for cows and horses); and finally, (iv) we used sheep equivalent (SE, animal ha^-1) based on the animal species size to standardize the current livestock density (for more details see [63]).

To characterize the understory, we used the point-intercept method [64] with 50 intercept points (every 0.5 m) along each transect. At each point, we recorded bare ground, debris and vascular and bryophyte species and calculated cover based on number of hits on each species along the transect: e.g. bare ground cover (BG, %), debris cover (>2.5 cm diameter) (DC, %), bryophyte cover (mosses and liverworts) (BC, %), and vascular plant cover including: tree regeneration less than 1.30 m height (RC, %), monocot plant cover of native species (MONO, %), monocot plant cover of exotic species (MONO-E, %), dicot plant cover of native species (DICO, %), and dicot plant cover of exotic species (DICO-E, %). We also survey the occurrence of degradation as indicator plant species at each sample point of each transect, including: Azorella caespitosa (AZCA), A. trifurcata (AZTR), Bolax gummifera (BOGU), which are mostly associated to open-lands, and Achillea millefolium (ACMI), Hieracium pilosella (HIPI), and Rumex acetosella (RUAC) which were more related to forest-lands. The first ones have been identified as early-seral species that settled after disturbances, and the second ones have been identified as species that reduce the quality of the understory quality for grazing purposes [9,53,63]. These data allowed us to determine plant species cover or species group cover (dicots, monocots, native and exotic) following Moore [65] and Correa [66].

For further analyses, the sampled plots were integrated into a geographical information system (GIS), and Sentinel 2 images (10 m pixel) were used to classify in: (i) environment types defined as open-lands (OL) and forest-lands (FL). Open-lands including dry grasslands (G) dominated by Festuca gracillima and Empetrum rubrum and wet grasslands (W) dominated by Juncus scheuchzerioides, Carex curta, C. macrosolen and Caltha sagittata [9,53]. Forest-lands were classified as closed forests (CF, >60% CC), open forests (OF, <60% CC), and edge forests (EF, 25 m strip of ecotone forests between open-lands and forest-lands) (Figure 2 and Figure 3A). (ii) Grazing impacts were defined as low (LG) and high grazing pressure (HG). Grazing impact (values from 0% to 100%) was calculated as the total animal use (native and domestic) relativized with the maximum value of sheep equivalents (SE, n ha^-1) recorded during sampling. Low grazing was considered when values were <20%, while high grazing was considered when values were >20% (Figure 3B). And finally, we considered (iii) the harvesting and fire impacts on the forest-lands, which were classified as unharvested areas (UH, natural forests without management or fire evidences), low harvested stands (LH, evidence of previous harvesting and CC >60%), high harvested stands (HH, evidence of previous harvesting and CC <60%), and stands with fire impacts (F), despite the intensity and kind of damage, detected during field sampling (Figure 2 and Figure 3C).

2.3. Statistical Analyses

One-way ANOVAs were used to compare environment types (OL-G, OL-W, FL-CL, FL-EF, FL-OF), grazing (OL-LG, OL-HG, FL-LG, FL-HG), and harvesting or fire impacts (FL-UH, FL-LH, FL-HH, FL-F), using Fisher test and Tukey test at p <0.05 and Statgraphics software (Statistical Graphics Corp., USA). These treatments and levels were used to analyze: (i) forest structure (CC, DR, DH, BAmax), (ii) topography and soil properties (SLO, SBD, SM, pH, SOC, SOM, SN, SP), (iii) understory plant cover (BG, DC, RC, BC, MONO, MONO-E, DICO, DICO-E), and (iv) indicator species cover (AZCA, BOGU, AZTR, HIPI, ACMI, RUAC).

Principal Component Analysis (PCA) was performed to evaluate patterns of plot distribution of environment types, impacts and indicator species. In the first analyses, we evaluate the plots by environment types (OL and FL) and impacts (grazing, harvesting, fires), as response to the complete pool of variables (four of forest structure, eight of topography and soil properties, and eight of understory cover = 20 variables). In the second group of analyses, PCA were employed to explore patterns in the cover of six degradation indicator plant species with the environment types and impacts. We selected correlation coefficients among columns to obtain the final cross-product matrices. Finally, we also performed the Multi-Response Permutation Procedure (MRPP) to analyze statistical differences among groups of plots for each comparison. These analyses were conducted in PC-ORD 5.0 [67].

ANOVA and PCA and MRPP were performed using the complete pool of plots and variables (n = 165) when analyzing environment types and grazing impacts, but only used FL plots (n = 111) when harvesting and fire impacts were considered. For the PCA of degradation indicator plant species, only plots with cover of 10% or higher for at least one species were employed (n = 66 for environment types and impacts, and = 49 for harvesting and fire impacts). Finally, we analyzed the relationship (average ± standard error) between crown cover (CC) and animal density (SE, livestock and guanacos) considering: (i) environment types (OL-G, OL-W, FL-CF, FL-EF, FL-OF), (ii) grazing impacts (OL-LG, OL-HG, FL-LG, FL-HG), (iii) harvesting and fire impacts in forest-lands (FL-UH, FL-LH, FL-HH, FL-F), and (iv) degradation indicator plant species cover (AZCA, BOGU, AZTR, HIPI, ACMI, RUAC) compared with the average cover for all the measured plots (TOTAL).

3. Results

3.1. Landscape Characterization

Sample points were classified according to the different environment types and human related impacts (Figure 3A), where 54 plots were classified as OL (33%) and 111 as FL (67%), and occurring in an intermingled landscape with grasslands and not in continuous patches. Open-lands included 29 plots of dry (54%, brown) and 25 plots of wet (46%) grasslands. Forest-lands included 51 plots of closed forests (46%), 31 plots of edge forests (28%), and 29 plots of open forests (26%). Grazing impacts varied greatly across the landscape, where the higher pressure was related to the open-lands and the forest stands located close to these environments (Figure 3B). Our sampling detected more plots with high (78%,) than low grazed grasslands (22%), and more plots with low (55%) than high grazed forests (45%). Finally, considering harvesting and fire impacts in FL (Figure 3C), we found more unharvested stands (51%) located far-away to the roads, than low (14%) or heavy (8%) harvested stands. Fires (27%) included more plots than harvesting, and occurred during different periods and events despite the distance to roads.

3.2. Changes of the Environment Characteristics across the Landscape

Forest structure, topography and soil properties showed significant differences among the different environment and impact types (Table 1 and Table 2). Forest structure showed differences among OL and FL levels, where forest stands presented higher crown cover and lower total direct radiation values (Table 1A). Besides, closed forests showed the highest values of forest structure (CC = 71.4%, DH = 8.5 m, BAmax = 51.1%), while edge and open forests did not present significant differences between them. Most of the topography and soil properties only showed differences between wet grasslands and other treatments (dry grasslands and forests) (Table 2A). Wet grasslands presented the lowest values of slope (SLO = 2.1°) and soil bulk density (SBD = 0.37 gr cm³), with the greater values of soil moisture (SM = 107.3%), soil organic carbon (SOC = 22.3%), soil organic matter (SOM = 55.5%) and soil nitrogen (SN = 1.3%). However, forest-lands presented the highest values of pH and soil phosphorus (closed > edge > open forests).

When we analyzed forest structure values (Table 1B), grazing impacts showed significant differences between land types (FL and OL) but not between intensity types. However, we found differences for both dimensions when we analyzed the topography and soil properties (Table 2B). As was expected, slope was higher on FL than OL, but grazing intensity increased at higher slopes. Besides, the livestock generated significant impacts over soil bulk density, and significantly increased the values on both OL and FL. Livestock also influenced the other studied soil properties, reducing the values according to impact intensity (low > high) and environment type (FL and OL). Some values were higher at OL (SM, SOC, SOM, SN), while pH and SP depended on the combination of both variables (e.g. these variables were higher at FL-LG and minimum at OL-HG).

Harvesting and fire impacts significantly influenced forest structure variables, and were slightly influenced by topography and soil properties (Table 1C and 2C). As was expected, high harvesting (CC = 35.5%, HD = 7.4 m, BAmax = 13.3%) and fires (CC = 38.1%, HD = 7.6 m, BAmax = 20.4%) showed the lowest values of forest structure variables, while low harvesting presented the higher values (CC = 77.3%, HD = 8.8 m, BAmax = 60.8%) due to the secondary forests presented abundant advanced regeneration. Unharvesting stands presented medium values, due to mature and old-growth forests presented open canopies and trees growing at different development and decay stages (Table 1C). Only one soil property showed significant differences, where fires increased (pH = 5.2) and harvesting decreased (pH = 4.7-4.9) the values compared to unharvested stands (pH = 5.0) (Table 2C).

3.3. Vegetation Cover

Understory vegetation cover showed significant differences among different environment and impact types (Table 3). Three variables presented the highest values in dry grassland areas (OL-G), including bare ground (BG = 16.1%), bryophytes cover (BR = 2.3%) and native dicots cover (DICO = 16.8%). Native monocot cover showed the highest values in wet grasslands (MONO = 67.0%) and the lowest values in closed forests (MONO = 25.1%). The exotic species cover (MONO-E and DICO-E) increased with tree crown cover (closed > edge > open forests), as well as at different land type (forests > open-lands), where closed forests presented the highest values of exotic monocot (MONO-E = 26.5%) and exotic dicot cover (DICO-E = 18.7%). Dry grasslands showed the lowest values of cover of exotic species, which indicates the greater ecosystem integrity of their original species assemblage. Finally, as was expected, closed forests presented the highest cover of debris (DC = 4.7%) and open forests presented the highest cover of regeneration (RC = 5.5%) (Table 3A).

Grazing intensity increased bare ground, both at open-lands and forest-lands. Also, debris cover was higher in the high grazing areas, probably related to higher impacts of grazing in the forest (e.g. harvesting generates higher debris and open canopies, which promotes higher grass development). In open-lands, grazing also occurred in areas with greater bryophyte cover (e.g. more humid areas) and higher native dicot cover. Besides, in open-lands the higher exotic species covers were found in areas with low grazing pressure, maybe because they are preferred for livestock during foraging. In the forest-lands higher grazing pressure was found in areas with higher exotic monocot cover (pastures) but lower native dicot cover (Table 3B). To fully understand these trends, it is necessary to understand the changes produced by the human related impacts (harvesting and fires) in the forest-lands. Fires promote higher native monocot cover, but with lower exotic monocot cover, while light harvesting generates better conditions for exotic species cover (MONO-E and DICO-E) and lower native monocot cover. Finally, high harvesting areas promote more bryophyte cover (BC = 1.8%) than other treatments (Table 3C).

3.4. Degradation Indicator Plant Species

Only four degradation indicator plant species showed significant differences considering different environment and impact types (Table 4). The two species Achillea millefolium (ACMI) and Rumex acetosella (RUAC) did not differ across types of environment, grazing harvest or fire regimes, but in general trends were higher at forest-lands. As was expected, Azorella caespitosa (AZCA) and Bolax gummifera (BOGU) showed the highest values of cover in open-lands, however AZCA was higher at wet grasslands (3.0%), while BOGU was higher at dry grasslands (3.4%). AZCA greatly decreased with grazing intensity (6.1% to 1.6%) in open-lands, while BOGU insignificantly increased with grazing at both open-lands and forest-lands (1.7% to 2.3% and 0.06% to 0.12%, respectively). Azorella trifurcata (AZTR) naturally grows in open-lands, but our samplings showed higher abundance in open forests (FL-OF and FL-EF > FL-CF and open-lands). This species decreased with grazing at both land types (6.7% to 3.9% in open-lands and 14.0% to 9.9% in forest-lands). Besides, AZTR significantly increased at high impacted forests (e.g. FL-F reaching to 22.5% cover). Finally, Hieracium pilosella (HIPI) increased from wet to dry grasslands (1.1% to 2.6%), and when tree canopy cover decreased (from 1% in CF to 7.1% in OF). However, the different impact types not significantly influenced over its cover. Trends showed marginal increases due to grazing (1.0% to 2.1% in open-lands and 2.0% to 4.1% in forest-lands), where the maximum covers were found at higher harvesting stands (7.9%).

3.5. Relationship among Environmental Variables and Impact Types

PCA showed different ordination patterns depending on the variables used to characterize the sampling plots by environment or impact types (Figure 4). Environment types were clearly split along the first axis explaining 31.5% of the variation, where the second axis separated within forested or open land types explaining 19.0% of the variation. Forest-lands (CF, EF, OF) were slightly overlapped with less dispersion than open-lands (Figure 4A). However, open forests and wet grasslands presented some overlap, generating a continuous gradient among the environment types. These dispersion trends can be associated with CC gradient. As was expected, closed forests (green dots) were positively related to forest structure variables (CC, BAmax, DH), but also with two soil properties (pH, SP) and three understory cover variables (DC, DICO-E, MONO-E). Besides, two soil properties (SLO, SBD) and two understory cover variables (RC, DICO) were associated to edge (grey dots) and open-forests (pale blue dots), and share some variables (TR, BG, BC) with wet grasslands (orange triangles). Four soil properties (SM, SOC, SOM, SN) and one understory variable (MONO) were associated to dry grasslands (orange triangles). Pairwise comparisons among treatments through MRPP showed significant differences among all groups (p <0.01). Grazing impacts also separated the environment types (OL in purple triangles and FL in blue dots) (Figure 4B), but the level of grazing impacts was strongly overlapped among treatments. Some variables can be associated to the grazing impacts: (i) forest structure (CC, BAmax), soil (SP, pH) and understory cover (DICO-E, MONO-E) can be related with low grazing in FL (light blue dots), while most of the understory variables (DC, RC, DICO, BG, BC), two soil properties (SLO, SBD) and one forest structure (DH) can be associated to areas with high grazing in FL (blue dots). Finally, BC, TR and MONO can be related to areas with high grazing in OL (purple triangles), while only soil properties (SOC, COM, SM, SN) can be related to areas with low grazing in OL (light purple triangles). Pairwise comparisons among treatments through MRPP showed significant differences for most of the groups (p <0.01), except between grazing levels in OL (p = 0.09) and grazing levels in FL (p = 0.07). In these two analyses (Fig 4A and 4B), axis 1 (Eigenvalue = 6.305, explained variance = 31.5%) was mainly influenced by CC, BAmax, TR, DH, SOC, SOM, SN, and MONO (Eigenvectors of 0.33, 0.30, -0.32, 0.31, -0.29, -0.29, -0.27, and -0.31, respectively). Axis 2 (Eigenvalue = 3.976, explained variance = 19.9%) was mainly related with SBD, SP, and BG (Eigenvector of -0.31, 0.32, and -0.32, respectively). The Eigenvalues in the PCA for axes 1 and 2 were significant (both presented p <0.001).

PCA that analyzed forest-land types clearly split the closed forests (green dots) with edge forests (grey dots) and open forests (pale blue dots) that presented a high overlapping (Figure 4C). Closed forests can be related to forest structure (CC, DH, BAmax) and understory cover variables (DC, DICO-E, MONO-E), while the edge and open forests con be related to the greater total direct radiation at ground level (TR) and the total monocot plant cover (MONO). Pairwise comparisons among treatments through MRPP showed significant differences among all groups (p <0.01). PCA that analyzed harvesting and fire impacts showed some trends, but with high overlapping among treatments (Figure 4D). Unharvested forests (green dots) and fires (red dots) presented a great dispersion at both axes, evidencing a great amplitude of environmental conditions that were naturally occurred across the landscape. Besides, some soil properties (SOC, SOM, SN, SBD) and understory variables (BC, BG, DICO) can be related with high harvested stands (purple dots), while some forest structure (CC, BAmax, DH) and understory variables (DC, MONO-E, DICO-E) can be associated to low harvested stands (light purple dots). Pairwise comparisons among treatments through MRPP showed significant differences for most of the groups (p <0.01), except for fires and heavy harvesting stands (p = 0.18). In these two analyses (Fig 4C and 4D), axis 1 (Eigenvalue = 4.915, explained variance = 24.6%) was mainly influenced by CC, BAmax, TR, MONO, MONO-E and DICO-E (Eigenvector of -0.41, -0.38, 0.39, 0.33, -0.31, and 0.30, respectively). Axis 3 (Eigenvalue = 2.068, explained variance = 10.3%) was mainly related with SLO, pH, SN, RC and DICO (Eigenvector of -0.40, -0.43, 0.31, -0.29, and 0.29, respectively). The Eigenvalues in the PCA for axes 1 and 3 were significant (both presented p <0.01).

PCA showed different ordination patterns, depending of the variables used to characterize the sampling plots by degradation indicator plant species cover (Figure 5). Environment types were clearly split between OL and FL, but with high overlapping among types (Figure 5A). Two indicator plant species (BOGU and AZCA) can be related to dry and wet grasslands, while the other species can be related to a wide range of forested environments. Pairwise comparisons among treatments through MRPP showed significant differences for most of the groups (p <0.01), except between grassland types (p = 0.44), and edge and open forests (p = 0.64). Grazing impacts also split the environment types (OL in purple dots and FL in blue dots) (Figure 5B), but once again, the level of grazing impacts was strongly overlapped among treatments. No clear associations were found between indicator plant species and grazing levels. Pairwise comparisons among treatments through MRPP showed significant differences for most of the groups (p <0.01), except between grazing levels in OL (p = 0.34) and grazing levels in FL (p = 0.09). In these two analyses (Fig 5A and 5B), axis 1 (Eigenvalue = 1.568, explained variance = 26.1%) was mainly influenced by AZTR, BOGU and AZCA (Eigenvectors of 0.70, -0.49, and -0.42, respectively). Axis 2 (Eigenvalue = 1.320, explained variance = 22.0%) was mainly related with RUAC and AZCA (Eigenvector of 0.67, and -0.44, respectively). The Eigenvalues in the PCA for axes 1 and 2 were non-significant, however axis 2 presented marginally significant differences (p = 0.12 and p = 0.06, respectively).

PCA that analyzed forest-land types related to degradation indicator plant species cover showed that closed forests (green dots) presented greater homogeneity compared to edge forests (grey dots) and open forests (pale blue dots) that presented a high dispersion (Figure 5C). Some species showed relationships with the closed forests (ACMI, HIPI), but the other species were related to particular canopy conditions of the open forests. Pairwise comparisons among treatments through MRPP showed non-significant differences (p >0.05), except between landscape types (OL and OF, p <0.01). PCA that analyzed harvesting and fire impacts showed some trends, but with high dispersion among treatments (Figure 5D), e.g. burnt stands showed greater relation with AZTR. In these two analyses (Fig 5C and 5D), axis 1 (Eigenvalue = 1.732, explained variance = 28.9%) was mainly influenced by AZTR, RUAC and BOGU (Eigenvector of 0.61, -0.56, and -0.49, respectively). Axis 2 (Eigenvalue = 1.245, explained variance = 20.8%) was mainly related with AZCA and BOGU (Eigenvector of 0.79, and 0.43, respectively). Pairwise comparisons among treatments through MRPP showed that half of them showed significant differences (p <0.01), and half showed non-significant differences, including closed forests compared to heavy (p = 0.92) and light harvested stands (p = 0.09), and between light and heavy harvested stands (p = 0.21). The Eigenvalues in the PCA for axes 1 and 2 were non-significant, however axis 1 presented marginally significant differences (p = 0.083 and p = 0.360, respectively).

The relationships between crown cover (CC) and animal density (SE, livestock and guanacos) synthetize the univariate and multivariate analyses (Figure 6), highlighting the coincidences between the two major studied impacts (CC as proxy of harvesting, and SE as proxy of animal uses). The environment types comparisons showed a gradient of CC, where FL-OF and FL-EF are closer to open-lands than FL-CF, and where SE not greatly changed, except for FL-EF which presented higher average values (Figure 6A). Interestingly, the lower and higher grazing intensity occurred at similar crown covers (Figure 6B). When we analyzed the impacts of harvesting and fires, some unexpected results were obtained. It was expected to be a clear trend between CC and SE, when CC increases, SE decreases, e.g. these trends can be observed between harvested and unharvested forests (FL-HH > FL-UH > FL-LH from higher values of SE to higher values of CC). Fires presented low values of CC, but also presented lower values of SE. Fires and FL-HH presented similar CC values, however, the harvesting areas presented higher SE levels. Probably, these differences showed some degradation processes (Figure 6C). Finally, the degradation indicator plant species showed a clear relationship between CC and SE (from lower to higher values), where BOGU < AZCA < HIPI-AZTR < RUAC < ACMI (Figure 5D).

4. Discussion

4.1. Changes in Forest Structure, Soil Properties and Understory Cover in the Landscape

A mosaic of vegetation types characterizes the ecotone of Southern Patagonia, where ñire forests often intermingle with diverse landscapes, including grasslands. Those landscapes are characterized by several environmental factors and the availability of different natural resources [9]. Our sampled plots were distributed throughout different environments that represent this landscape. Our approach considered both natural conditions, such as transition zones and soil moisture levels, and human impacts influences like extensive grazing, which often lead to additional activities like harvesting and fire aimed at enhancing forage availability in forested areas [28].

The differences in forest structure attributes were directly related to transition zones. Open and edge forests exhibited low canopy cover and high total radiation at the understory level, while close forests displayed high dominant height, aligning with observations in other ñire forests [53]. This direct relationship showed the importance of transition zones, where changes influencing the availability of critical resources (e.g. light and moisture) for understory plant growth. Although we anticipated these differences in forest structure might induce changes in soil properties, we found that the most significant fluctuations occurred primarily in wet grasslands characterized by high soil moisture conditions, organic carbon and organic matter, distinguishing them from other areas. While most soil properties remained relatively consistent across different forest structures [58,68], our study revealed that pH and soil phosphorus values increased with crown cover. Besides, the availability of limitant resources significantly affected understory species cover. The impact of these factors was significant, where soil properties and forest structure played a crucial role in shaping the plant community. These findings align with previous research that underscores the influence of environmental factors and resource availability on the establishment and persistence of plant species [69,70]. In our study, we found that grasslands were better preserved than forests, showing the highest cover values for native species and bryophytes. Besides, consistent with the findings of other studies, our results also revealed a correlation between changes in forest structures, soil phosphorus availability and the increase of exotic species cover [16].

4.2. Environmental Changes and Indicator Plants Related to Human Impacts

In ecosystems under the pressure of human impacts, new environmental conditions and resources pose novel challenges for both resident and invading species. The establishment and expansion of plant species depend on how closely these new conditions align with historical baselines [71]. Several studies have demonstrated how harvesting and fire impact over light availability, and how grazing affects soil resources, which in turn can influence the cover of understory plant species [9,28]. Our study revealed that high grazing pressure significantly affected soil properties, leading to increase soil bulk density and reduce soil carbon, organic matter, nitrogen and pH across the different land types [72]. Consequently, grazing had a substantial influence on species cover, resulting in decreased native species cover in open lands and forests, while the response of exotic species cover varying based on different combinations of factors. On the other hand, harvesting and fire significantly increased solar radiation, while fire directly affected soil conditions, raising pH levels [73]. These impacts facilitated the establishment of exotic and native understory species that increased under similar conditions to those found in more open forested areas [9,53,73].

Understanding how human impact modifies natural structures that facilitate the invasion of understory species is critical for assessing the magnitude of these impacts, especially in ecosystems where productivity and ecological resilience are highly affected [55]. In fact, early detection of indicator species is crucial for mitigating the long-term impacts of new species introductions or expansion of already established species [16,39]. In ñire forest-lands, several studies have linked the presence of specific native and exotic plant species to human impacts, making them excellent indicators of environmental and forage productivity degradation [44,53]. Two native species of our study, Bolax gummifera and Azorella caespitosa, were more prevalent in open-lands, particularly in dry grasslands. Grazing intensity significantly influenced the presence of these species, where A. caespitosa decreased with higher grazing intensity, while B. gummifera increased where soil became more acidic. Additionally, A. trifurcata, primarily associated with open-lands, showed increased cover in open forests and high-impact forested areas [42,55]. The increase of B. gummifera and A. trifurcata are associated with soil erosion conditions [44]. In fact, A. trifurcata displayed grazing response characteristics, such as rosette growth forms, which are more tolerant to trampling compared to erect forms, leading to increased competition for resources and space [74].

Apart from changes in local factors, landscape-level factors such as fragmentation also played a role in facilitating the invasion and establishment of invasive species [23]. In the context of ñire forests, human land use particularly practices like grazing and harvesting, increases landscape fragmentation, making these ecosystems more susceptible to invasions [21,22,39]. Our results showed an increase of cover of exotic species in forests subjected to high pressure of grazing and harvesting. Most of those species were intentionally introduced to improve native grasslands (e.g. Rumex acetosella) over a century ago [75], however their cover increased after disturbances, particularly fire and overgrazing in forest areas [76]. Moreover, Hieracium pilosella varied with wetness and tree canopy cover, where cover increased as areas became drier (e.g. dry grasslands) and with reduced tree canopy cover (e.g. open forests). Although this species is considered the most significant threat to profitable livestock farming [52,53], our study did not find significant differences in species cover in response to the impacts analyzed across different environments.

4.3. Relation between Invasive Understory Species and Forest Harvesting

Invasive species that take advantage of disturbances or reduced resilience in impacted forest pose significant ecological and economic challenges. Their presence has adverse effects on critical ecosystem functions including primary productivity, nutrient cycling and tree regeneration [14,16,39]. To effectively address these challenges and mitigate the associated ecological and economic consequences, it is imperative to gain a deep understanding of the resilience thresholds within natural ecosystems and how these thresholds are altered by human activities. This knowledge forms the foundation for developing more sustainable forest management and conservation strategies. Although invasive plant issues are often experienced by many landowners and managers at the regional level, the practical management of invasive plants typically takes place at the stand scale [18]. In this context, the development of early detection tools becomes crucial, as they can assist owners and forest managers in identifying alien plants in their early stage of establishment [77]. This approach is particularly important for species that are recognized as invasive in other countries, such as H. pilosella in Europe. Additionally, adaptive management practices in the framework of silvopastoral systems, can be effectively employed to regulate inter-species competition. These systems help to reduce the establishment of invasive species while promoting the growth of native species [28,32]. The removal of disturbances (e.g. reducing animal density) is known to benefit the native species and results in a gradual reduction in invasive species succession over time [13]. However, it is essential to recognize that the trajectory of plant community succession can follow multiple paths when disturbances and invasive species play significant roles [78]. In this regard, long-term monitoring of plant invasions in forests becomes a necessary tool for unraveling the patterns of invasion longevity (stage) and the factors influencing the direction of these changes [13].

5. Conclusions

We can conclude that ñire forest landscapes involve the interaction of various factors, including environmental characteristics, species distributions, and their ecological interactions, alongside the effects of both natural and human disturbances. We found considerable variation in forest structure, soil properties and vegetation cover across different environments at forest landscape level. Human impacts exhibited profound effects on different variables and land types, emphasizing the importance of managing intensities. The use of indicator plant species allows to identify potential ecological threshold related to human impacts and the establishment of species linked to ecological and economic degradation, e.g. Bolax gummifera and Azorella trifurcata were associated with high grazing pressure in grasslands and fire in forested areas, while Rumex acetosella and Achillea millefolium, typically associated with forested areas, were related to high harvesting pressure and fire impacts. Although the influence of these impacts is well-documented, our study provides a deep understanding of the environmental patterns at forest landscape level and how human impacts can alter these factors affecting understory cover and facilitating the establishment of invasive species.