Bryophyte Conservation in Managed Boreal Landscapes: Fourteen-Year Impacts of Partial Cuts on Epixylic Bryophytes

Clear cut harvest simplifies and eliminates old growth forest structure, negatively impacting biodiversity. Partial cut harvest has been hypothesized (1) to have less impact on biodiversity than clear cut harvest, and (2) to encourage old growth forest structures. Long-term studies are required to test this hypothesis as most studies are conducted soon after harvest. Using epixylic bryophytes as indicators, this study addresses this knowledge gap. Fourteen years after harvest, we examined changes in epixylic bryophyte community composition richness and traits, and their microhabitats (coarse woody debris characteristics and microclimate) along an unharvested, partial cuts and clear cuts harvest treatment in 30 permanent plots established in the boreal black spruce (Picea mariana) forests of northwestern Quebec, Canada. Our results were compared to those of an initial post-harvest study (year 5) and to a chronosequence of old growth forests to examine species changes over time and the similarity of bryophyte communities in partial cut and old growth forests. Coarse woody debris (CWD) volume by decay class varied among harvest treatments with partial cuts and clear cuts recording lower volumes of early decay CWD. The epixylic community was richer in partial cuts than in mature unharvested forests and clear cuts. In addition, species richness and overall abundance doubled in partial and clear cuts between years 5 and 14. Species composition also differed among treatments between years 5 and 14. Furthermore, conditions in partial cut stands supported small, drought sensitive, and old growth confined species that are threatened by conditions in clear cut stands. Lastly, over time, species composition in partial cuts became more similar to old growth forests. Partial cuts reduced harvest impacts by continuing to provide favorable microhabitat conditions that support epixylic bryophytes. Also, partial cut harvest has the potential to encourage old growth species assemblages, which has been a major concern for biodiversity conservation in managed forest landscapes. Our findings support the promotion of partial cut harvest as an effective strategy to achieve species and habitat conservation goals.


INTRODUCTION
Over the past two decades, partial cut harvest has gained global research interest due to its proposed biodiversity benefits compared to clear cut harvest (Franklin et al., 1997;Lindenmayer et al., 2012). Residual stands with variable retention of merchantable stems left after partial cut harvest are hypothesized to serve as refugia and enhance species recovery after harvest (Franklin, 1993;Drever et al., 2006). Also, residual stands could ensure a continual supply of deadwood, promote the development of old growth structural features, and the presence of old growth adapted species, which is currently not achieved with clear cut harvest (Harvey et al., 2002;Lilja et al., 2005).
Consequently, many studies have examined the effects of partial cut harvest on biodiversity (Baker and Read, 2011;Fenton et al., 2013;Mori and Kitagawa, 2014;Kuuluvainen et al., 2019;Fačkovcová et al., 2019). Numerous short-term studies have demonstrated the ability of partial cuts to generate favorable habitat conditions that support many forest species (Vanderwel et al., 2007;Caners et al., 2013;Fenton et al., 2013;Baker et al., 2016). This has ignited a lot of interest in long-term studies to better understand species recovery after partial cut harvest (Baker et al., 2015;Fedrowitz et al., 2014), as longitudinal studies in permanent plots can provide strong evidence of species recovery (Dynesius, 2015). However, very few long-term partial cut harvest studies have resurveyed species assemblages in the same sites (Halpern et al., 2012;Bartels et al., 2018;Wu et al., 2020). Furthermore, given that the loss of old growth habitat and species is a major conservation concern (Vaillancourt et al., 2009), it is important to examine the ability of partial cut harvest to encourage old growth species assemblages. However, except for a few short-term studies (Man et al., 2008;Boudreault et al., 2013;Bose et al., 2015), no long-term studies have assessed this.
Bryophytes constitute an important component of forest plant diversity, playing vital ecological roles as they: influence the nutrient cycle, provide habitat for micro-organisms and invertebrates and seed beds for tree regeneration (Do¨bbeler, 1997;Glime, 2011;Rudolphi and Gustafsson, 2011). The poikilohydric nature of bryophytes make many forest bryophytes highly sensitive to reduced humidity, thus, they are very sensitive to forest harvest, which alters forest microclimate through canopy removal (Fenton et al., 2003;Hylander, 2005). Clear cut harvest has been shown to have a long-lasting impact on bryophytes, particularly epixylic bryophytes (deadwood-living bryophytes), as epixylic bryophytes require a constant supply of decaying wood, and a humid and shady microclimate, which are lacking in regenerating clear cut stands (Söderström, 1988b). The unique physiology and habitat requirements of epixylic bryophytes make them well suited to study the ability of partial cut harvest to support bryophyte recovery, and to encourage old growth species assemblages. While some studies (Arseneault et al., 2012) have shown richer epixylic bryophyte communities in regenerating partial cut stands compared to clear cut stands, others (Halpern et al., 2012;Vanha-Majamaa et al., 2017) have shown little or no difference even 10 years after harvest. Further long-term studies are required to enhance our understanding on epixylic bryophyte recoveries in partial cut stands.
This study examined the efficacy of partial cut harvest to maintain epixylic bryophytes 14 years after harvest. We resurveyed epixylic bryophytes and their microhabitats in permanent plots established in the boreal black spruce (Picea mariana) forests of northwestern Quebec, Canada. Specifically, we (1) Examined changes in microhabitat conditions (CWD characteristics and microclimate) along an unharvested, partial cut, clear cut harvest treatments; (2) Examined epixylic bryophyte species composition, richness and functional traits along the harvest treatments. Also, we compared our results (year 14) to that of an initial post-harvest study (year 5) to examine species changes over time; and (3) Examined partial cut ability to encourage old growth species assemblages by comparing species composition in years 5 and 14 to that of a chronosequence of old growth forests. We hypothesized that (1) CWD abundance and volume in partial cuts will be higher in year 14 than in year 5, whereas CWD abundance and volume will decrease in clear cuts between years 5 and 14; (2) Year 14 epixylic bryophyte community diversity will be higher in partial cuts but not in clear cuts compared to the results of year 5; (3) Partial cut epixylic communities will be similar to communities of unharvested and old growth forest stands.

Study Location and Design
This study took place in northwestern Québec. Fire constitutes the main disturbance type and black spruce (Picea mariana) dominates the forest (Bescond et al., 2011). The landscape is characterized by low topographic relief and is prone to paludification as a result of poorly drained clay soils. Mean annual temperature is 0.7 • C, and mean annual precipitation is 890 mm (Environment Canada, 2019). We studied two sites (Fenelon and Gaudet) in an experimental partial harvest network. The average age of the study sites before harvest was more than 180 years since the last stand-replacing fire. Harvests were carried out in 2003 and 2004, an initial study was completed 5-6 years post-harvest (Arseneault et al., 2012), and this study was conducted 13-14 years post-harvest (in 2017). Each site is made up of one block (each > 50 ha) of three treatment types: clear cut harvest that preserves advance regeneration but removes all merchantable stems (DBH > 9 cm), partial cut harvest with variable retention and an unharvested control. Within each silvicultural treatment, 17,400 m 2 permanent plots were randomly established.
A total of 30 permanent plots were selected for this study, 10 per harvesting treatment, divided among the two sites (Supplementary Figure 1). We used the same ten plots per clear cut and control block as used by the initial post-harvest study (year 5), except in the control plots in one site (Fenelon) because the plots had been harvested. Four unaffected Fenelon control plots that were not selected in the 2009 study were used for our study (year 14). In the case of partial cut blocks, only plots considered a silvicultural success (i.e., there was a net growth in residual tree volume ten years after harvest; Moussaoui et al., 2020) were selected. This is because sites with a high residual tree mortality rate will decline in merchantable tree volume, making them silviculturally unsuccessful and hence not recommended for forestry practice as residual stands are intended for harvest in the next partial cut rotation. Two-thirds (67%) of the total partial cut plots were considered silviculturally successful whiles onethird (33%) were silviculturally unsuccessful (Moussaoui et al., 2020). Eight of ten partial cut plots in both sites differed from those of the year 5 study because they did not qualify. Seven of the ten selected partial cut plots were below 25% harvest retention while the remaining three were 30, 45, and 70% harvest retention. Overall, twelve plots in year 14 out of the total thirty plots differed from that of year 5. Four control plots and eight partial cut plots.

Microhabitat and Bryophyte Sampling
Our sampling method mostly followed the initial post-harvest study (Arseneault et al., 2012). We measured substrate (CWD) characteristics and canopy openness within each plot (11.28 m radius, 400 m 2 ). CWD characteristics included abundance, volume, decay class, diameter size and length, measured along a 23 m line transect in the north-south direction (Figure 1). The length and diameter of each piece of CWD were measured at the point of intersection along with tree species (where possible) and decay class. Decay class was estimated by observing the physical characteristics of CWD and was categorized from 1 (fresh material) to 5 (well decomposed) based on Hunter (1990). CWD abundance was calculated per count on the line transect for each plot (Marshall et al., 2000), and CWD total and by decay class volume for each plot was calculated using Van Wagner (1978) formula.
The first five CWD with a minimum diameter of 9 cm on each line transect (Andersson and Hytteborn, 1991) were selected for epixylic bryophyte community sampling. For each selected CWD, the total length, minimum and maximum diameter, percent of the length directly in contact with the ground, and maximum distance from the ground were measured. Belts of quadrats were randomly established at three different positions (ends and middle) on each selected CWD for bryophyte sampling. At each belt, diameter and decomposition class were evaluated to determine the intra-piece variability of decay class and diameter size. Bryophytes were sampled (Supplementary Figure 1) at each belt in three quadrats (totaling nine quadrats per CWD and forty-five quadrats per plot) at systematically placed positions: the top of the CWD (5 cm × 10 cm), the side of the CWD (5 cm × 10 cm), and the ground directly next to the CWD (10 × 10 cm). Larger quadrats were used on the forest floor because of the bulky nature of forest floor species compared to species growing on logs. In addition, this ensured the same amount of sampled area on CWD and forest floor (Arseneault et al., 2012). We also examined the remainder of the CWD to ensure capture of species that were not found in the quadrats. Species were identified in the field when possible, all other species were sampled and brought to the laboratory for identification. Nomenclature followed Faubert (2012Faubert ( , 2013Faubert ( , 2014. In terms of canopy openness, unlike the initial post-harvest study, which used a densiometer, hemispheric photos were used to measure canopy openness. Hemispheric photos were taken over each 5 × 10 cm quadrat on top of each CWD. Photos were taken on a cloudy day at the height of CWD using a circular fisheye camera lens at each quadrat belt. Images were analyzed using Gap Light Analyzer (GLA Ver. 2.0) (Frazer et al., 1999).
FIGURE 1 | Coarse woody debris volume by decay class per harvest treatment at the plot level (n = 30). Control (CT), Partial cut (PC), Clear cut (CC) and Decay class 1 to 5 (D1, D2, D3, D4, and D5). The dark horizontal lines in the box and whisker plots represent the median, with the box representing the 25 and 75th percentile and the whiskers the 5 and 95th percentiles and the dots represents outliers.

Chronosequence of Old Growth Forests
To test the ability of partial cut harvest to advance forest succession to a later stage, we took advantage of a chronosequence dataset compiled by Barbé et al. (2017) in old growth forests located in northwestern Quebec. These are unharvested and unmanaged forests that have remained unburned for 50-200 years as determined by dendrochronology. In the chronosequence, bryophytes were sampled on all available microhabitats in three 50 m 2 plots in 44 sites (Chaieb et al., 2015;Barbé et al., 2017). From this, a data set containing only the coarse woody debris microhabitats was compiled for our study.

Statistical Analysis
To address objective 1, microhabitat variables [substrate (CWD) characteristics and canopy openness] were compared among treatment types. CWD volume by decay class and total volume per plot were analyzed using linear mixed models [LME; R package "nlme" (v. 3.1-131), R Core Team, 2017], while CWD abundance at the plot level was analyzed using the Kruskall-Wallis test as transformed data were skewed. CWD diameter, length and decay class were analyzed at the log level using LME models. All data measured at the belt level (CWD characteristics, canopy openness) were averaged for a mean measure per CWD. For CWD decomposition class, the class closest to the mean was used. Canopy openness data were square root transformed and analyzed at the log level using LME models.
In all LME models, plots nested within sites were used as random factors. Significant differences between treatments were explored using Turkey HSD [R package "multcomp" (v. 1.4-1), R Core Team, 2017] tests for LME models and the dunn test (R package "dunn.test" v. 1.3.5, R Core Team, 2017) for the Kruskall-Wallis test with significance levels of p ≤ 0.05.
For objective 2, epixylic community changes among harvest treatments were examined in terms of species richness, composition and functional traits. For richness, species occurring more than once in the three belts of a CWD were considered only once in total richness. Species richness on CWD was based on only the quadrats that were on the CWD. Species richness data were analyzed at the CWD level using LME models. We assessed species frequency by counting the number of occurrences of each species in each harvest treatment (i.e., the overall abundance of epixylic species). LME models were developed to assess the relative importance of different habitat variables in explaining species richness. Microhabitat variables were divided in two groups: substrate characteristics and canopy structure. Substrate characteristics were CWD length, diameter, decay and percent contact with ground. Canopy structure variables were canopy openness and harvest treatment type. All environmental variables had an inflation factor [VIF (< 10)] and were included in the models. In total, 19 candidate models were developed including 1 global model, 6 models for each variable category alone, 11 models that tested biologically relevant combinations of different categories of variables, and the null model (Supplementary Table 1). Models were ranked based on the Akaike's Information Criterion corrected for small sample size (AICc) using the R package AICcmodavg (v. 2.1-1). Model averaging was carried out for models with AICc weight < 2.0 (Burnham et al., 2011). We then calculated average parameter estimates, and standard errors for top ranking models. R 2 for each variable in the best models was calculated to evaluate their relative contribution in the models, using the function r.squaredGLMM of the MuMIn package in R (v.1.42.1).
Furthermore, a comparative analysis was carried out with the data set from the post-harvest study (Arseneault et al., 2012) to assess bryophyte community changes with increasing time after harvesting disturbance. As above, species data on CWD was based on only the quadrats that were on the CWD. A ttest (R package, R Core Team, 2017) was used to compare epixylic species richness between treatments of the two studies. NMDS was used to visualize patterns in the species compositional data. A PERMANOVA was conducted to determine whether the species composition differed among treatments of year 14 and across time (i.e., between treatments in years 5 and 14). We used the adonis function in the vegan package (Oksanen et al., 2015) with Bray-Curtis dissimilarities, and 999 permutation runs. Plot in site was specified as a random effect factor to take into account the nested structure of our data.
To test partial cut ability to recreate old growth species assemblages for objective 3, we compared partial cuts species composition in years 5 and 14 to that of a chronosequence of old growth forests. We compared the year 5 data to the chronosequence with a PERMANOVA and used Bray-Curtis dissimilarities with 999 permutation runs. A similar analysis was then run comparing year 14 data to the chronosequence. NMDS was used to visualize species compositional patterns.

Effects of Harvest Treatments on Microhabitat
Mean number of CWD per plot did not differ among harvest types, however, it was reduced in year 14 compared to year 5 ( Table 1). Mean decay class values in year 14 were significantly higher in partial cut and clear-cut stands than in control stands. CWD volume did not differ among the harvest types, however, CWD volume of decay class 1 and 2 were significantly lower (or absent) in partial cut and clear-cut stands compared to control stands (Figure 1)

Species Richness Change Among the Harvest Treatment Overtime
Overall, 59 species were recorded in this study, 27 true mosses, 24 liverworts, and 8 sphagna (Supplementary Table 3). By harvest type, 44 species (20 true mosses, 17 liverworts, and 7 sphagna) were found in unharvested control stands, 50 species (25 true mosses, 21 liverworts, and 4 sphagna) in partial cut stands, and 39 species (18 true mosses, 14 liverworts, and 7 sphagna) in clear cut stands. Of these, three species were observed only in control stands, eight only in partial cut stands, and three only in clear cut stands. Eighteen epixylic species were found including 13 liverworts and 5 true mosses, with 12, 17, and 13 epixylic species found in unharvested control, partial cut and clear-cut stands, respectively. Four epixylic species were observed only in partial cut stands (Calypogea neesiana, Chiloscyphus profundus, and Lophozia longidens), and one only in control stands (Cephaloziella spinigera).
Epixylic species richness and total species richness were significantly higher in partial cut stands compared to control and clear-cut stands. However, no significant differences were observed between control and clear-cut stands (Figure 2). Compared to the post-harvest study, epixylic species richness (Figure 2) and overall species abundance (Supplementary Figure 2) was significantly higher in year 14 in partial cut and clear-cut stands than in year 5. Overall, total species richness did not differ significantly between the two years in control stands however partial cut and clear-cut stands were richer in year 14 compared to year 5 (Figure 2). In addition, epixylic species occurred mostly on CWD in year 14 compared to year 5, where most epixylic species occurred on small pieces of wood on the forest floor (Supplementary Figure 3).
For species functional traits, liverworts, mosses, generalists and terricolous species were richer in control and partial cut stands compared to clear cut stands but did not differ significantly from control stands in the case of liverworts and generalist species (Supplementary Figure 4). Richness of medium and smaller sized species were higher in partial cut stands than in control and clear cut stands but did not differ significantly from control stands in the case of small sized species (Supplementary  Figure 4). In contrast, large sized and bog species were significantly richer in control stands compared to partial cut and clear cut stands.
Microhabitat variables influenced by forest harvest (Table 1) were the drivers of the differences in epixylic richness found among treatments in year 14. Akaike model selection (Supplementary Table 1) indicated that two models were equivalent and model-averaged estimates ( Table 2) of these models indicated that epixylic richness was influenced by canopy openness, CWD decay class and diameter size, as was found in year 5. CWD decay class and diameter size positively influenced and canopy openness negatively influenced epixylic richness. Harvest treatment also affected epixylic richness indicating that other variables associated with harvest were important that were not accounted for in our models.

Change in Species Composition Among the Harvest Treatment and Over Time After Harvesting
PERMANOVA results indicated differences in species composition among harvest treatments of year 14 (F = 3.6876, R 2 = 0.0478, p < 0.001). The first NMDS axis (Figure 3) separated species by substrate type with epixylic species (Crossocalyx hellerianus, Dicranum fuscescens, Ptilidium pulcherrimum) on the right, generalist species (Aulacomnium palustris, Pleurozium schreberi, Ptilidium ciliare) and terricolous species (Rhytidiadelphus triquetrus, Dicranum polysetum) in the center and bog species (Sphagnum fuscum, S. magellanicum, and S. russowii) on the left. Control stands were dominated by bog species, partial cut stands by epixylic species and clear-cut stands by generalist and terricolous species. The second NMDS axis separated species along a moisture gradient. Species common in wet areas (Sphagnum magellanicum, S. russowii) on the top, moist areas in control stands; (Frullania oakesiana, Nowellia curvifolia) in the center in partial cut stands, and dry areas below (Pleurozium schreberi, Ptilidium ciliare) in clear cut stands. The relationship between microhabitat variables and the species composition pattern was also explored by passively correlating the variables with the plot positions ( Figure 3B). The first axis FIGURE 2 | Mean species richness per harvest treatment in years 5 and 14 (n = 30 plots). Significant differences (p ≤ 0.05) are indicated by different letters following the ranking (a > b > c). Different upper-case letters indicate significant differences between the same treatment of the two studies; different lower-case letters indicate significant differences between all treatments of the same year.
was correlated with CWD characteristics (length, decay class and diameter size) and the second axis was correlated with canopy openness.
Furthermore, PERMANOVA results indicated differences in species composition between treatments of years 5 and 14: clear cut stands (F = 4.6387, R 2 = 0.0548, p < 0.001), partial cut stands F = 3.6839, R 2 = 0.0402, p < 0.001), and control stands (F = 2.8942, R 2 = 0.0318, p < 0.002). The NMDS graph ( Figure 3D) showed a shift in species composition of clear cut stands from sphagnum and/or bog species in year 5 to generalist and terricolous species in year 14. Also, species composition in partial cut stands shifted toward epixylic species in year 14. However, species composition overlapped in the control stands of both studies.

Old Growth Species Assemblages
To assess partial cut ability to encourage old growth species assemblages, we compared CWD species composition in partial cut stands of both studies to a chronosequence of old growth forests. PERMANOVA results indicated no significant differences in species composition between partial cut stands and chronosequence stands: year 5 partial cut vs. chronosequence (F = 3.9774, R 2 = 0.0889, p < 1) and year 14 partial cuts vs. chronosequence (F = 3.7758, R 2 = 0.0803, p < 1). The NMDS graph did not provide a clear temporal pattern of species composition changes on logs from the chronosequence data. Species composition in year 5 partial cut stands was separated from that of old growth forests on the horizontal axis while year 14 partial cut species composition overlapped with that of old growth forest (Figure 4).

DISCUSSION
Our results demonstrated that partial cut harvest supports higher diversity and encourages old growth species assemblages. This study examined temporal trends in epixylic species richness in partial and clear-cut forests and found that over the time span covered by this study, partial cut stands continued to provide favorable habitat conditions that maintained higher epixylic richness compared to control and clearcut stands. Species composition in partial cut stands was more similar to that found in old growth forests in year 14 compared to year 5. As the majority of our plots had less than 25% retention, our results suggest that even low levels of retention in partial harvests can enhance postharvest recovery of bryophytes compared to clear-cuts. These results agree with previous studies (Arseneault et al., 2012;Bartels et al., 2018) but disagree with other studies that reported slow recovery of bryophytes in partial cut stands with species composition comparable to that of clear-cut stands even ten years after harvest (Halpern et al., 2012;Vanha-Majamaa et al., 2017).
Our first hypothesis on changes in CWD availability over time was supported for clear cut stands and refuted for partial cut stands as CWD abundance decreased in year 14 compared to year 5 in both treatments due to decomposition and lack of CWD input. While it is unsurprising that no new CWD was recruited in the clear cut, selected partial cut plots in year 14 recorded less tree mortality compared to year 5 (Moussaoui et al., 2020). Low retention levels in the year 14 partial cuts may have accounted for lower volumes of newly recruited CWD.
Our second hypothesis was partially supported, as epixylic species richness increased in both partial cut and clear-cut stands in year 14 compared to year 5. Habitat conditions improved in both partial cut and clear-cut stands due to high volumes of well decayed CWD and higher moisture conditions due to canopy closure from tree regeneration and growth. As has been found in other studies in different systems, bigger and well decayed CWD supported high epixylic richness as bigger CWD stays longer on the forest floor due to slow decomposition, giving more time for epixylic species to colonize the habitat (Söderström, 1988a;Andersson and Hytteborn, 1991). Bigger CWD also physically separates epixylic species from fast growing forest floor species thereby reducing CWD overgrowth and competition (Dynesius et al., 2010;Arseneault et al., 2012). CWD in advanced decay stages are generally decorticated with soft and spongy wood, which makes them better able to absorb water and help regulate forest floor moisture for epixylic bryophytes (Rambo, 2001;Haughian and Frego, 2017). Thus, high volumes of well decayed CWD recorded in both partial cut and clear-cut stands of year 14 partly explains the increase in epixylic richness over time.
Additionally, differences in decay rate associated with diameter size may explain the switch of the occurrence of most epixylic species from smaller woody debris on the forest floor in year 5 to CWD in year 14. Small harvest generated pieces were at a more advanced decay class compared to the harvest generated CWD in year 5. Consequently, most epixylic species were found on these pieces at that time (Arseneault et al., 2012). However, the small pieces disappeared from the forest floor quickly due to both rapid decomposition and overgrowth by forest floor species such as Pleurozium scherberi and Sphagnum spp. Considering potential epixylic bryophyte dispersal limitation (Söderström, 1987;Söderström and Jonsson, 1989) and their competitive disadvantage with forest floor species, smaller woody debris might only favor species with rapid growth and frequent reproduction. In year 14 the harvest generated CWD had advanced in decay stage and its longer residence time and bigger dispersal target allowed more species to colonize this habitat compared to year 5. Even though higher volumes of advanced CWD were recorded in both partial cut and clear-cut stands, the observed differences in epixylic community between the two harvest types highlights the importance of moisture conditions as reported by Söderström (1988a) and others.
Canopy openness had a negative effect on epixylic bryophytes. More open canopy increases forest floor temperatures and reduces moisture conditions (Fenton and Frego, 2005;Hylander, 2005). This dries up CWD making them an unsuitable substrate for these poikilohydric organisms (Arseneault et al., 2012;Fenton et al., 2013). This is evident in lower epixylic richness in clear cut stands in both Year 5 and Year 14. More open canopy in clear cut stands, favored the abundance of generalist and large-sized species such as feathermosses that are more tolerant to drought and competition (Caners et al., 2013;Boudreault et al., 2018). However, tree regeneration observed in clear cut stands created patches of shaded floor that improved moisture conditions and consequently created suitable substrates for epixylic species establishment resulting in a richer epixylic community on CWD in year 14 clear cut stands compared to those of year 5. Despite the increase in epixylic richness in clear cut stands over time, drought sensitive and old growth indicator species like Nowellia curvifolia and Blepharostoma trichophyllum (Boudreault et al., 2018) were absent in clear cut stands. This agrees with our third hypothesis and Söderström (1988b) assertion that drought sensitive epixylics are the most threatened species in managed landscapes.
Furthermore, the closed and open canopy generated by residual and regenerating trees in partial cut stands created heterogeneous microclimates on the forest floor. This heterogenous forest floor microclimate accounted for richer bryophyte community (including epixylic, liverworts and small-sized species), as different species have different moisture requirement due to their functional traits. The presence of richer specialized species in partial cut stands conformed to our second hypothesis. Additionally, the heterogenous microclimate conditions coupled with abundance of advanced decay CWD mimicked the complex habitat conditions observed in old growth forest. This partly explains old growth species assemblages including the presence of old growth indicator species (Nowellia curvifolia and Blepharostoma trichophyllum) in partial cut stands. This result also confirms our third hypothesis in relation to old growth species assemblages in partial cut stands. However, that was not the case in 2009 as species assemblages on logs fell outside the range of natural variability. This might be attributed to habitat degradation due to high tree mortality after harvest observed in these plots (Arseneault et al., 2012;Moussaoui et al., 2020). Whereas, in our case, stands in the selected plots experienced less tree mortality and there was a net growth in residual tree volume ten years after harvest (Moussaoui et al., 2020).
Continual availability of suitable substrates can also affect species richness as Siitonen (2001) describes for saproxylic species. As a consequence of the dynamic nature of CWD, epixylic bryophytes maintain a colonist life habit as they must colonize new suitable substrates as the old substrate becomes unsuitable because of competition and loss via decay and overgrowth. The lower volumes of early decay classes observed in partial cut and clear-cut stands put in doubt the future availability of substrate for epixylic bryophytes. In the case of clear-cut harvest, regenerating stands with no large trees will not provide a reliable source of CWD for epixylic bryophytes in the future and substrate availability will decrease in the near future due to the high proportion of well decayed CWD. Thus, the increase in epixylic richness observed in clear cut stands in our study might not last. Unlike clear cut, residual trees in partial cut stands provide a potential source of CWD supply.
The changes in partial cut plots between years 5 and 14 may have contributed to the observed differences in epixylic community. High tree mortality in Year 5 partial cut plots resulted in high volumes of CWD, but the opened canopy resulting from tree mortality exposed CWD to the direct effect of sun rays. Consequently, CWD dried up and became less suitable for epixylic species colonization (Andersson and Hytteborn, 1991;Fenton et al., 2003). Thus, habitat conditions in some of Year 5's partial cut plots were less favorable for epixylic bryophytes compared to selected plots in Year 14 which experienced less tree mortality.

CONCLUSION
As an emerging alternative to clear cut harvest, understanding both the immediate and long-term benefits of partial cut harvest is crucial for its adoption. Our study highlights the ability of partial cut harvests to reduced impacts on forest biodiversity by continuing to provide favorable microhabitat conditions, which supported epixylic bryophytes. Also, even low harvest retention (<25%) levels, facilitated faster post-harvest recovery of bryophytes, as compared to clear cut. This was also consistent with reports on the ability of partial cut harvests to offer an effective strategy for biodiversity conservation compared to traditional clear cuts in managed forest landscapes (Bescond et al., 2011;Fenton et al., 2013;Bartels et al., 2018). However, lower volumes of newly recruited coarse woody debris raise concerns about the deadwood delivery potential of partial cut harvests. Deadwood input should therefore be considered in implementation strategies to ensure continual persistence of epixylic bryophytes and deadwood living organisms in general.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding author/s.

AUTHOR CONTRIBUTIONS
JO-N, AL, and NF conceived and designed the study. JO-N collected and analyzed the data. All authors contributed to manuscript writing and gave final approval for publication.

ACKNOWLEDGMENTS
We thank Julie Arseneault for her technical assistance in the field and during bryophyte identification. Also, we thank Manuel Desrochers, Samuel Roy, and Sylvain Lemieux for their field assistance and Danielle Charron and Julien Raynald for field logistics support.