Compared to Wildfire, Management Practices Reduced Old-Growth Forest Diversity and Functionality in Primary Boreal Landscapes of Eastern Canada

Large primary forest residuals can still be found in boreal landscapes. Their areas are however shrinking rapidly due to anthropogenic activities, in particular industrial-scale forestry. The impacts of logging activities on primary boreal forests may also strongly differ from those of wildfires, the dominant stand-replacing natural disturbance in these forests. Since industrial-scale forestry is driven by economic motives, there is a risk that stands of higher economic value will be primarily harvested, thus threatening habitats, and functions related to these forests. Hence, the objective of this study was to identify the main attributes differentiating burned and logged stands prior to disturbance in boreal forests. The study territory lies in the coniferous and closed-canopy boreal forest in Québec, Canada, where industrial-scale logging and wildfire are the two main stand-replacing disturbances. Based on Québec government inventories of primary forests, we identified 427 transects containing about 5.5 circular field plots/transect that were burned or logged shortly after being surveyed, between 1985 and 2016. Comparative analysis of the main structural and environmental attributes of these transects highlighted the strong divergence in the impact of fire and harvesting on primary boreal forests. Overall, logging activities mainly harvested forests with the highest economic value, while most burned stands were low to moderately productive or recently disturbed. These results raise concerns about the resistance and resilience of remnant primary forests within managed areas, particularly in a context of disturbance amplification due to climate change. Moreover, the majority of the stands studied were old-growth forests, characterized by a high ecological value but also highly threatened by anthropogenic disturbances. A loss in the diversity and functionality of primary forests, and particularly the old-growth forests, therefore adds to the current issues related to these ecosystems. Since 2013, the study area is under ecosystem-based management, which implies that there have been marked changes in forestry practices. Complementary research will be necessary to assess the capacity of ecosystem-based management to address the challenges identified in our study.

Large primary forest residuals can still be found in boreal landscapes. Their areas are however shrinking rapidly due to anthropogenic activities, in particular industrial-scale forestry. The impacts of logging activities on primary boreal forests may also strongly differ from those of wildfires, the dominant stand-replacing natural disturbance in these forests. Since industrial-scale forestry is driven by economic motives, there is a risk that stands of higher economic value will be primarily harvested, thus threatening habitats, and functions related to these forests. Hence, the objective of this study was to identify the main attributes differentiating burned and logged stands prior to disturbance in boreal forests. The study territory lies in the coniferous and closed-canopy boreal forest in Québec, Canada, where industrial-scale logging and wildfire are the two main stand-replacing disturbances. Based on Québec government inventories of primary forests, we identified 427 transects containing about 5.5 circular field plots/transect that were burned or logged shortly after being surveyed, between 1985 and 2016. Comparative analysis of the main structural and environmental attributes of these transects highlighted the strong divergence in the impact of fire and harvesting on primary boreal forests. Overall, logging activities mainly harvested forests with the highest economic value, while most burned stands were low to moderately productive or recently disturbed. These results raise concerns about the resistance and resilience of remnant primary forests within managed areas, particularly in a context of disturbance amplification due to climate change. Moreover, the majority of the stands studied were old-growth forests, characterized by a high ecological value but also highly threatened by anthropogenic disturbances. A loss in the diversity and functionality of primary forests, and particularly the old-growth forests, therefore adds to the current issues

INTRODUCTION
Human-induced climate change, large-scale biodiversity loss and the global disappearance of intact ecosystems are major challenges faced by humanity in the 21st century (Steffen et al., 2015;IPCC, 2018;Watson et al., 2018a). Forest ecosystems play and will continue to play a key role in the management of these issues. They are indeed important habitats for biodiversity, and their functions, such as carbon sequestration or water provision, provide many ecosystem services that have a crucial value for both indigenous as well as industrial societies (Gauthier et al., 2015a;Watson et al., 2018b). However, forests have been under increasing pressure from human activities over the centuries. Only about 23.5% of current forest ecosystems are still intact (Potapov et al., 2008) and can therefore be considered as "primary, " i.e., large, continuous forested landscapes undisturbed by human activities (Kormos et al., 2018). Emphasis has been placed over recent decades on management approaches closer to natural dynamics to ensure the sustainable management of forest ecosystems (Kuuluvainen, 2002;Gauthier et al., 2009;Puettmann et al., 2009). In this context, primary forests are models to be studied thoroughly to develop relevant management models capable of addressing contemporary challenges.
Remnant primary forests situated in the boreal biome, mainly in Canada and Russia (Mackey et al., 2014;Potapov et al., 2017), together constitute a significant part of the Earth's primary forested land. Their area is however decreasing rapidly due to anthropogenic activities (Morales-Hidalgo et al., 2015;Potapov et al., 2017;Grondin et al., 2018). Old-growth forests (i.e., forests relatively old and driven by secondary disturbances; Oliver and Larson, 1996;Wirth et al., 2009) are often abundant in primary boreal areas (Shorohova et al., 2011). It is important to stress here that "primary forest" and "old-growth forest" do not represent the same ecological concepts, although these terms are often used interchangeably (Wirth et al., 2009). A primary forest may contain stands of all ages (e.g., forests that burned recently or a long-time ago in a landscape never disturbed by anthropogenic activities), while an old-growth forest may be a previously managed forest that has been abandoned for several centuries. Previous research highlighted that industrialscale forest management in primary forests led to a degradation of old-growth forest surfaces (Östlund et al., 1997;Cyr et al., 2009;Grondin et al., 2018), changes in tree species composition (Bouchard and Pothier, 2011;Boucher and Grondin, 2012;Kuuluvainen et al., 2017), landscape homogenization and fragmentation (Schmiegelow and Mönkkönen, 2002;Haeussler and Kneeshaw, 2003) and a decrease in deadwood richness Moussaoui et al., 2016). Moreover, recent studies have emphasized the wide diversity of forest types (i.e., stands defined by specific structural attributes and/or tree species composition) in boreal landscapes, in particular once the old-growth stage is reached (Shorohova et al., 2009;Martin et al., 2018;Moussaoui et al., 2019). These results underscore that boreal forests are complex webs of changing structures driven by secondary disturbances (i.e., disturbances that are not stand replacing, such as low severity insect outbreaks or partial windthrows) variable in nature, severity and temporality (Shorohova et al., 2009;Martin et al., 2018Martin et al., , 2019).
Yet, little is known about how forest management may impact the diversity of forest types in primary landscapes compared to natural disturbances. Previous research has mainly focused on changes in landscape age structure and/or tree composition (Östlund et al., 1997;Cyr et al., 2009;Boucher et al., 2017a). Others observed a posteriori if protected and/or unmanaged forests were representative of the surrounding areas (Lõhmus et al., 2004;Joppa and Pfaff, 2009;Svensson et al., 2020). Stands of a similar age-class and tree composition can, however, be defined by different structural characteristics and may provide different habitats for biodiversity (Rheault et al., 2009;Martin et al., 2018). As forest structure and composition changes over time, the different forest types that can be observed in boreal landscapes are interrelated (Martin et al., 2018). The disappearance of specific forest types may for this reason subsequently impact others, and consequently the associated habitats. For example, the peak in wood volume generally comes before the peak in deadwood richness and diversity (Brassard and Chen, 2006). The harvesting of forests at optimal timber volumes can therefore lead to a regression of deadwood habitats in the long term (Stokland et al., 2012). Forests left undisturbed by human activities (protected or not) are also generally defined by a lower economic value (Joppa and Pfaff, 2009;Belote and Aplet, 2014;Sabatini et al., 2020). This implies that in addition to landscape rejuvenation, homogenization and fragmentation, forest management may also lead to the rarefaction of specific forest types and habitats in logged areas, probably those with the highest economic value.
The coniferous boreal landscapes of Québec, Canada are relevant areas to evaluate the impacts of anthropogenic activities on primary forests. Indeed, industrial-scale forest management started relatively recently (around the 1960's; Boucher et al., 2017b) and many primary forests remain. Fire is the main type of natural stand-replacing disturbance (i.e., disturbance of high severity that reinitiates forest succession) in this region and is still an important driver of the landscape dynamics Gauthier et al., 2015b). Insect outbreaks and windthrows are the other main agents of natural secondary disturbances in this territory but are rarely stand-replacing (Shorohova et al., 2011;De Grandpré et al., 2018). It has been demonstrated that industrial-scale forest management is not necessarily an analog to fire in terms of landscape and stand scale impacts (Haeussler and Kneeshaw, 2003;Boucher et al., 2017b;Kuuluvainen and Gauthier, 2018). Much of the controversial impact of forest management in boreal territories is due to the overuse of largescale clearcuts with a rotation inferior to the onset of the oldgrowth stage (Haeussler and Kneeshaw, 2003;Kuuluvainen and Gauthier, 2018;Martin et al., 2020a). To address this issue, the Province of Québec has been following the principles of ecosystem-based forest management since adopting them in 2013 (Gouvernement du Québec, 2018); the aim is to emulate natural disturbance regimes . Within this framework, primary forest landscapes become references for evaluating management effectiveness. The impact of logging activities in this area should be as close as possible to that of natural disturbances, in particular fire. The comparison of the characteristics of logged and burned stands over recent decades would therefore help to evaluate how close forest management currently comes to following natural disturbance regimes. In addition, this will contribute to enriching our understanding of the ecological quality (e.g., in terms of time since the last fire, structural attributes and tree species composition) of remnant primary forests in managed landscapes.
Since the 1960's, Québec's government has performed several field inventories on a large spatial scale, including in the primary boreal landscapes [Ministère des Forêts de la Faune et des Parcs (MFFP), 2016b]. Many of these field plots have therefore been logged and burned since being surveyed. In this way, plotlevel data collection makes it possible to compare the impact of management practices and forest fires on Québec's primary boreal forests by examining the characteristics of stand samples that predate the disturbance. The objective of this study is thus to identify the main structural and compositional attributes that differentiate burned and logged stands in coniferous boreal forests of Eastern Canada. We hypothesized that in comparison to fire, logging affects stands with a higher economic value (e.g., high merchantable wood volume or tree density), thereby threatening specific primary forest types (i.e., stands defined by specific structural attributes and tree species composition). This study will help clarify and shed new light on current issues arising from the gap between traditional forest management and a natural disturbance regime. The results obtained could serve as a baseline for defining sustainable forestry strategies in the new context of ecosystem-based management in Eastern Canada. On a larger scale, this study will also yield a more precise portrait of how anthropogenic activities have altered forest landscapes in regions where primary forests are now scarce.

Study Area
The study area covers an area of 412,400 km 2 , representing the limit of the bioclimatic domain of the coniferous and closed-canopy boreal forest (black spruce (Picea mariana Mill.)feather moss bioclimatic domain) in Québec, Canada (Figure 1). Its latitude ranges between the 49 • N and 52 • N parallels, while its longitude ranges from the western to the eastern limit of the Province of Québec. Black spruce forests dominate this territory, accompanied by balsam fir (Abies balsamea (L.) Mill.), jack pine (Pinus banksiana Lamb.), paper birch (Betula papyrifera Marsh.), and trembling aspen (Populus tremuloides Michx.) (Ministère des Forêts de la Faune et des Parcs (MFFP), 2018). Mean temperature ranges from −5 to 1 • C and annual precipitation (rain and snow) from 700 to 1,200 mm/year (Saucier et al., 2009). The territory is divided in two bioclimatic subdomains which each have specific features with regard to climate, physiography, natural disturbances and forest composition (Figure 1, Saucier et al., 2009). A drier, more continental climate and shorter fire cycles (often <200 years) characterize the western subdomain. Black spruce forests with jack pine are well-represented. In contrast, a wetter, more maritime climate and longer fire cycles (generally >200 years) are observed in the eastern subdomain. Mixed black spruce-balsam fir forests are more abundant. Fires are mainly stand-replacing in the study area (Shorohova et al., 2011) and their susceptibility increases along a south-north latitudinal gradient, while forest productivity decreases (Gauthier et al., 2015b). For this reason, a northern limit for commercial forests has been assessed in 2005 (re-evaluated in 2016 with almost no changes) and logging is prohibited beyond this limit (Jobidon et al., 2015). Logging activities started in the middle of the 20th century at the southern edge of the boreal forests and has progressed northward since (Boucher et al., 2017b). Prior to this period, most of the study area was primary forest and around half of its surface could still be considered primary forest in 2016 (Potapov et al., 2017). The vast majority of the forests that were logged in the studied territory during recent decades were therefore harvested for the first time. Due to the low wood volume and the remoteness of the forests, logging is essentially done with heavy machinery. Short-rotation (70-100 years) clearcuts are the preferred silvicultural treatment in these forests. If natural regeneration is lacking after the harvest, plantations can be conducted to ensure forest renewal. Protection of natural regeneration and soil improved progressively over time, for example by limiting machine movement within the harvested areas.

Québec's Decadal Surveys of Temporary Plots
The Ministère des Forêts, de la Faune et des Parcs du Québec (Québec's Ministry of Forests, Wildlife and Parks; MFFP) has performed four decadal forest surveys since the 1960s, and the fifth is still ongoing (final data will be available in 2027). These surveys include the sampling of 400 m 2 circular temporary plots (i.e., plots that are sampled during only one survey) all over the province, in order to obtain a representative description of Québec's forests. Temporary plots are defined following a linear transect, with a distance between each plot generally ranging between 100 and 400 m. Transect locations were defined using a stratified sampling design that aimed to characterize the FIGURE 1 | Map of the study territory. The insert map indicates its location in Canada. different stand types (age class and species composition) present in this region. The results of these surveys are thus expected to be representative of the studied landscapes [Ministère des Forêts de la Faune et des Parcs (MFFP), 2016b]. Diameter at breast height (DBH) was measured for each merchantable tree (DBH ≥ 9 cm) in the 400 m 2 plot or in a 40 m 2 subplot for saplings (DBH > 1 cm, DBH < 9 cm and height ≥ 1.3 m), and their species was identified. Sampled trees and saplings were divided into 2 cm DBH classes and according to tree species. The following structural data were then computed for each species/DBH class by the MFFP: density (n/ha), basal area (m 2 /ha), and wood volume (m 3 /ha). In addition, one to three study trees representative of the plot were selected in each temporary plot among the merchantable trees. The height of these study trees was measured, as was their minimum age using a core sample taken at 1 m from the ground. For our analysis, we computed MFFP temporary plot data based on nine attributes ( Table 1). We also used topographic data from the United States Geological Survey and National Aeronautics and Space Administration (USGS/NASA) Shuttle Radar Topography Mission (Farr et al., 2007) to compute the altitude and slope of each temporary plot using ArcGIS [version 10.7; Environmental Systems Research Institute (ESRI), Redlands, California].
For this study, we first retained only the results of temporary plot surveys that were publicly available for the study territory, i.e., those sampled during the second (1980-1991), third (1992-2004), and fourth (2005-2016) decadal forest surveys. Second, we retained temporary plots defined as mature forests [i.e., canopy height > 7 m; Ministère des Forêts de la Faune et des Parcs (MFFP), 2016a]. This selection was based on the fact that some data, such as those related to study trees, are not collected by the MFFP if these conditions are not met. Moreover, all temporary plots are sampled in forests considered as economically viable, i.e., stands that will reach a wood volume of at least 30 m 3 /ha 120 years after the last fire [Ministère des Forêts de la Faune et des Parcs (MFFP), 2016b]. Unlike logging, fire burns both productive and non-timber productive forests (Madoui et al., 2010(Madoui et al., , 2015. By analyzing only productive stands, we avoid this bias due to the differences in "behavior" between fire and logging. We then used public maps of forest fires and logging provided by the province of Québec to identify the temporary plots that had been disturbed by stand-replacing fires or logging after sampling. We only selected transects for which all the temporary plots were subsequently disturbed by the same type of standreplacing disturbance (i.e., fire or logging) during the same time period. Finally, we only selected transects that were below the northern limit for commercial forests. Indeed, logging activities We chose to perform the analysis at the transect scale because some temporary plot data (e.g., stand height or tree age) rely on a low tree sample. Hence, we used an approach similar to that of a previous study (Bouchard et al., 2008) that also grouped plot values at the transect scale using similar data in the coniferous boreal forest of Québec. This approach has the advantage of reducing variation associated with local forest conditions and emphasizing large-scale processes, e.g., at the scale of the bioclimatic domain. Transect attributes were thus equal to the mean of the corresponding temporary plot attributes, except for estimated stand age, for which we retained only the highest value. This choice is justified by the low number of cored trees per temporary plot, as well as the coring methodology (i.e., core sampled at a height of 1 m), which may underestimate tree age in the case of juvenile tree suppression (McCarthy and Weetman, 2006;Marchand and DesRochers, 2016;Martin et al., 2019). In this way, we considered the age of the oldest tree in the transect as the best available indicator of the time since the last stand replacing disturbance at the transect scale. In the study area, 423 transects met our criteria (total number of temporary plots: 2,352), with 201 burned and 222 logged. These areas were burned by at least 79 different fires that occurred between 1986 and 2014, or disturbed by logging activities that took place between 1985 and 2016. The transects studied were thus representative of the stand-replacing disturbance regime in the study territory over the last three decades. The mean number of temporary plots per transect was equal to 5.7 ± 1.5 plots for burned transects and 5.4 ± 1.6 plots for logged ones. The mean time elapsed between mature forest sampling and subsequent disturbance was equal to 11.4 ± 6.6 years for burned transects and 8.8 ± 6.7 years for logged ones. Compared to the temporal scale of forest dynamics, we considered this a relatively short lapse of time (Bergeron and Harper, 2009;Martin et al., 2018;Portier et al., 2018a). Consequently, we assumed that stand structure and composition at the time of sampling were similar to those at the time of disturbance.

Most Discriminant Transect Attributes Between Burned and Logged Forests
To test our hypothesis that, compared to fire, forest types with higher economic value (e.g., high merchantable wood volume or tree density) are affected more by logging, we first aimed to identify the most discriminating transect attributes according to the type of disturbance (i.e., fire or logging) using permutational analyses of variance (PERMANOVA) (Anderson, 2001). Each transect attribute was analyzed separately and without transformation. Analyses were performed with 10,000 permutations using Euclidean distances. As mixed black sprucebalsam fir forests are more frequent in the eastern black sprucefeather moss bioclimatic subdomain and black spruce forest is more abundant in the western black spruce subdomain, the potential influence of bioclimatic subdomain should be tempered. Hence, nested PERMANOVAs were performed for each transect attribute, with the disturbance type (burned or logged) as fixed effect and the bioclimatic subdomain as random effect (western or eastern subdomain). We then used the random forest classification algorithm (Breiman, 2001) as a complementary analysis to determine which transect attributes better discriminated burned and logged transects. The mean decrease accuracy (MDA) and the mean decrease Gini (MDG) values obtained from a random forest tree analysis are efficient indicators to identify the most relevant variables for random forest classification. High MDA and MDG values for a given variable indicate a high ranking for an efficient classification (Han et al., 2017). The dependent variable was the disturbance type and the independent variables were the transect attributes, as well as the bioclimatic subdomain. Preliminary analyses indicated that discriminative ability (constant out-of-bag error value) of the random forest algorithm reached its peak at 5,000 trees created and two nodes per tree subdivision.

Identification of the Most Logged and Burned Forest Types
To identify if predominantly burned or logged forests corresponded to distinct forest types defined by homogeneous structural and compositional attributes, we used the Ward's linkage method (Ward, 1963), with preliminarily standardized data and Euclidean distances. The transect attributes retained for clustering of forest types were those that showed significant differences according to disturbance type on PERMANOVA analysis and that showed high MDA and MDG values on the random forest test. We aimed to achieve a balance between a number of forest type clusters high enough to be representative of landscape structural diversity, yet small enough to permit ecological interpretation. To better highlight the distribution of clusters according to the disturbance type, we also carried out a principal component analysis (PCA), based on the standardized transect attributes used for the clustering of forest types. We then compared transect attributes among the forest types using PERMANOVA, with 10,000 permutations and Euclidean distances. When the results of the test were significant, we performed a pairwise PERMANOVA (Martinez Arbizu, 2017), with 10,000 permutations, Euclidean distances and Bonferroni adjustment for pairwise comparisons. The proportion of burned and logged transects between forest types was also compared using pairwise chi-square tests (Mangiafico, 2016) with the Bonferroni adjustment for pairwise comparisons.

Differences in the Structural and Compositional Attributes of the Logged and Burned Stands
Wood volume, basal area, tree density, stand height, and jack pine proportion in logged transects significantly differed from those that had burned (p < 0.001), as did sapling density, broadleaved species proportion and slope (p < 0.05) (Figure 2). Wood volume, basal area, tree density, stand height and broadleaved proportion were, respectively, 45.6, 33, 20, 19.5, and 41.5% higher in logged stands compared to the burned ones. Logged stands were thus denser and more productive than burned stands, resulting in a higher wood volume. In contrast, sapling density, jack pine proportion and slope were, respectively, 10.5, 42.1, and 13.4% lower in logged stands compared to those burnt. Disturbance type did not significantly influence estimated stand age, black spruce proportion, balsam fir proportion, or altitude. We also observed higher proportion of broadleaved species in logged forests in comparison to the burnt.
The out-of-bag error value of the random forest classification was equal to 28.8%, implying that around three-quarters of the burned or logged transects presented highly distinctive characteristics. Stand height was the transect attribute presenting higher MDA and MDG values, followed by wood volume, basal area, tree density and black spruce proportion, respectively (Figure 3). Balsam fir proportion, jack pine proportion, broadleaved proportion, estimated stand age, sapling density and bioclimatic subdomain were the transect attributes presenting the lowest MDA and MDG, but with differences in ranking order between the two metrics. Stand height was thus the main transect attribute discriminating logged from burned stands, followed by wood volume.

Clustering of Logged and Burned Stand Structures and Composition
Stand height, wood volume, basal area, and tree density were the transect attributes used to define the forest types, as they showed significant results on the nested PERMANOVA as well as was the highest MDA and MDG values. We determined that five forest type clusters were the optimal balance between statistical significance and ecological interpretability. The four transect attributes considered for the clustering were sufficient to explain 94.6% of the variance between sites on the first two axes of the PCA (Figure 4). The number of transects was well-balanced among the forest types, ranging from 56 to 130 transects per forest type (Forest type 1: 56 transects, Forest type 2: 130 transects, Forest type 3: 82 transects, Forest type 4: 95 transects, Forest type 5: 60 transects, mean = 84.6 transects/forest type). Almost all the transect attributes presented significant differences within the forest types at p < 0.001, except for estimated sapling density, estimated stand age and altitude (p = 0.003, p = 0.002 and p = 0.009, respectively). This implies that each forest type represented a specific forest structure and composition, as well as specific topographic attributes ( Figure 5).
The forest types were distributed along an increasing gradient from Forest type 1 to Forest type 4 or Forest type 5 for wood volume (39.1-155.7 m 3 /ha), stand height (11.6-18.2 m), basal area (8.8-26.9 m 2 /ha), tree density (642.1-1,419 trees/ha), balsam fir proportion (4.9-19%), broadleaved species proportion (4.3-17.5%), and slope (2.1-4.1%). However, wood volume, basal area and tree density were lower for Forest type 5 than Forest type 4 (respectively, 141.5 m 3 /ha, 21.8 m 2 /ha and 854.1 trees/ha for Forest type 5; 155.7 m 3 /ha, 26.9 m 2 /ha and 1,419 trees/ha for Forest type 4), but its stand height was superior (18.2 m for Forest type 5 and 16.1 m for Forest type 4). For the other transect attributes, differences between the forest types were less striking, with the exception of jack pine proportion, which was significantly higher for Forest type 1 (35%) than for the other forest types (8.5-15.7% for Forest types 5 and 3, respectively). Sapling density showed an opposite gradient from that previously observed, with the highest values for Forest type 1 (3390.8 saplings/ha) and the lowest values for Forest type 5 (2370.8 saplings/ha). Forest type 2 was significantly older than Forest type 1 (161.9 and 130.7 years, respectively), with no clear difference found for the other forest types. Similarly, Forest type 3 was defined by a higher proportion of black spruce (80.7%) than Forest types 1, 4 and 5 (68.1, 67.2, and 65.2%, respectively). Finally, the altitude of Forest types 2 and 5 (450.8 and 457.6 m, respectively) was superior to that of Forest types 1 and 4 (414.3 and 405.7 m, respectively), although their values were relatively close. Overall, the five forest types represented specific structural and compositional attributes: open black spruce-jack pine forests, of low productivity (Forest type 1), open black spruce forests (Forest type 2), dense black spruce forests (Forest type 3), dense mixed forests (Forest type 4), and highly productive mixed forests (Forest type 5).
The frequency of burned and logged transects differed markedly between the forest types. We observed an inverse gradient starting from Forest type 1 (mainly burned, rarely logged) to Forest type 5 (rarely burned, mainly logged) (Figure 6). The relative frequency of burned transects was equal to 91.1, 59.2, 46.3, 25.3, and 18.3% from Forest type 1 to Forest type 5, respectively. Conversely, the relative frequency of logged transects was equal to 8.9, 40.8, 53.7, 74.7, and 81.7% from Forest type 1 to Forest type 5, respectively. An even proportion of burned and logged transects was therefore reached only for Forest type 3, and, to a lesser extent, Forest type 2 and 4. Overall, burned and logged forests were defined prior to disturbance by significantly different structures and compositions in Eastern Canadian boreal forests. Some forest types were particularly targeted by logging activities while others seem to have been noticeably avoided. Red dots indicate transect attributes that presented significant differences based on the disturbance type on the nested PERMANOVA test, white dots non-significant results and gray dots indicate transect attributes that were not analyzed with the nested PERMANOVA test. Prop., proportion.
FIGURE 4 | Biplot of the PCA, presenting the transects studied according to the transect attributes considered for clustering, the forest types and the type of disturbance, on the two first axes of the PCA (PC1 and PC2, respectively). Numbers in brackets indicate the percentage of variance explained by each axis of the PCA, red crosses indicate the location of the transect attributes.

DISCUSSION
The impact of forest harvesting on the study area between 1985 and 2016 differed significantly from that of fire. Logging activities mainly targeted the most productive primary forests, defined at the time of disturbance by higher stand height, wood volume, basal area and tree density than the burned stands. Some specific forest types were logged to much higher extent than they were burned, while others were consistently not selected for logging, hence supporting our hypothesis. As a consequence, logging was not simply another disturbance in addition to fire in managed boreal territories. Forestry practices acted as a divergent disturbance that may alter boreal forest type diversity and lead to habitat depletion as well as a loss in landscape functionality. Ecosystem-based management, that began to be applied 3 years before the end of the study period, should therefore aim to reduce the gap observed between former forestry practices and natural disturbance dynamics. This would ensure the sustainability of anthropogenic activities in boreal forests. It is in particular necessary to balance harvests along the productivity gradient highlighted in this study, with a lower logging pressure on the most productive stands. Continuous cover forestry should also be promoted as much as possible in place of clearcuts to maintain or restore old-growth forest attributes.

Divergent Impacts of Logging Compared to Fire May Threaten Boreal Old-Growth Forest Functionality and Diversity
Burned and logged transects were distributed along an inverse forest type gradient: low productive forest types were almost exclusively burned (Forest type 1, open black spruce-jack pine forests, low productivity) while highly productive forest types were almost exclusively logged (Forest types 4 and 5, mixed forests, dense or highly productive). Differences in disturbance frequencies (burned vs. logged) were less striking for the intermediary forest types (Forest types 2 and 3, black spruce forests, open or dense) (Figures 5, 6). These results are consistent with previous studies in boreal forests, showing that remnant natural forests are often found in less productive environments (Schröter et al., 2014;Martin et al., 2020a;Sabatini et al., 2020). A steeper topography also often explains the absence of anthropogenic disturbances (Joppa and Pfaff, 2009;Belote and Aplet, 2014;Svensson et al., 2020). We found in the study area the opposite trend, probably because the topography of the boreal landscapes of eastern Canada is relatively smooth and a slightly steeper slope improve stand productivity due to a better drainage (Martin et al., 2018;Laamrani and Valeria, 2020).
The majority of the studied transects were also old, with an estimated stand age exceeding 100 years for 80.3% of the studied transects (79.1 and 81.5% of the burned and logged transects, respectively) ( Figure 5). This indicates that most of the stands studied exceeded the commonly used Eastern Canadian old-growth age threshold for the boreal forest, i.e., 100 years (Bergeron and Harper, 2009), independently of the disturbance type. The methodology used to estimate minimum tree age gives however little detail on the age structure of the forests studied (e.g., even-aged or uneven-aged, presence of survivors of old disturbances). An uneven-aged, multicentury-old forest, for example, may indeed be characterized by relatively young dominant trees (Martin et al., 2020b,c), whereas some young, even-aged forests may have survivors from the previous fire (Smirnova et al., 2008). This predominance of old trees is nevertheless consistent with Eastern Canadian boreal forest ecology, in which old-growth forests are dominant because of relatively long (>200 years) fire cycles (Cyr et al., 2009;Gauthier et al., 2015b;Grondin et al., 2018). Trees that survived the last fire are probably rare, as fires are generally stand-replacing in the study area (Shorohova et al., 2011). Consequently, our study provides precious insights about the impacts of forest management over the last 30 years not only in primary forests, but also on boreal old-growth forests in Eastern Canada. Specifically, this research highlights that the structural and compositional attributes of remnant old-growth forests in managed areas may differ strongly from those observed in natural landscapes. It is also likely that their ecological quality, e.g., in terms of large tree density or deadwood volume, is significantly lower in managed areas.
The concentration of logging activities in the denser, most productive old-growth stands raises important concerns about the sustainability of past management practices in the study territory. This pattern indeed strongly contrasts with that of fire. More generally, these results raise questions about how remnant primary forest are represented in the study area. Boreal old-growth forest types range from dense and closed canopy to more open structures (Shorohova et al., 2009;Martin et al., 2018;Portier et al., 2018a). Our results support these findings, highlighting the variety of structure and composition of the forests that generally exceeded 100 years of age in the study area. This diversity in structure and composition fosters the emergence of multiple habitats, and as a result, forest biodiversity differs significantly from one old-growth structure to another (Desponts et al., 2004;Rheault et al., 2009;Fenton and Bergeron, 2011;Janssen et al., 2011). The most productive forests can be defined by a higher variety of habitats than the less productive ones (Gillman and Wright, 2006;Niedziałlkowska et al., 2010;Martin et al., 2018). For example, some forest species depend on large pieces of dead wood or a complex vertical structure of the canopy (Rheault et al., 2009;Stokland et al., 2012;Schowalter, 2017). The degradation of the most productive old-growth types (i.e., Forest types 4 and 5, mixed forests, dense or of high productivity) in harvested areas may therefore threaten the associated biodiversity due to the disappearance of specific habitats.

Remnant Forests May Be Defined by Lower Resistance and Resilience in Managed Areas
The differences in the characteristics of stands prior to disturbance raise questions about the resistance and resilience (respectively, the ability of a system to absorb small perturbations and prevent them from amplifying into large disturbances, and the ability to recover following disturbance; Perry and Amaranthus, 1997) of remnant forests to future disturbances. Old-growth stands defined by an open canopy, meaning a lower wood volume, are common in boreal landscapes because of the action of secondary disturbances such as windthrows or insect outbreaks of low to moderate severity (Kuuluvainen et al., 2014;Khakimulina et al., 2016;Martin et al., 2020c). Similar characteristics of structure, composition and topography characterized Forest type 2 in comparison to Forest types 3 and 4, but we also observed a lower tree density, basal area and volume. It is thus possible that Forest type 2 grouped together recently disturbed forests, explaining this more open structure despite the similarities with more productive forest types. A relatively open structure caused by secondary disturbances is often temporary, as shade-tolerant regeneration quickly reach the canopy and fill the gaps (Pham et al., 2004;Aakala et al., 2007;Martin et al., 2020d). Yet, compound disturbances (i.e., a succession of multiple disturbances of low and moderate severity in a short period of time) may overwhelm stand resistance and resilience. This dynamic has been empirically demonstrated in research on partial or selective cuts in the boreal forest: an excessive cutting rate can lead to significant post-harvest mortality because of subsequent windthrows or tree stress, jeopardizing the future of the stand (Bose et al., 2014;Montoro Girona et al., 2019). For this reason, old-growth stands that have recently undergone secondary disturbances can be considered as more vulnerable to compound disturbances than those that have not been significantly disturbed for a long time. In our study, Forest type 2, which seems to represent relatively productive and recently disturbed forests, was the second most avoided forest type by logging activities. Remnant forests within managed areas may thus contain a higher proportion of recently disturbed stands in comparison to the pre-harvest landscape.
Natural disturbances also leave important biological legacies (e.g., deadwood, seed bank, undisturbed patches) that help the ecosystems to recover (Franklin et al., 2000;De Grandpré et al., 2018). The most productive forests in particular provide a higher deadwood volume and diversity when they burn compared to the less productive (Moussaoui et al., 2016). For this reason, biological legacies present in Forest types 1 and 2 after fire may not be surrogates for those present in more productive forests, due to differences in abundance and quality. Logged forests are also generally characterized by a lower deadwood volume and the scarcity of large pieces of deadwood in comparison to naturally disturbed forests (McRae et al., 2001). Merchantable timber is indeed removed from the stands and some important structural elements, such as large trees, take a long time to reappear after the harvest. It is therefore unlikely that the logging of the most productive forest types provides satisfactory biological legacies that can compensate their harvest.
The characteristics of Forest type 1 and to a lower extent Forest type 2, i.e., low stand height, tree density, basal area or wood volume combined with a low slope, can indicate stands undergoing paludification (Simard et al., 2007;Fenton and Bergeron, 2011;Martin et al., 2018), or stands that experienced recurrent forest fires of moderate to high severity (Smirnova et al., 2008;Côté et al., 2013). Such stands are more likely to present regeneration failures, i.e., a regeneration too limited in density and/or growth to settle a closed-canopy forest. In paludified stands, the organic layer may be too thick and too wet to burn completely. This can result in a layer of residual organic matter that limits seedling establishment and growth after the disturbance (Greene et al., 2007;Terrier et al., 2014;Jayen et al., 2017). Recurrent fire can also reduce seed availability, resulting in poor post-fire regeneration (Smirnova et al., 2008;Pinno et al., 2013). The resulting lack of combustibles may then prevent subsequent fires from reaching the severity necessary to allow the remaining seeds to germinate (Miquelajauregui et al., 2016). The discrepancy between fire and logging dynamics observed in this study may therefore cause marked stand depletion and regeneration failures in landscapes with many logged areas. More generally, the lower resistance and resilience of remnant forests raises questions about their capacity to maintain the functions of primary boreal landscapes.

Differences in Logging and Fire Drivers Explain Their Divergent Impacts on Boreal Landscapes
The differences between logged and burned stands observed in this study probably result from the divergent dynamics and drivers of logging activities and forest fires. Logging in boreal forests is an industrial process, hence profitability and technical feasibility are its primary drivers (Perry, 1998;Puettmann et al., 2009;Halme et al., 2013). For these reasons, the superior wood volume observed in logged transects compared to the burned ones was expected, as harvesting less productive stands could be an unprofitable endeavor. In contrast, fire is a more stochastic disturbance, driven by a complex combination of climatic factors as well as biotic and abiotic landscape characteristics (Bélisle et al., 2016;Erni et al., 2018;Portier et al., 2018b). Some forest stands may thus be burned repeatedly in a relatively short time period while others may remain unburned for centuries, even millennia (Bergeron et al., 2002;Couillard et al., 2016;Kneeshaw et al., 2018). Based on analysis and radiocarbon dating of charcoals from forest mineral soils, Couillard et al. (2016) also underscored that flat to undulated reliefs of lowlands are defined by a shorter fire cycle than hilltops in Eastern Canadian boreal landscapes. According to a more recent study (Martin et al., 2018), conducted in the same territory as Couillard et al. (2016), the transect attributes of Forest type 1 were close to that of even-aged and old-growth stands situated in riverside lowlands, defined by sandy deposits and where jack pine (a species adapted to recurrent fire) is abundant. Forest type 5 structural characteristics were close to those of stands situated on hilly reliefs, where balsam fir (a species not adapted to fire) is wellestablished. Poorly drained lowlands were rarely considered in our study, as they are generally not economically viable. This can explain the predominance of jack pine and fires in the lowlands studied, as the sampled sites were mainly situated on dry areas. In contrast, soils have lower drainage on hills, and the slightly higher elevation as well as the denser canopy of these forests can also cause later snowmelt, providing better protection from fire. Forest type 1 therefore represented the sites most prone to fire while Forest type 5 contained the less prone, with the other forest types representing different points of the gradient.
The characteristics of the forest types identified also broadly correspond to the different potential vegetation composition (i.e., the theoretical tree-species composition at the end of forest succession) that can be observed in the coniferous and closed-canopy boreal forest of Eastern Canada along a toposequence (Blouin and Berger, 2004;Saucier et al., 2009;Grondin et al., 2013). Hence, the five forest types identified in this study capture the diversity of forest structure, composition and habitat characteristics distributed from valley to hilltops that characterize the study area (Figure 7). This corresponds globally to the succession of the three main potential vegetation compositions: black spruce-feather moss, black spruce-balsam fir and balsam fir-white birch, respectively. These results imply that the productive and accessible part of the boreal forests of eastern Canada may now be driven by two stand-replacing disturbance types (fire or logging), that are also distributed inversely along the toposequence.

Implications for Ecosystem-Based Management in a Context of Global Change
The findings of our study underscore that the conservation of boreal old-growth forests is not only a matter of quantity, but also of quality. They also mainly illustrate the impact of forest management in Québec's boreal forests prior to the implementation of ecosystem-based management Gouvernement du Québec, 2018). Indeed, while ecosystem-based forest management principles have been in effect in Québec since 2013, the last logging activities considered in this study took place in 2016; the shift in forest management strategy in Québec is thus too recent to be visible in our data. Given this new management paradigm, the limitations of past forest management, demonstrated here, are particularly important to consider. Extensive criticism has been leveled at short logging rotation cycles, compared to longer fire cycles, as well as the predominance of clearcut in boreal forest treatments. This is because of the negative impacts of clearcut logging on forest diversity and functionality (Bergeron et al., 2002;Haeussler and Kneeshaw, 2003;Boucher et al., 2017b).
To address these issues, management plans in Québec now aim to maintain at least 30% of the historical mean of old-growth forest surface areas in the territories where they are applied, with the potential vegetation composition as an indicator of forest diversity [Ministère des Forêts de la Faune et des Parcs (MFFP), 2016a]. Low-productivity and/or more vulnerable oldgrowth forest types (i.e., Forest types 1 and 2) can be present in all the main potential vegetation composition of the coniferous boreal forest (Figure 7), which calls into question the efficacy of current management strategies. It is possible that logging activities will continue to harvest the most productive old-growth stands in the different potential vegetation composition and thus will perpetuate the loss of specific habitats, for example those related to deadwood or complex vertical structures (Lassauce et al., 2011;Stokland et al., 2012;Schowalter, 2017).
Recent studies have emphasized the ecological importance of irregular shelterwood treatments in sufficiently productive boreal stands (Fenton et al., 2013;Franklin et al., 2019;Kim et al., 2021). While more productive stands should nevertheless be better represented in protected and/or unmanaged areas (Joppa and Pfaff, 2009;Price et al., 2020;Sabatini et al., 2020), promoting continuous cover forestry in these stands is also likely to help maintain necessary habitat for many species.
Fire cycles are expected to shorten in the circumboreal zone over the coming decades Boulanger et al., 2013;Oris et al., 2014). However, the impact of climate change on old-growth forest surfaces in boreal landscapes will likely be relatively limited in comparison to that of logging in Eastern Canada (Fall et al., 2004;Bergeron et al., 2010Bergeron et al., , 2017. In this study, we highlighted that the forest types generally avoided by logging activities are also the most susceptible to being burned. If the trend observed in this study continues, it is likely that the remnant primary forests in managed landscapes will also be most vulnerable to the fires that are expected to increase due to climate change. We may see an additive impact of climate change and forest management, reinforcing the current threats on boreal forests. Given this scenario, forest management strategies should consider these risks to better adapt to this emerging context. This is, for example, an opportunity to apply restoration practices on less productive forests (e.g., drainage, plantations) to improve their resistance and resilience. Promoting continuous cover silviculture as an alternative to clearcutting in the most productive forests would also help maintain the ecosystem functions and services they already provide.

CONCLUSION
This study is the first to give a general overview of the impact of more than 30 years of industrial-scale forest management on primary, generally old-growth, boreal forests at the scale of the bioclimatic domain in Eastern Canada. By directly observing the effect of logging on landscapes initially undisturbed by human activities, these results can help to estimate how forest management has influenced the characteristics of remnant primary forests, including outside the study area. The divergent impacts of logging in comparison to fire in Eastern Canadian boreal forests may significantly threaten ecosystem functionality and the associated biodiversity. The higher rate at which the most productive old-growth stands are logged may have resulted in remnant landscapes characterized by a higher proportion of recently disturbed or less productive forest types. This implies that remnant landscapes may be defined by lower resistance or resilience, higher fire susceptibility and greater risk of regeneration failure. Further, logging activities may reinforce landscape vulnerability in a context where the natural disturbance regime is amplified due to climate change. To address this and other issues, ecosystem-based management began to be applied in the study area at the end of the period considered here. This change in practice occurred too recently, however, for its effectiveness to be assessed in this study. Our results can nevertheless serve as a baseline to verify whether the effect of logging approaches that of fire in this new management paradigm. Future research will be thus necessary to evaluate if ecosystem-based management effectively protects remnant primary forests in the boreal landscapes of eastern Canada. Specifically, the diversity of boreal old-growth forests in terms of forest types should be directly addressed in management planning to guarantee that remnant forest structural and compositional diversity remains similar to that of the preindustrial landscape in logged areas.

DATA AVAILABILITY STATEMENT
Publicly available datasets were analyzed in this study. This data can be found at: https://figshare.com/articles/dataset/Dataset_ Martin_et_al_2021_Frontiers_in_Forests_and_Global_Change/ 13626737.

AUTHOR CONTRIBUTIONS
MM and HM conceived the ideas and designed methodology. MM collected, prepared, analyzed the data, and led the writing of the manuscript. M-CL and PG supervised the analysis. MM, PG, YB, HM, and M-CL interpreted the results. All authors contributed critically to the drafts and gave final approval for publication.