Using Naturalness for Assessing the Impact of Forestry and Protection on the Quality of Ecosystems in Life Cycle Assessment

: A novel approach is proposed to evaluate the impact of forestry on ecosystem quality in life cycle assessment (LCA) combining a naturalness assessment model with a species richness relationship. The approach is applied to a case study evaluating different forest management strategies involving concomitantly silvicultural scenarios (plantation only, careful logging only or the current mix of both) combined with an increasing share of protected area for wood production in a Qu é bec black spruce forest. The naturalness index is useful to compare forest management scenarios and can help evaluate conservation needs considering the type of management foreseen for wood production. The results indicate that it is preferable to intensify forest management over a small proportion of the forest territory while ensuring strict protection over the remaining portion, compared to extensive forest management over most of the forested area. To explore naturalness introduction in LCA, a provisory curve relating the naturalness index (NI) with the potential disappeared fraction of species (PDF) was developed using species richness data from the literature. LCA impact scores in PDF for producing 1 m 3 of wood might lead to consistent results with the naturalness index but the uncertainty is high while the window leading to consistent results is narrow.


Introduction
Green building relies on quantitative tools such as life cycle analysis (LCA) to evaluate the environmental impact of alternative design choices, among them the supply of different building materials [1], to support environmentally sound choices in the context of decisionmaking based on scientific evidence [2]. For harvested wood products, the impact on ecosystem quality is influenced by the forest management strategies, silvicultural practices and management intensity [3,4]. Impact assessment of anthropic activities on ecosystem quality in LCA has been a subject of research for more than 20 years [5,6]. According to the LCA framework, the impact corresponds to the marginal change in ecosystem quality related to a given land use during a certain period of time (considered here on an annual basis), multiplied by the area required to produce the functional unit [7], corresponding here to 1 m 3 of wood used in building structure. These impacts are assumed to grow linearly with the quantity produced [2], whereas land use intensity should be related to impact using non-linear functions [8].
Anthropic activities and especially habitat change related to land use and land use change has caused losses of biodiversity and contributed to an unprecedented rate of example, Turner et al., 2019 [5] found an LCA impact score lower for a forest management system based on pine plantations to produce softwood, compared with a native forestry system relying on partial cuttings, though the former was recognized as the one having the most negative impact on ecosystem quality. Another case using a SAR-based approach obtained results which were not always aligned with the existing scientific knowledge on biodiversity and forest management [24]. This situation is problematic as it can lead to inaccurate conclusions about management practices [24,33]. As stated by Michelsen et al., 2014 [34], "when different approaches give rise for different recommendations, the usefulness for decision-making will be questioned". Therefore, it is important to examine the passage from an assessment of the impact on the ecosystem quality over a given territory to the LCA impact evaluation, which calculates the impact per area required to produce one functional unit (i.e., 1 m 3 of wood for building structure).
We propose an approach to evaluate the potential impacts of forest management intensity on ecosystems based on the concept of naturalness. A naturalness assessment model has been designed for this purpose for boreal forests [35]. Many authors have proposed to use the concepts of naturalness and hemeroby to evaluate impacts of land use, such as forestry, on the quality of ecosystems in LCA [14,15,[36][37][38]. Naturalness is defined as "the similarity of a current ecosystem state to its natural state" [39], whereas hemeroby expresses "distance to nature" in landscape ecology [15]. The use of these concepts provides indexes to be used as guiding principles for forest management overcoming the challenge of gaps in obtaining empirical data on biodiversity. Côté et al., 2019's [35] naturalness assessment model relies on forest inventory maps and data to compare the anticipated outcome of different forest management strategies, including protection over a given territory [40].
This study aims at exploring the integration of the naturalness assessment results in a characterization model to evaluate ecosystem quality impacts of different forest management scenarios within the Life Cycle Analysis (LCA) framework, in order to verify its relevance for forest management decision-making. The specific objectives are to:

1.
Establish the relationship between naturalness and species richness and expand the model developed by Côté et al., 2019 [35] to express results in potential disappeared fraction of species (PDF); 2.
Evaluate the effects on naturalness and PDF of increasing the share of protected area combined with different silvicultural scenarios and compare different forest management strategies in the Québec boreal forest; 3.
Provide an example of LCA impact score calculation based on a naturalness assessment transformed in PDF to evaluate impacts on the quality of ecosystem of harvested wood products supplied from different forestry management practices in the LCA context.

Test Area
The test area is the same that has been used in Côté et al., 2019 [35] and corresponds to a public forest formed by three forest management units (FMU) located in the western black spruce feathermoss bioclimatic sub-domain, near the locality of Chibougamau in Northern Quebec (Appendix A: Figure A1).
The area arrangement considered for the 3FMU is presented in Figure 1. The share of protected area was applied on the available forested area (forested minus excluded), and practices related to wood procurement were applied over a given proportion of the productive area. Silvicultural scenario components must encompass 100% of the productive forested area. The excluded forested area [41] and the protected area were ascribed a natural evolution based on historical reference data (Appendix A: Table A1). Typology of the territory: silvicultural components of the scenario are applied over the productive area and protection corresponds to the portion of forest area not excluded, nor used for wood production. Specified areas are those from the test area: FMU no. 2663, 2665 and 2666 [41,42]. Scenario coding: 60CL100-40PL70-54.7p: 60% of the productive area in careful logging with a rotation of 100 years, with 40% of the productive area in plantation on a rotation of 70 years, combined with 54.7% of the forest area in strict protection © Sylvie Côté, Laval University, 2021.

Forest Management Scenarios and Strategies Tested
Three components were involved in the tested forest management scenarios: strict protection and silvicultural scenarios of careful logging and plantation. Strict protection was defined as protected areas belonging to IUCN category I, II or III [43]. The practices considered for silvicultural scenarios applied in the productive area were careful clearcut logging (CL), which corresponds to clearcut with the protection of pre-established regeneration and soil as required by law for clearcut operations in Quebec [44] and forest plantation of indigenous species mechanically released (PL). To determine results at the landscape level, each silvicultural scenario, corresponding to the sequence of silvicultural interventions over time, were applied on a constant basis over a given proportion of the productive area. However, this exercise is highly theoretical, considering that in reality, the different scenarios are adjusted through time. The impact of future fires was not taken into account. Ecological effects of careful clearcut logging and plantations are detailed in Appendix A.
Three silvicultural scenarios were considered: (1) CL only: 100% in careful logging on a 100-year rotation (coded 100CL100); (2) Current Mix: 60% CL on a 100-year rotation combined with 40% PL on a 70-year rotation (coded 60CL100-40PL70); (3) PL only: 100% plantation on a 70-year rotation (coded 100PL70). The two 100% scenarios evaluate the maximal theoretical effect of each practice. Forest management scenario coding (Figure 1) details the percentage of productive area of each silvicultural component (CL or PL) specifying the related rotation length (PL70: plantation with a 70-year rotation) completed by the percentage of protected forest area (e.g., 17.9p: 17.9% of the forest area in strict protection).
A forest management strategy corresponds to a plan of action or a policy designed to ensure a major aim, here defined as sustainable timber production. Two different forest management strategies were considered: (1) variable wood production over a fixed productive area, and (2) fixed wood production over a variable productive area completed with strict protection. The first strategy implies that each silvicultural scenario was applied over the current productive area, considering the current proportion in protection (17.9% of the forest area) and varying wood production. The second forest management strategy considers intensifying wood production to increase area available for conservation but aims to produce the same volume from each of the three tested silvicultural scenarios. The area available for protection was computed as follows: forested area minus excluded area minus productive area necessary to produce the expected volume ( Figure 1).
The effect of the level of protection was evaluated by the reassessment of the three silvicultural scenarios with an increasing proportion of forest area in protection, from 0% up to 90%, by steps of 5%. We used the proportion of forest area to control forest protection objectives so as to not take into account unproductive areas that occupy a significant proportion of the territory (18% of the terrestrial area). We assumed that conservation plans considered ecosystem irreplaceability and vulnerability [45].

Naturalness Assessment
The impact of forest management scenarios involving several forestry practices on ecosystem quality was evaluated using the Naturalness Assessment Model developed for the Picea mariana-feathermoss ecological domain of Quebec [35]. The naturalness assessment of forest management scenarios considered that the scenarios were being applied over the whole territory, as described in Côté et al., 2020 [40].
The Côté et al., 2019 model [35] evaluates five naturalness characteristics: landscape context, forest composition, structure, dead wood and regeneration process ( Table 1). The model uses indicators of condition (closed forests for context; cover type and late successional species groups for composition; old forests and irregular stands for structure), which are assessed against historical values, and pressures indicators resulting from practices (anthropization, modified wetlands, clearcuts on wetlands, exotic species, companion species, horizontal structure, dead wood and regeneration process). The evaluation is realized in two steps (see Appendix A for details). First, a partial naturalness index for each condition indicator (condition_pni: pni in lower case) is determined as a function of empirical curves relating measures (percentage of area) to condition_pni ( Figure A2) and naturalness degradation potentials (NDP) are evaluated through tables relating percentage of forest area by practices to NDP factors (Appendix A: Tables A2 and A3). Second, the partial naturalness for each naturalness characteristic (characteristic_PNI: PNI in capital letters) is calculated using the corresponding equation as per Table 1. The final result, the naturalness index (NI), is obtained from the arithmetic mean of the five characteristic_PNI. For the full description of the method, see Côté et al., 2019 [35].

Expanding Naturalness to Species Richness Assessment
In LCA, it is proposed to express the impact of land use on ecosystem quality by the potential disappeared fraction of species (PDF) [11]. The naturalness assessment model used in this study was expanded to estimate the anticipated effects on species richness using global data. To do so, the five naturalness classes (by NI interval of 0.2) based on the conceptual model defined in [35] were related to species richness data as per Figure 2, using data associated with land use characteristics corresponding to the pressure level (i.e., primary vegetation-minimal, light or intense use; secondary vegetation-young, intermediate or mature-minimal, or light/intense use; plantation forest-minimal, light or intense). Data were taken from the PREDICTS data base [46] and built on land use characteristics corresponding to different pressure levels based on 485 studies across the globe with available data of any taxon [3]. Species richness assessed in that study corresponds to the within-sample number of differently named taxa recorded at a given site. Values provided are expressed relatively to an un-impacted baseline (score of 100%), corresponding to complete naturalness. Species richness values from [3] were considered to be the biodiversity potential and correspond to (1-PDF). The PDF value for a given NI was computed using a linear estimation between established points of correspondence. The resulting curve ( Figure 2) is provisional considering the global scale, the restricted representation of boreal forests in the database used and the necessity to quantify empirically the relationship between naturalness and PDF. This curve was used to calculate an impact score for LCA based on PDF derived from naturalness (NI). The NI-PDF curve was given a concave shape ( Figure 2) showing a decreasing PDF with a rising naturalness index. This corresponds to a growing negative effect (i.e., more species loss) as the pressure resulting from forest management for wood production rises, which is consistent with cumulative impacts related to multiple stressors, and effects accumulation over large areas and over time [47].
We underline the high level of uncertainty of the NI-PDF curve resulting from the high variability of the available biodiversity data. Therefore, it should not be used as a basis for decision-making. However, we performed a sensitivity analysis of the curve setting using the mean, to explore what would be the effects of variations in corresponding values of PDF on the forest management scenarios assessment using LCA impact scores. The curve building procedure is detailed in Appendix B and the alternative curve settings used for sensitivity analysis are detailed in the Supplementary Materials.

Life Cycle Assessment Framework of Land Use
According to the land use impact assessment framework in LCA, the impact score (IS) on ecosystem quality due to the land occupation by a given forest management was calculated as follows [7]: where ∆Q corresponds to the change in ecosystem quality, here assessed against to the natural situation, A corresponds to the area of land occupied (ha) and T the time of occupation (years). A × T typically represent the area-time of occupation related to the functional unit, e.g., the ha × years required to produce 1 m 3 of harvested wood product. ∆Q corresponds to the characterization factor and is computed by a given impact assessment methodology with its specific category indicator units. The theoretical framework considers a sudden transformation of ecosystem quality from a land use to another [10]. In reality, the initial transformation related to forestry is gradual at the scale of the whole productive area and will span over the length of the first cutting cycle [35]. Land occupation for forestry scenario evaluation corresponds to the state of the ecosystem after the implementation of the production scenario over the entire productive area [40].
For our study, ∆Q was expressed as the potential disappeared fraction of species (PDF) (a dimensionless value between 0 and 1) and A × T was the terrestrial area required annually to produce 1 m 3 of wood (TAR) and corresponds to the life cycle inventory (LCI) data representing the flow for land occupation [7]. The TAR was based on yield hypothesis (see Appendix A for details) and the sensitivity of yield was evaluated by varying the respective scenario yields by ±20%, for plantation and careful logging. The area needed each year for the production of 1 m 3 was calculated over the terrestrial area as follows: where the yield related to each scenario component, expressed in (m 3 /(ha productive × yr)), where yr = 1, was multiplied by the productive area to calculate the annual harvest in (m 3 /yr), which was then divided by the terrestrial area of the 3FMU to calculate the annual harvest related to the total terrestrial area in (m 3 /(ha terrestrial × yr)). For LCA, we need the area required for the production of 1 functional unit defined here as 1 m 3 of wood.
Applying the rule of three, the inverse (1/x) of the annual harvest related to the terrestrial area corresponds to the TAR (terrestrial area required annually to produce 1 m 3 of wood) and is expressed in ((ha terrestrial × yr)/m 3 ), where yr = 1 and m 3 = 1. Therefore, the LCA impact score (IS) corresponds to: where yr = 1 and m 3 = 1.

Naturalness Assessment Results on the Test Area
Considering a fixed productive area with varying wood production, the CL only scenario results in the highest naturalness, followed by the Current Mix and the PL only (Table 2; supported by detailed results in Supplementary Materials, sheet "results_current_ protection"). In comparison, the PL only silvicultural scenario has a positive effect on composition, but negative effects for the four others naturalness characteristics. A forest management strategy aimed at producing a maximum of wood over a given productive area results in a low level of naturalness. With 17.9% of the forest area in protection, regeneration through careful logging would produce 502,100 m 3 /yr, leading to a naturalness index in the lower portion of the semi-natural class (NI = 0.470). With the same protection level, regeneration with indigenous species plantations would produce more wood (1,629,025 m 3 /yr) but lead to an altered naturalness state (NI = 0.301).
Considering a fixed production of wood (502,100 m 3 /yr), it would be possible to produce the same amount of wood with PL only over 28.3% of the forest area with the remaining 71.7% of the forest area in strict protection. Such a strategy would lead to a naturalness index in the higher portion of the near-natural class (NI = 0.746). Applying the current regeneration mix (60% careful logging and 40% plantation) over 45.3% of the forest area while protecting the remaining 54.7% would lead to a naturalness index in the lower portion of the near-natural class (NI = 0.649). According to the conceptual model of the naturalness assessment model, there is a high probability of species persistence associated with the near-natural class ( Figure 2).
The naturalness assessment results, expressed as a function of the level in strict protection, show that the level of naturalness rises with an increase in the proportion of protected forest area ( Figure 3; see detailed results in Supplementary Materials, sheet "NI-PDF_a"). The naturalness index of CL only without protection (NI = 0.414) is almost equivalent to the naturalness index resulting from PL only combined with 35% of the forest area in protection (NI = 0.408), or from the Current Mix combined with 20% of protection (NI = 0.417).

Transformation of NI to PDF
The transformation to PDF required for life cycle assessment has the following effects on the NI ( Figure 3): (1) it reverses tendencies to express damages (with NI representing the state of "what is remaining": PL only < Current Mix < CL only, while PDF represents "what is lost": PL only > Current Mix > CL only); (2) it narrows and reduces markedly the range of values (NI: 0.23-0.92; PDF: 0-0.16); (3) when increasing the protection level from 0% onwards, it reduces more rapidly the differences between silvicultural scenarios; and (4) it sets all effects at 0 when above a certain level of protection. We adjusted part of it by introducing infinitesimal values for PDF when NI lies in the near-natural class, but this does not affect the results and conclusions.

Life Cycle Assessment Results
Life cycle assessment results and related key input data are presented in Table 3 for the two forest management strategies, resulting from the combination of each of the three silvicultural scenarios assessed with the protection level required to achieve the objective of the strategy. Silvicultural scenario ranking results are summarized in Table 4 for each LCA parameter.   For the strategy over a fixed productive area considering 17.9% of the forest area protected, the scenario ranking based on the PDF is the opposite of the TAR-related ranking ( Table 4). The strategy producing annually a fixed volume of 502,100 m 3 of wood, with the CL only scenario, requires 651,610 ha of productive forest, which leaves 149,856 ha available for protection, corresponding to the current level of protection of 17.9% (Table 3). The same production of wood can be obtained with the Current Mix scenario combined with 54.7% of the forest area in protection or with the PL only scenario combined with 71.7% in protection. Therefore, producing the same volume of wood (502,100 m 3 ) over the same terrestrial area (1,072,294 ha) for the three management scenarios (each combining one silvicultural scenario with different percentage of area in strict protection) leads to the same value for the TAR.
For the strategy over a fixed area, we observe that the impact scores are highly correlated to the PDF characterization factors and lead to the same ranking of management scenarios (Table 4), but show less difference between the scenarios compared with the PDF (Table 3). For the fixed wood production strategy, it is not possible to discriminate between Current Mix and PL only based on PDF and LCA impact score results because of zero values.
The evaluation of silvicultural scenarios considering a protection gradient is illustrated in Figure 4. It presents the modeled relationships between the impact scores (IS) and the potential disappeared fraction of species (PDF), the terrestrial area annually required (TAR) to produce 1 m 3 of wood and the proportion of the forest area in protection. The PDF characterization factor is illustrated along the protection level using the solid chart bars associated with the left axis. The TAR is illustrated along the protection level using the empty chart bars associated with the right axis. The impact score, which corresponds to the multiplication of the PDF with the TAR, is illustrated with the lines associated with the left axis.
The results along the protection gradient show a nonlinear relationship between the level of protection and the TAR (Figure 4). The shape corresponds to the inverse function of the land productivity resulting in the area required to produce 1 m 3 of wood (as per Equation (2)). The TAR results display a wide range varying between 0.53 and 9.04 ha_ter × yr/m 3 for PL only, between 0.92 and 15.47 ha_ter × yr/m 3 for the Current Mix, and between 1.74 up to 29.4 ha ter × yr/m 3 for CL only (see Supplementary Materials for detailed results). The PDF results expressed as a function of the level in protection show a decrease in the PDF as the proportion of protected forest area in protection increases ( Figure 4). For CL only, the PDF results varies between 0 and 0.033, between 0 and 0.070 for the Current Mix, and between 0 and 0.1638 for PL only. The TAR and the PDF show opposite effects (Figure 4), as well as opposite scenario ranking ( Table 4). The impact score ranking corresponds to the PDF one.

Sensitivity Analysis
A sensitivity analysis of the values used to set up the provisional NI-PDF curve was performed by applying alternative correspondence values for PDF (indicated in red in Figure 5 and S1) to evaluate their effects on the resulting impact scores. The setting and detailed results for each alternative curve are included in the Supplementary Materials, the EXCEL file sheets correspond to the letter associated with the curve tested. The representations of these results displayed accordingly to Figure 4 are included in the Supplementary Materials, sheet "Figures_curves_a_to_f".
The analysis shows high sensitivity to the values used to generate the curve relating potential disappeared fraction of species (PDF) to the naturalness index (NI). A small change in correspondence values, comprised within the confidence interval of the PDF data, can affect not only the magnitude of differences between scenario impact scores (Figure 5b,c) but can lead to different results (Figure 5d,e). Based either on the naturalness index or on the PDF the scenario ranking in decreasing order of damage is PL only > Current Mix > CL only. This ranking is consistent with scientific knowledge on ecological impacts of silvicultural treatments. The use of the provisional curve (curve (a)) leads to impact scores displaying the same ranking. However, using the curve (d), with 0 and 5% levels of protection, the impact scores would rank PL only > CL only > Current Mix. Using the curve (e), the impact scores with a 0 to 20% protection level would rank CL only > PL only ≈ Current Mix. With more than 25% of protection, they would rank PL only > Current Mix > CL only (see Supplementary Materials, sheet "NI-PDF_e" for detailed results).  We also observe that the use of the impact score to compare the level of protection related to a scenario would lead to different results depending on the NI-PDF curve setting. For example, the protection level where losses reach near zero values when using curves (c) and (d) is different from the level obtained with the provisional curve.
The second parameter involved in the impact score calculation is the TAR associated to the scenario, which results from the division by the yield. The effect on the impact scores of a variation of ±20% of the scenarios yield for PL or CL are illustrated in Figure 6. Results are affected by the difference between the scenario's productivity. Larger differences between the basic scenarios can lead to questionable conclusions at low levels of protection (Figure 6b,c).

LCA Impact Assessment
The trial use of naturalness in LCA considering gradient of protection sheds light on the LCA model behavior. The two parameters involved in the LCA impact score calculation, the potentially disappeared fraction of species (PDF) and the terrestrial area required annually to produce 1 m 3 of wood (TAR), showed opposite non-linear effects that have to be properly put in balance in order to lead to impact score results consistent with existing scientific knowledge about the impacts of forest management. The transformation in PDF using the provisional curve provided a range of values and a distribution allowing to offset the effect of the TAR and led to impact scores correlated with the PDF rankings (i.e., PL only showing the greater impact followed by Current Mix and CL only). As the biodiversity data necessary to calibrate and confirm empirically the NI-PDF curve are not available for our study area, the use of the proposed NI-PDF curve could be hazardous with the current state of knowledge. The sensitivity analysis showed that the window for obtaining consistent results is narrow. There are risks of inaccurate results if data on species richness loss do not have the mathematical characteristics allowing to compensate for the effect of the area required (TAR), itself determined by the inverse of the yield function (1/x). When the quality indicator does not vary appropriately with the variation of the TAR, it can produce misleading conclusions. This supports raised doubts about the possibility to implement the non-linear relationship between LCI data related to land use productivity (e.g., land occupied/transformed in area × year) and impact on biodiversity (e.g., potentially disappeared fraction of species locally) in LCA models that typically scale impacts up or down using simple proportionality [8]. This also supports worries expressed by Teixeira et al., 2016 [8] about the use of LCA biodiversity models. It highlights the need to look beyond the characterization factor and the necessity to analyze the effect on the impact score of introducing an indicator into the LCA model, in regard to the effective protection of the safeguard subject (i.e., the ecosystem quality).
Turner et al., 2019 [5] also obtained diverging results from LCA impact scores when comparing with existing knowledge on the effects of a given type of forest management; they linked this with issues related to scaling/weighing, cumulative effects and with the reference benchmark. Our results confirm the relevance of scaling and weighing indices from ecological models for their inclusion into the LCA model as ecosystem quality indicator. When evaluating the environmental impact of anthropogenic production systems, LCA assumes a linear relationship between the impact and the quantity produced, which does not reflect the shape of cumulative effects [2]. However, our results indicate that the mathematical framework imposed by the TAR calls for the use of an indicator reflecting somehow the cumulative effects, because it has the power to compensate for TAR's preponderance within the LCA framework. Interestingly, most biodiversity indicators or weighing approaches proposed in the land use LCA context do reflect a cumulative effect. It is the case in our study with the potentially disappeared fraction of species along the alteration gradient (Figure 2), as it is often found in the literature expressed through the use of some kind of exponential/logarithmic relationship between habitat degradation and species losses. It can be seen in Chaudhary et al., 2015 [48] with the SAR, in Rossi et al., 2018 [38] between hemeroby level and partial biodiversity scores, in Lindqvist et al., 2016 [49] to express the contribution of a given parameter to a biodiversity potential function or in Turner et al., 2019 [5] for weighing results from expert opinions.
In the perspective of using LCA results for forest management decision-making, the present study showed difficulties for both the protection level and the silvicultural scenarios. The transformation of NI in PDF required for LCA involved values equal or near 0 for PDF in the higher portion of the protection gradient to offset the non-linear portion of the TAR. This impaired the possibility to properly consider the effects of the level of protection on the quality of ecosystems. As for the evaluation of silvicultural scenarios, the transformation in PDF reduced the difference between scenarios more rapidly along the protection gradient. The sensitivity analysis showed that the LCA impact scores could not be used as comparative basis to associate the level of protection related to a given management scenario as it would correspond to an artefact of the NI-PDF curve. Moreover, the sensitivity analysis showed that the LCA impact scores were sensitive to the PDF value and to productivity, as well, and there are risks of curve inter-crossings (i.e., changing ranking of silvicultural scenarios) which produce results that are difficult to interpret. Another point is that the use of an LCA impact score for decision-making in sustainable forest management has the potential to misinform as it does not take into account the crossing of ecological thresholds. An ecological threshold corresponds to a point where a small marginal change in the quality of ecosystems can cause a large response leading to an irreversible reorganization of the ecosystem [50,51]. LCA does not aim to evaluate the shifting state of ecosystems or biodiversity [2] and most of the time we do not know exactly which level corresponds to an ecological threshold until it has been crossed. Furthermore, for a given ecosystem, biodiversity thresholds can be numerous [52] and cannot be expressed on a functional unit basis. This situation supports the recommendation of using site-specific evaluations for decision-making about land use and management practices [24].
In any case, best management practices or strategy identification must be aligned with local evaluation over a given territory, as we are managing at the scale of hundreds of thousands of hectares of forests, which are finite over the planet. It would not be environmentally sound to take a direction towards more production at the expense of the overall ecosystem quality. Site-specific and/or territorial studies should be used to support decisions [24] in order to prevent inaccurate conclusions about the forest management strategy and practices and their effects on landscape [8].
The available LCA model calls for another round of critical evaluation of the underpinning of the land use aspects in LCA modelling, considering the mathematical features related to the non-linear relationships involved [48]. Moreover, questions still remain about the proper method to consider ecosystem quality in LCA. As recommended by Teixeira et al., 2016 [8], LCA modelling should strive for a better coupling between results and policy decisions and existing strategic plans.

Targeting the Level of Protection Using Naturalness
Limitations imposed to the incorporation of naturalness assessments into the LCA model does not question the use of naturalness assessment as a tool to assist decisionmaking related to forest management scenarios and strategies.
The area set aside for strict protection is an important component of a forest management strategy to prevent or limit biodiversity losses [30][31][32]. The question of how much is enough has been a source of questioning for years [53,54].
With the hypothesis that strictly protected forests have the characteristics of natural forests, the increase in protection contributes to the improvement of overall naturalness, and to the decrease in the potential loss of species due to forest management. Protection targets could be estimated using the relation between the naturalness index and the potential impacts on biodiversity ( Figure 2). For example, aiming a management target seeking to limit species loss (i.e., NI ≥ 0.5), the proportion in protection should be 25% of the forest area with the CL only scenario (Figure 3 and Supplementary Materials, sheet "NI-PDF_a"). It should be set at 35% with the Current Mix scenario (which represents twice the current level of protection for the forest management units (FMU) under study), and 50% for the PL only scenario. This level of protection is consistent with the level recently proposed as a "big bold conservation target" [32]. However, according to the naturalness conceptual model (Figure 2), if the objective is to avoid species loss, the naturalness level target should aim for the lower limit of the near-natural class (NI ≥ 0.6). Such a level is reached with 45% of protection with the CL only scenario, 50% of protection with the Current Mix scenario or 60% of protection with the PL only scenario (Figure 3 and Supplementary Materials, sheet "NI-PDF_a"). This 60% level corresponds to the proportion of area set aside estimated by Framstad et al., 2002in Michelsen 2008 to be necessary to maintain biodiversity within boreal forests. The current protection level for the three FMUs under study (17.9% of the study forest area) combined with the Current Mix of regeneration methods (60% in careful logging and 40% in plantation) is expected to lead to an ecosystem at the limit of being altered (NI = 0.407). This level is likely to produce some biodiversity losses. These losses could be limited by doubling the proportion of forest area in strict protection to reach a level of 35%.

Considering Protection in Forest Management Strategies Using Naturalness
A forest management strategy producing wood over a given productive area shows a lower naturalness level if wood is produced through plantation of indigenous species when compared with careful logging ( Table 2). In contrast, a forest management strategy producing a given amount of wood leads to a higher naturalness if wood is produced through plantation of indigenous species with a revolution of 70 years over a small portion of the area combined with strict protection of the remaining fraction, compared with production over the whole current productive area through careful logging with a revolution of 100 years combined with the actual level of protection. Based on the naturalness assessment results, it would be better for the environment to intensify forest management over a small proportion of the forest territory while ensuring strict protection over the remaining portion, compared with an extensive forest management applied over the vast majority of the forested area. This confirms the pertinence of a zoning approach such as the "Triad strategy" [55,56]. Obviously, this conclusion holds true only if the area not used for wood production is permanently set aside for strict protection through an effective strict protected area system [57], and if the enhanced production does not involve plantation of exotic species. However, the implementation of scenarios involving increased production to compensate for protection is complicated by the necessary delay to obtain wood products from plantations and the related question of wood procurement in the meantime [58,59].
Such approaches require a paradigm shift in policy-making [60] and necessitates prompt action to ensure conservation of primary forests.

Future Research
High priority should be given to studying the effects of forest management on biodiversity, especially species loss in boreal forests. This should be put in relation with research on the relationship between naturalness and biodiversity. The reliability of the NI-PDF curve set up could be improved using an updated subset of the PREDICTS database corresponding to boreal forest with the appropriate definition of land use intensity. Further empirical data should be sought to reduce uncertainty and confirm the relevance of using species richness data in the LCA model.
The use of the naturalness concept could be promoted, by adapting the model for more bioclimatic domains and performing assessment over more forest management units (FMU) and at the bioclimatic domain scale. To expand our approach to ecosystem quality evaluation in LCA, the method used for naturalness assessment of forestry should be adapted for other land uses in order to allow land-use comparisons based on the ecosystem quality. To do so, the naturalness assessment should be integrated into the whole alteration gradient [35], along with the other land uses by designing a hemeroby assessment model based on the habitat characteristics.
To improve the robustness of the LCA model, it is of outmost importance to investigate the mathematical effects of non-linear parameters in order to verify if and how issues related to scaling and non-linear effects in the LCA model could be properly addressed.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Naturalness Assessment of the 3 FMU and Related Data
The test area corresponds to three Forest Management Units (FMU) (nos. 2663, 2665 and 2666) covering a total area of 1,305,200 ha [41]. The localization of the test area is shown in Figure A1. The study area is currently subject to the first cycle of harvest. Figure A1. Test area localization.

Appendix A.1. Historical Data of the Test Area
Historical data are shown in Table A1. Reference data for cover type (CT), late successional species groups (LS), old forests (OF) and irregular stands (IR) were taken from the first Quebec forest inventory, corresponding to the period from 1965 to 1974, using Quebec's SIFORT system (tessellation of provincial forest inventory maps) [61], from which 6% of harvested areas and other anthropic disturbances were removed. To conform to the integral application of the forest management scenario, the reference value for CF (closed forest: >40 years old) corresponds to the (total forest area minus the forest area in old forests), spread evenly over the younger age classes (0-100), then adding the corresponding area for 0 to 40 years.

Appendix A.2. Ecological Effects of Careful Logging and Plantations
Careful clearcut logging represents a significant and abrupt change in habitat conditions for forest wildlife [62]. At the stand level, clearcutting induces modifications in vegetation structure and composition, microclimate and ground substrate [62][63][64]. It also resets the vegetation succession to earlier stages [65]. At the landscape level, even aged management (regenerated through clearcuts) modifies the age structure (distribution of age classes) of the forest by reducing the proportion of old forests to an extent related to the rotation length [63]. This precludes the formation of characteristics related to old growth forests [66], thus affecting species associated with old-growth habitat [60]. Impacts on biodiversity are mainly influenced by the rate and the spatial pattern of harvesting and the type and intensity of subsequent management practices [63].
Plantation corresponds to an intensive form of forest management involving a suite of anthropogenic manipulations of forest regeneration including site preparation, tree planting (including sometimes genetic selection of the planted trees) and various treereleasing and tree-thinning treatments [64]. Plantation forests support generally lower diversity compared with natural forests [64,67]. Negative effects are particularly important when exotic species are used [29,58,67], which is not the case in the study area. Plantation of indigenous species can have positive effects by contributing to forest composition restoration, for example [68], but still there are potentially negative effects related to modification of genetic pools, site preparation, shorter rotation, repeated thinning and homogenization of conditions [29,66,67]. Cumulative effects on biodiversity of intensive forest management are likely to increase over time as a result of reduced species diversity, reduced amount of dead wood and enhanced frequency of high disturbance over large areas [66]. By producing more wood from intensively managed areas, plantations could be used to reduce indirectly pressure on natural forests [58,67,68]. To be efficient, such an approach requires that forest conservation objectives are permanent and guaranteed by rigorous legislation [60,66], which can be obtained with the creation of strictly protected areas [43].

Appendix A.3. Scenario Naturalness Assessment
The territory used for the analysis covers the whole area included in the perimeter of the FMU (without cutting tessell in SIFORT maps), including the surrounding strict protected areas (IUCN categories I to III) associated with these units. The percentage of forested area over the territory of analysis was calculated for measures of forest condition (CT, LS, OF and IR) and the percentage of terrestrial area over the territory of analysis for context measures (CF, W_CC, Wm, ant).
Age structure related to each scenario under sustainable production has been used to evaluate closed and old forests. For example, for a rotation of 100 years, 10% of the productive area will be in each 10 years class (0-10 years up to 90-100 years), and for a rotation of 70 years, 14.2857% of the productive area will be in each 10 years class (from 0-10 years up to 60-70 years).
Current data, corresponding to the 2011 to 2013 period, were taken from the fourth inventory program [42]. Clearcuts on wetlands and anthropization levels were set to current values using data from ecoforest4 map [42] (% of terrestrial area: W_CC = 7.88%, ant = 0.74%) and kept constant in all scenarios. Percentages of forested area by origin and by age class considering silvicultural treatments has been used to determine CT and LS under 20 years and over 20 years old, in careful logging (CL) and plantations (PL). Hypotheses used for condition indicators (CT, LS, IR and OF) are shown in Table A2 (see Côté et al., 2019 [35] for details). Values for CF (closed forests over 40 years old) are determined based on age structure resulting from scenario components application. The natural degradation potential (NDP) factors used are shown in Table A3 and correspond to those used in Côté et al., 2019 [35]. The curves elaborated for pni evaluation specifically for the three FMUs are presented in Figure A2.
The naturalness assessment of the Current Mix scenario forecasting 60% of CL on a 100-year rotation with 40% PL on a 70-year rotation with the current level of 17.9% of the forest area in protection (60CL100-40PL70-17.9p) is provided as an example for the detailed description of the naturalness assessment method (File: SM_Results; Natural-ness_assess_ex.xls).

Appendix A.4. Scenario Yield Hypotheses
For the protection, obviously there was no yield. For the careful logging, the yield estimation of 0.77 m 3 /(ha productive × yr), was based on the sustainable harvest calculated in 2014 [41] (502,100 m 3 /yr) over the productive area (651,610 ha), considering a rotation of 100 years. For the plantations, an expected yield of 2.5 m 3 /(ha productive× yr) was based on tables for black spruce plantations having a station quality index of 6 (IQS = 6) with 2000 or 2500 stems/ha has been used [69]. For the 60% careful logging with 40% plantation scenario, the yield corresponds to the weighed sum of the productivity of both scenarios (1.462 m 3 /ha productive × yr).
The main source of uncertainty related to the estimation of the yield is for plantations, as no plantation has ever reached maturity in the study area and the plantation yield is more likely to be overestimated due to unaccounted mortality over the lifespan of the forest plantations [70]. Moreover, for the PL only scenario, the mean expected yield should be lower as plantations covering 100% of the productive area should be lower than the productivity of plantations concentrated on the best sites.

Appendix B. NI-PDF Curve Setting
The points of correspondence used to build the NI-PDF curve were established as follow. According to the naturalness model's conceptual framework [35], there is no specie loss (i.e., PDF = 0) in the natural class (NI = 1-0.8), and there is a very low probability of losses (i.e., PDF ≈ 0) in the near-natural class (NI = 0.6-0.8). The PDF has been set to 0.0000001 for NI = 0.7 and to 0.00001 for NI = 0.6 to consider a risk of marginal losses of sensitive species while allowing the detection of improvements in the near-natural class. The losses of very sensitive species begin in the semi-natural class (NI = 0.4-0.6), but PDF are still limited. In another territory located in a balsam fir forest, a naturalness index around 0.5 could be put in relation with decline of few sensitive species [40]. Unfortunately, these impacts have not yet been fully quantified, but they indicate that the PDF should be different from 0 when the PDF reaches that level. The PDF has been set to 0.001 for NI = 0.5 and alternatives values were tested thru sensitivity analysis. The correspondence values for NI ≥ 0.5 represents expert opinions proposed by the authors and have to be validated. For the near-natural class (NI 0.6-0.8), we introduced infinitesimal values for PDF to verify the possible effects of such introduction in the model, in the perspective of opening the reflection surrounding the consideration of effects on biodiversity when there are no species losses. However, the results for that portion on the protection gradient were not adjusted accordingly in the text (results rounded to the third digit), nor used in for interpretation. The results without the introduction of the infinitesimal section for the near-natural class are provided in the sensitivity analysis (curve f). To link PDF to NI, the age structure (area by age classes) resulting from the application of the 100% scenarios over the whole area without protection (excluded portion not considered) has been used to estimate the proportion of the area in each development stage of secondary forests used by Newbold et al., 2015 [3] (Table A4). The naturalness evaluation for the scenario involving CL only without protection (100CL100-0p) results in NI = 0.414 (see Supplementary Materials, results using the provisional curve with 0% protection), which can be roughly associated with the lower bound of the semi-natural class (NI = 0.4). Gradual application of clearcuts over the entire productive area will lead to secondary forests, spread over the various development stages depending on their age (young: 0-20 years, so 20% of the area subject to CL only; intermediate: 21-70 years, so 50% of the area in CL only; mature: 71-100 years, so 30% of the area in CL only). Using the species losses ((100 species richness data) ÷ 100) from Newbold et al., 2015 [3] for secondary forests lightly and intensively used (Table A4), weighted by the proportion in each development stage associated with their age, gives a mean of 0.038. Therefore, we associated the lower bound of the semi-natural class (NI = 0.4) with a PDF of 0.038. Losses of sensitive species are measurable in the altered class (NI = 0.2-0.4) and important losses of species can be associated with the very altered class (NI < 0.2). The naturalness evaluation for the scenario involving PL only of indigenous species without protection (100PL70-0p) gives a result of NI = 0.231 (see Supplementary Materials, results using the provisional curve with 0% protection), which can be roughly associated with the lower bound of the altered class. Therefore, we associated the lower bound of the altered class (NI = 0.2), to a PDF of 0.192 using the species richness data for plantation with minimal use from Newbold et al., 2015 [3], considering that plantations are of native species, not fertilized nor cleared with herbicides. The PDF at NI = 0.3 has been set to 0.1 to smooth the curve. Finally, tests of extremes hypothetical scenarios involving exotic species lead to a naturalness index lying in the very altered class [35]. We associated the middle of the very altered class (NI = 0.1), as no forestry use will ever create a biological desert, with the species richness data for intensively used plantation (PDF = 0.394). The last part of the curve (below NI = 0.2) has no effect in this study, since no scenario in the context of forestry modeled in this study is likely to lead to such degree of alteration. The levels used to estimate the end of the curve, for the altered and very altered levels, are consistent with the proportion of species associated with dead wood which are recognized to be sensitive to forest management. According to Klamerus-Iwan, et al., 2020 [71]: "It is now estimated that 20-40% of organisms in forested ecosystems depend, during some part of their life cycle, on wounded or decaying woody material from living, weakened, or dead trees". Species richness data for forests plantations intensively used equals 60.6 [3] (corresponding to a loss of 39.4% of species), which is very close to the maximal loss of 40% associated with forestry [71]. In Finland, where intensive forestry using indigenous species is practiced, nearly 20% of all species are extinct, threated or near threated [28]. This value is also very close to the PDF of 0.192 used to estimate the curve at NI = 0.2. Alternative correspondence values are tested in the sensitivity analysis.