Balsam Fir and American Beech Influence Soil Respiration Rates in Opposite Directions in a Sugar Maple Forest Near Its Northern Range Limit

Conifers and deciduous trees greatly differ in regard to their phylogenetics and physiology as well as their influence on soil microclimate and chemical properties. Soil respiration (Rs) in forests can therefore differ depending on tree species composition, and assessments of the variation in Rs in various forest types will lead to a more thorough understanding of the carbon cycle and more robust long-term simulations of soil carbon. We measured Rs in 2019 and 2020 in stands of various species composition in a sugar maple forest near the northern range limit of temperate deciduous forests in Quebec, Canada. Seasonal variations in soil temperature had the largest influence on Rs, but conditions created by the stands also exerted a significant effect. Relative to the typical sugar maple-yellow birch forest (hardwoods), Rs in stands with >20% of basal area from balsam fir (mixedwoods) was increased by 21%. Whilst, when American beech contributed >20% of litterfall mass (hardwood-beech stands), Rs was decreased by 11 and 36% relative to hardwoods and mixedwoods, respectively. As a whole, Rs was significantly higher in mixedwoods than in other forest types, and Rs was significantly higher in hardwoods than in hardwood-beech stands. Sugar maple and American beech at the study site are near their northern range limit, whereas balsam fir is near its southern limit. Rs in mixedwoods was therefore higher than in hardwoods and hardwood-beech stands due to high root activity in the presence of fir, despite colder and drier soils. We estimated that root respiration in mixedwoods was more than threefold that in hardwoods and hardwood-beech stands. The lower Rs in hardwood-beech stands compared to hardwoods points to the lower soil temperature as well as the poor quality of beech litter (low decomposability) as indicated by a generally lower heterotrophic respiration. Other than soil temperature, regression models identified mixedwoods, soil water potential and Mg2+ activity in the soil solution as important predictor variables of Rs with about 90% of its variation explained. Our study shows the benefits of combining forest-specific properties to climatic data for more robust predictions of Rs.

Conifers and deciduous trees greatly differ in regard to their phylogenetics and physiology as well as their influence on soil microclimate and chemical properties. Soil respiration (R s ) in forests can therefore differ depending on tree species composition, and assessments of the variation in R s in various forest types will lead to a more thorough understanding of the carbon cycle and more robust long-term simulations of soil carbon. We measured R s in 2019 and 2020 in stands of various species composition in a sugar maple forest near the northern range limit of temperate deciduous forests in Quebec, Canada. Seasonal variations in soil temperature had the largest influence on R s , but conditions created by the stands also exerted a significant effect. Relative to the typical sugar maple-yellow birch forest (hardwoods), R s in stands with >20% of basal area from balsam fir (mixedwoods) was increased by 21%. Whilst, when American beech contributed >20% of litterfall mass (hardwood-beech stands), R s was decreased by 11 and 36% relative to hardwoods and mixedwoods, respectively. As a whole, R s was significantly higher in mixedwoods than in other forest types, and R s was significantly higher in hardwoods than in hardwood-beech stands. Sugar maple and American beech at the study site are near their northern range limit, whereas balsam fir is near its southern limit. R s in mixedwoods was therefore higher than in hardwoods and hardwood-beech stands due to high root activity in the presence of fir, despite colder and drier soils. We estimated that root respiration in mixedwoods was more than threefold that in hardwoods and hardwood-beech stands. The lower R s in hardwood-beech stands compared to hardwoods points to the lower soil temperature as well as the poor quality of beech litter (low decomposability) as indicated by a generally lower heterotrophic respiration. Other than soil temperature, regression models identified mixedwoods, soil water potential and Mg 2+ activity in the soil solution as important predictor variables of R s with about 90% of its variation explained. Our study shows the benefits of combining forest-specific properties to climatic data for more robust predictions of R s .

INTRODUCTION
Soil respiration (R s ) is the second largest CO 2 flux in forests after plant respiration and in turn, it has a high potential to modify atmospheric CO 2 levels (Schlesinger and Andrews, 2000;Raich et al., 2002). It is constituted of heterotrophic (R h ) and autotrophic (R a ) respiration. The R h component comes from the respiration of the soil faunal and microbial communities that physically breakdown and biochemically decompose dead organic matter, whereas the R a component comes from root and rhizosphere respiration. Soil temperature is the main abiotic control of R s in forests (Davidson and Janssens, 2006;Subke and Bahn, 2010). The temperature sensitivity of R s is often expressed as Q 10 − it is the proportional change in respiration with a 10 • C increase in soil temperature (Curiel Yuste et al., 2004). Soil moisture also interacts with soil temperature to influence R s in forests (Cisneros-Dozal et al., 2006;Moyano et al., 2012). Soil temperature and moisture may vary seasonally and across years, thereby modifying the temperature sensitivity of R s in forest ecosystems (Davidson et al., 1998). Other factors also influence the seasonal patterns of R s in forests such as variation in root activity during the growing season and fresh litter additions in the fall Lavigne et al., 2004;Prévost-Bourré et al., 2010).
Changes in forest structure and species composition that lead to a modification in soil conditions are thus expected to influence R s . Soil respiration tends to increase with stand age due to root growth and respiration (Subke et al., 2006), but this trend can be absent in forests where crown closure reduces soil temperature (e.g., Gough et al., 2005). Forest species composition may also influence R s by affecting: (1) soil chemical properties (e.g., pH, N, and P availability), (2) root biomass, respiration and phenology, (3) litter inputs, quality and decomposability, and (4) microbial communities and activity (Raich and Tufekcioglu, 2000). Side-by-side comparisons of R s between deciduous and coniferous stands generally suggest higher fluxes under deciduous trees (Tewary et al., 1982;Weber, 1985Weber, , 1990Raich and Tufekcioglu, 2000;Curiel Yuste et al., 2004;Fahey et al., 2005), but lower R s under deciduous (Lee et al., 2010) or no difference in R s between forest types (Reiners, 1968;Raich and Potter, 1995;Davidson et al., 1998;Borken et al., 2002;Fernández-Alonso et al., 2018) have also been reported.
The generally lower quality of conifer litter and R h relative to deciduous tree species was a common explanation for the lower R s measured in coniferous stands, whereas differences in phenology and physiology between tree species explained no change or higher R s in coniferous stands.
Worldwide, temperate deciduous forests were estimated to hold 7-8% of global terrestrial C stocks (Pregitzer and Euskirchen, 2004). Cool temperate deciduous forests in eastern Canada offer high C sequestration potential because they are among the most productive ecosystems in the country due to their relatively long growing season and sufficient precipitation falling during the growing season (Kurz and Apps, 1999;Stinson et al., 2011). Also, there is little fire disturbance or forest harvesting (mostly uneven age management) that could release large amounts of C back to the atmosphere (Saucier et al., 2009a). Fahey et al. (2005) provided R s estimates for various deciduous and coniferous stands in the cool temperate Hubbard Brook Experimental Forest, New Hampshire. However, there is no study that investigated the variations in R s near the northern limit of temperate deciduous forests in northeastern North America. These forests are associated with common tree species changes due to the presence of deciduous species to the south and coniferous (boreal) species to the north. Furthermore, northeastern North American deciduous forests are characterized by replacement patterns between sugar maple (Acer saccharum) and American beech (Fagus grandifolia). This replacement pattern is also observed at their northern limit (Collin et al., 2017). Light, water and nutrient availability as well as reproduction rates govern the replacement patterns between the two species. These factors can be modified due to natural disturbances such as ice storms, disease and herbivory (Arii and Lechowicz, 2002;Bohn and Nyland, 2003;Nyland et al., 2006), forest management practices (Nolet et al., 2008;Bannon et al., 2015) and atmospheric acid deposition (Bailey et al., 2004;Long et al., 2009). Water limitations for tree growth are suspected to increase in frequency and intensity under climate change in eastern North America (Gustafson and Sturtevant, 2013), whereas sugar maple is likely more sensitive than American beech to water stress (Nolet and Kneeshaw, 2018). Overall, the recent increase in the dominance of beech in the understory of maple stands suggests that beech stands will occupy forest land in the future (Hane, 2003;Gravel et al., 2010). Similar to conifer litter, the high recalcitrance of American beech litter (Melillo et al., 1982;Côté and Fyles, 1994) has the potential of reducing R s by slowing down litter turnover and thus decreasing R h (Bowden et al., 1993;Cisneros-Dozal et al., 2007).
Climate change is causing physiological constraints that have the potential of impacting the survival, growth and distribution of plant species globally (Rosenzweig et al., 2008;Chen et al., 2011), including northeastern North American tree species (Beckage et al., 2008;Charney et al., 2016;Sittaro et al., 2017;Iverson et al., 2019;Boivert-Marsh and de Blois, 2021). In addition to understanding how tree growth and thus forest C uptake will be affected by climate change, analyses of how forest species compositions affect R s are needed for more robust predictions of temporal and spatial changes in ecosystem C pools and atmospheric CO 2 (Raich and Tufekcioglu, 2000;Buchholz et al., 2014). As sugar maple, American beech and conifers differ in regard to their physiology, biogeochemistry (e.g., litter quality/decomposition) and phylogenetics (angiosperm vs. gymnosperm), R s is expected to vary depending on the proportion of these tree species within a forest. The objective of this study was to assess R s across a range of plots that captured variations in the abundance of conifers and American beech within a forest dominated by sugar maple near the northern range limit of temperate deciduous forests in eastern North America. Mostly due to their lower litter qualities compared to maple litter, it was hypothesized that an increase in the abundance of conifers and beech would decrease R s in this maple forest.

Study Site
The study was conducted at the Station de Biologie des Laurentides (SBL) of Université de Montréal in St. Hippolyte, Quebec (Figure 1). The SBL is found within the northern limit of the maple-yellow birch (Betula alleghaniensis) bioclimatic domain of the lower Laurentians. The maple-yellow birch domain is the northernmost deciduous forest domain in Quebec (Saucier et al., 2009b). A mosaic of tree species is found at the site. It is composed mostly of sugar maple, red maple (Acer rubrum), American beech, yellow birch (Betula alleghaniensis), white birch (Betula payrifera), largetooth aspen (Populus grandidentata), eastern white cedar (Thuja occidentalis), white pine (Pinus strobus), and red spruce (Picea rubens) (Savage, 2001). Because SBL is only 35 km south of La Macaza, i.e., the doorway to temperate mixedwoods at that specific longitude in Quebec (Figure 1), stands dominated with tree species that are typically boreal [e.g., balsam fir (Abies balsamea) and white spruce (Picea glauca)] are not uncommon at the site. Also, stands exhibiting sugar maple regeneration failures and expansion of beech are frequent at SBL (Collin et al., 2017). The mean annual temperature at SBL simulated with the BioSIM model (Régnière and Bolstad, 1994) between 2003 and 2013 was 4.9 • C, mean degree-days were 2845, mean days without frost were 153 and mean precipitation was 1,270 mm, with about 30% falling as snow. Soils were developed from glacial till made mostly of anorthosite and felsic rocks of the Precambrian Shield (Bélanger et al., 2012). They are classified as Orthic Humo-Ferric and Ferro-Humic Podzols with a sandy loam texture and the forest floor is characterized by a moder humus form (Soil Classification Working Group, 1998).

Experimental Design and Stand/Plot Characterization
In 2018, eight stands with four 3 × 3 m plots within each stand (total of 32 plots) were selected in three zones distributed within a 18 ha area. Stands were selected to capture as much variation as possible in regard to tree species composition. Stands were a minimum of 0.5 ha. The four plots within each stand were delineated in the center of the stand (within a maximum of 15 m between the most distant plots) under a very similar tree species composition. A series of ecological variables were collected at the stand level between 2017 and 2020, allowing for a detailed classification of tree species composition. These included (1) basal area by species based on trees with a stem diameter at breast height above 9 cm (Table 1) and (2) litterfall mass by species  -beech is hardwood-beech, LAI is leaf area index, Leaves is all leaves, Needles is mostly Abies balsamea (balsam fir), AS is Acer saccharum (sugar maple), Betula is Betula alleghaniensis and Betula papyrifera (yellow and white birch, respectively), FG is Fagus grandifolia (American beech), AR is Acer rubrum (red maple), PG is Populus grandidentata (large-tooth aspen), and AP is Acer pensylvanicum (striped maple). The sum of Leaves and Needles or the sum of Needles and AS, Betula spp., FG, AR, PG, and AP equals 100%.
( Table 2). At the center of each stand, litterfall was collected in a plastic bin (0.9 cm × 0.7 cm × 0.5 cm) that was perforated at the bottom and filled with silica sand (depth of about 5 cm) to drain the bin and thus prevent the trapped litter to immerse in water between samplings. Litter was collected each fall and dried in an oven at 65 • C for 48 h before being weighted by species. Leaf area index was also measured in the center of each plot on a sunny morning of mid-August 2020 when the canopy was fully developed using a CI-110 Plant Canopy Imager (CID Bio-Science, Camas, WA). Stands were classified as mixedwoods when basal area of balsam fir contributed to more than 20% of total basal area (i.e., stands 1, 4, and 7, Table 1). However, basal area by species may not be fully representative of deciduous litterfall type in the plots. For example, the high American beech sapling density in some stands was not captured in the basal area measurements, despite that beech litterfall is high (i.e., stands 3 and 8, Table 2). Because R h is a significant component of R s in temperate deciduous forests [25-35% of R s according to Bowden et al. (1993) and Cisneros-Dozal et al. (2007)], we used litterfall mass by species to separate between hardwoods and hardwood-beech stands. In this respect, hardwood-beech stands were identified as having at least 20% of their total litterfall mass from beech (stands 3, 5, and 8), whereas hardwoods had less than 20% (stands 2 and 6, Table 2). In the case of stand 2, there was no presence of beech.
Soil temperature and water availability at a 10 cm depth were measured every 15 min in each plot from two temperature sensors (Spectrum Technologies, United States) and two water potential sensors (Irrometer 200SS-5, Watermark, United States), all connected to the same WatchDog 1650 Micro Station data logger (Spectrum Technologies). Plant Root Simulator (PRS) probes (Western Ag Innovations, Canada) were used to assess ionic activity (i.e., NO 3 − , NH 4 + , H 2 PO 4 − , Ca 2+ , Mg 2+ , K + , Al 3+ , Fe 3+ , and Mn 2+ ) in the soil. Four pairs of cation and anion probes were carefully inserted vertically into the forest floor at random locations within each plot. They were installed in early June of 2019 and 2020 and collected 6 weeks later. The PRS probes represent a dynamic measurement of ions flowing through the soil over time compared to conventional soil extraction methods that provide a measurement of soil nutrient availability at a particular point in time. The probes are now frequently used in forest ecology research (Hangs et al., 2004;Bilodeau-Gauthier et al., 2013). Once extracted from the soil, the probes were cleaned with deionized H 2 O and stored in the fridge in zipseal bags until analysis. Probes were eluted for 1 h with 0.5 M HCl. NH 4 -N and NO 3 -N were determined colorimetrically by continuous flow analysis (Autoanalyser III, Bran and Luebbe, United States), whereas other ions were determined by inductively coupled plasma atomic emission spectroscopy (Optima 3000-DV, PerkinElmer, United States).

Gas Sampling, Laboratory Analysis, and Flux Calculations
Four circular bases were installed in each plot in July 2018. The bases were set out in a square shape pattern with a distance of 1.5 m between them. They were fabricated from high-density PVC 15 cm inside diameter pipes (i.d., surface area of 177 cm 2 ) with smooth 6 mm thick side walls. The bases were cut to a length of 12 cm and were inserted at a soil depth of 6 cm. The flux chambers were fabricated from the female connection end of the pipes (i.d. of 15.5 cm, surface area of 189 cm 2 ) with a double seal lockedin gasket, thus securing a tight seal with the base. The chamber was closed with a 4 mm thick Teflon sheet and the exterior was wrapped with reflective insolation to avoid overheating. The combination of the base and the length of the chamber results in a headspace volume of 4.56 L. The risk of pressurization of the flux chamber was eliminated by (1) venting with a Tygon tube (4 mm i.d.) that passed through the top of the chamber and (2) avoiding to sample under windy conditions. Gas sampling was done from an injection site (Bung Interlink IV, Baxter, United States). Debris were removed only when they obstructed the placement of the chamber on the base.
In 2019, gas sampling started in June and finished just before soil freezing in November 2019. It was done in all the plots within each stand and at five different dates (n = 160). Despite the large number of plots (4 plots per stand, 8 stands, and 32 plots in total), sampling was done within a day (∼6 h). In 2020, we opted for a higher sampling frequency but did not sample all stands or all four plots in each stand each time. We sampled two plots in all stands on June 9, June 22, July 9, July 22, August 26, and September 23, all plots in stands 1, 2, 3, and 4 (i.e., block 1) on June 15, June 30, July 21, and August 11, and one plot per stand on October 22. A rotation of plots was done for a better assessment of spatial variability within each stand. In total, 160 gas samples were collected at 11 different dates in 2020, but 4 samples were discarded due to sampling errors (n = 156).
Gas sampling was done as soon as the flux chamber was deposited on the base (t 0 ). To do so, 5 ml of gas sample was withdrawn from the chamber using a 50 ml polypropylene syringe equipped with a 25-gauge 5/8-inch needle. The gas samples collected from each chamber (4×) were injected in the same pre-evacuated (ca. 0.005 atm.) 12 ml Exertainer R vial (LabCo, United Kingdom). The same procedure was repeated at 4 (t 4 ), 8 (t 8 ), 16 (t 16 ), and 24 (t 24 ) min. The gas samples were stored under a positive pressure of approximately 1.7atm, i.e., respectively 20 ml of headspace gas into a 12 ml vial, to minimize any gaseous exchange with atmospheric air. The sampling scheme (i.e., pooling gas samples) followed the recommendation of Arias-Navarro et al. (2013) as a means to better capture spatial variability of soil gas efflux within the plots.
To provide an estimate of the contribution of the autotrophic (R a ) and heterotrophic (R h ) components to R s , we used a root exclusion approach similar to Kelting et al. (1998). In mid-May 2020, we created four 1 × 1 m root-free plots by severing all the roots to a depth of 40 cm. This depth was considered to capture most of the roots because Lajeunesse (1990) measured that 70% of root mass at SBL is found within the first 25 cm of soil, i.e., the foret floor and upper podzolic B horizon. We had one plot in stands 1, 2, 3, and 4, being, respectively, mixedwood, hardwood, hardwood-beech and mixedwood. All four plots were delineated in proximity of the other sampling plots. There was barely any vegetation in those plots and any new growth was removed by clipping. We used the same sampling procedure as described above. No barrier was used to further limit root development because CO 2 efflux measurements were planned to be performed in 2020 only. Sampling of these plots started on June 15, exactly 1 month after trenching, and continued on June 30, July 21, August 11. This sampling was done at the same time as the sampling scheme described above. For each forest type and sampling date, root respiration (R root ) was calculated as the difference between R s (i.e., mean of the standard 3 × 3 m plots) and soil respiration in the root-free plots (1 × 1 m plot), whereas R h was assumed as the net flux from the trenched plots. The percent contributions of R root and R h to R s were also calculated. Soils can become wetter in the trenched plots due to the absence of water uptake by roots and this can affect soil respiration patterns when they exhibit extremes in moisture levels (Hanson et al., 2000). In addition, soil respiration from the trenched plots could include an added flux of CO 2 because of a possible increase in the decomposition of dead roots [rhizosphere (or rhizomicrobial) respiration], especially a few months after trenching (thus the use of R root and not R a ). This means that R root could be an underestimation of root respiration, whereas R h could be overestimated. Bowden et al. (1993) provided a strong argument that the contribution of decomposition of dead roots in trenched plots is small, at least in the short term. However, because we do not have estimates of root decomposition in the trenched plots, we need to interpret the data with care. The estimates were therefore used as a means to exhibit the main contrasts in R root and R h patterns among the three forest types studied.
Carbon dioxide and CH 4 analysis of gas samples was carried out within 48 h of completion of each sampling campaign. Cavity-ring down (CRD) spectroscopy (G2201-i Isotopic Analyzer, Picarro, United States) and ultra-zero air as the carrier gas were used for analysis. Before analysis, water vapor from gas samples was selectively removed using a monotube Nafion TM gas dryer (MD-050-12-2, PermaPure, United States).
Raw CO 2 fluxes were estimated by fitting t 0 , t 4 (exp. 2 only), t 8 , t 16 , and t 24 data to a linear model using the HMR package in R (Pedersen et al., 2010). We used a linear model for all flux samples for two main reasons: (1) it was the best fit for over 95% of samples, and (2) it provides more consistent flux results for short enclosure times and in cases of a concave flux curve, i.e., the most common non-linear pattern for the remaining 5% of samples (Kutzbach et al., 2007;Görres et al., 2014;Kandel et al., 2016). Selecting a single linear model thus avoided calculating differences in CO 2 fluxes due to differences in models and not CO 2 fluxes per se. Air temperature and pressure were used to adjust raw gas fluxes (Rochette and Bertrand, 2007). The well drained soils resulted in small CH 4 sinks and therefore, CH 4 fluxes were not investigated any further.

Data Analysis
To confirm the influence of soil temperature on R s , we first fitted single exponential models of R s with soil temperature using the equation R s = a × e (b ×soiltemp) , where a and b are parameters to be estimated, e is the base of the natural logarithm (2.71828) and soil temp is measured soil temperature in • C. Models were fitted for each forest type and for all forest types combined (general model) for both years combined. We calculated R 2 from the linear regression between measured vs. predicted R s . The Q 10 value for each model was obtained with the equation Q 10 = e (b × 10) . The residuals of the general model were then modeled using a forward stepwise regression with the following independent variables: soil water potential, total basal area, LAI, and NO 3 − , NH 4 + , H 2 PO 4 − , Ca 2+ , Mg 2+ , K + , Al 3+ , Fe 3+ , and Mn 2+ activity in the soil solution as measured by PRS probes. Forest types were also tested after they were converted into dummy variables, i.e., mixedwoods (0) vs. other forest types (1), hardwood-beech stands (0) vs. other forest types (1), and hardwoods (0) vs. other forest types (1). Only the first three selected variables were kept in the model. This procedure may leave out some variables that are not necessarily unimportant. However, our purpose was to find a small set of logical predictor variables that did an adequate job of prediction while being relevant to biological theory (Sokal and Rohlf, 2012). All predictor variables individually improved the fit between the observed and predicted values at P < 0.01. The marginal sum of squares (type III) was also used to measure the predictive information contained in the predictor variables individually while considering the other predictor variables in the model. The variance inflation factor (VIF) was used to detect any co-linearity between the independent variables in the model. No VIF exceeded a value of 2.8. Transformation of the raw data was not necessary to assure equal variance and normality in the distribution of the residuals. Predicted R s was plotted against measured R s and the R 2 of the linear regression was calculated as an indication of the capacity of this two stage modeling approach to explain the variation in R s .
To assess the effect of forest types on R s , we used two linear mixed effect models with forest type as the fixed factor in both models. In the first model, sampling date was used as the random factor, whereas soil temperature was used as the random factor in the second model. We proceeded this way to consider that soil temperature differences may be produced by the forest types and their canopies and, in turn, this could lead to differences in R s . A square root transformation was used in both models to normalize the residuals, whereas predicted values were computed for subsequent display. To assess the effects of forest types on soil temperature and water potential, we first computed daily averages from May 15 to October 15 in 2019, and from May 15 to August 31 in 2020, i.e., the data that we could secure while avoiding logistical problems associated with snow cover in both years. Technical problems with dataloggers reduced the continuity of the data in 2020 and thus, we did not use time-series data for September and October for further analysis. With these data, we used linear mixed effect models with forest type as the fixed factor and sampling date as the random factor. A logarithmic transformation was used to warrant normality of the residuals. Finally, one-way ANOVA was used to detect differences in R root and R h between forest types as well as their respective contributions (in %) to R s . Block identity was not tested in any of the linear mixed effect models and ANOVAs because it was not possible to obtain all forest types within each zone and thus, zones should not be considered as a blocking structure per se. Significant differences between forest types detected with these tests were further depicted with post hoc multiple comparisons using the Tukey's honest significant difference (HSD) test. Linear mixed models and ANOVAs were, respectively, conducted with the R statistical software (version 4.0.4) and the "Analysis" module in SigmaPlot 12.0.

RESULTS
Soil respiration rates varied from as low as 21 mg CO 2 m −2 h −1 in hardwood-beech stands on November 8, 2019, to as high as 726 mg CO 2 m −2 h −1 in mixedwoods in July 9, 2020 (Figure 2). The large variation in R s was associated to seasonality, increasing from June to July-August and decreasing again to the lowest values just before snow in late October of 2020 or early November of 2019 (Supplementary Figure 1). The seasonal effect was more apparent in 2020 because of a higher sampling frequency. Relatively strong exponential relationships were found between soil temperature and R s for each forest type and all forest types combined (general model) (Figure 2). Model parameter estimates, Q 10 and R 2 are shown in Table 3. The modeled curve for mixedwoods was consistently above all other models. Conversely, the modeled curve for hardwoodbeech stands was below all other models, although the regression lines between hardwoods and hardwood-beech stands crossed each other at a soil temperature of 15.5 • C. Soil temperature alone explained 59% of the variance in R s rates (see general models in Table 3). The predictive ability was higher for hardwood-beech stands and mixedwoods (64-65%) and lower for hardwoods (55%). Computed Q 10 values were 2.26 for hardwoods, 2.44 for mixedwoods, 2.58 for all forest types, and 2.97 for hardwoodbeech stands ( Table 3).
The forward stepwise regression model of the residuals of the general exponential model (i.e., soil temperature vs. R s ) explained an additional 19% of the variation in R s ( Table 4). The dummy variable representing mixedwoods (0) vs. other forest types (1) was selected as the first predictor variable in all models and explained between 11% of the variation in R s . Soil water potential was selected as the second predictor variable (5.3%), whereas Mg 2+ activity in the soil solution was selected as the last predictor variable (3.0%). In the end, the linear relationship between measured and predicted R s over the two study years suggest that the combination of both models (i.e., exponential model followed by multiple linear model) captured 89% of the total variation in R s . A graphical illustration between measured and predicted values shows that this approach underestimated lower R s values and overestimated higher R s values (Figure 3). FIGURE 2 | Exponential relationships between soil temperature and soil respiration rates (R s ) for each forest type individually and for all forest types combined (general model). Parameter estimates and Q 10 of equation R s = a × e (b ×soiltemperature) are presented in Table 3.
Over the two study years, R s in mixedwoods was 21 and 36% higher than hardwoods and hardwood-beech stands, respectively (Figure 4). The difference in R s between hardwood and hardwood-beech stands was 11%. Mixed model analysis with sampling date as a random variable suggest a significant difference in R s between all forest types, whereas mixed model analysis using soil temperature as a random variable suggest a significantly higher R s in mixedwoods than in hardwoods and hardwood-beech stands, but not between hardwoods and hardwood-beech stands ( Table 3).
Soil temperature exhibited a seasonal variation similar to R s , increasing from June to July-August and decreasing again to < 10 • C in October and < 4 • C in November (Supplementary Figure 2). Soils reached higher temperatures in 2020 than in 2019. In 2019, soil water potential reached a maximum in early August (Supplementary Figure 2). In 2020, a series of early heatwaves and drought events explain the unusually high soil TABLE 3 | Parameter estimates and Q 10 values of exponential models (Figure 3) between soil temperature respiration rates for each forest type individually and for all forest types combined (general model). temperature and very high water potential recorded in June.
Overall, soil water potential in 2020 reached higher values than in 2019. Linear mixed effect models detected significantly lower soil temperature in mixedwoods, followed by hardwoodbeech stands, and then hardwoods, whereas soil water potential was significantly lower in hardwood-beech stands, followed by hardwoods and then mixedwoods (Supplementary Figure 2). Root-free plots were used to estimate R root and R h in the three forest types and their contributions (%) to R s . Root respiration in mixedwoods was approximately threefold higher than in hardwoods and hardwood-beech stands, whereas R root in hardwood-beech stands was 22% higher than in hardwoods ( Table 5). The R h component in mixedwoods and hardwoodbeech stands generated a smaller CO 2 flux than hardwoods by about 26 and 11%, respectively. Due to the low degrees of freedom of the ANOVAs, however, the only significant difference was found between R root of mixedwoods and R root of hardwood and hardwood-beech stands. Root respiration contributed about 20% of R s in hardwoods, 25% in hardwood-beech stands, and 50% in mixedwoods ( Table 5). The remaining CO 2 flux thus originated from the heterotrophic component. The contributions of R root and R h to R s were, respectively, significantly lower and higher in hardwoods and hardwood-beech stands than in mixedwoods.

DISCUSSION
The seasonal variation in soil temperature in this sugar maple forest near the northern distribution limit of temperature deciduous forests was conducive to relatively strong exponential Adj. R 2 is adjusted R 2 and P is the P-value, SEE is the standard error of estimate, and the last three columns are the predictor variables (i.e., Mixedwoods is a dummy variable [i.e., mixedwood (0) vs. other forest types (1)], soil water potential, and magnesium activity in the soil solution. Full lines report performance and parameters of the models, whereas delta R 2 , partial P and SSmarg [i.e., marginal sum of squares (type III)] provide an assessment of the predicting role of each variable to the model.
FIGURE 3 | Measured soil respiration rates (R s ) vs. predicted R s after the two stage modeling approach, i.e., exponential of the raw R s data and soil temperature ( Figure 2 and Table 3) followed by a multiple linear regression of the residuals ( Table 4). The solid line is the linear regression between measured vs. predicted R s , whereas the dashed line is the 1:1 line.
relationships between soil temperature and R s for all forest types. This result was expected because, at all spatial scales, soil temperature is the main abiotic control of R s and its associated processes (Kutsch et al., 2009) and the key variable for predicting R s in the context of climate change (Subke and Bahn, 2010). However, further statistical analyses of our data indicated that, over the 2 years of measurements at SBL, tree species composition has had a substantial influence on R s as well. First, forest stands where balsam fir contributed at least 20% of basal area (i.e., mixedwoods) significantly increased R s , whereas forests where American beech contributed at least 20% of litterfall (i.e., hardwood-beech stands) significantly decreased R s . These results contradict our initial hypothesis that conifers would lower R s in this typically maple dominated forest, but they validate our hypothesis that beech would lower R s .

Influence of Mixedwoods on R s
Most studies that investigated the influence of forest type on R s were conducted at spatially broad experimental scales and thus emphasized the direct effect of temperature and forest productivity across landscapes instead of forest type per se. Because R s is largely dependent on the translocation of photosynthates produced aboveground to the roots (R a component, Högberg et al., 2001), the expected trend is for R s to increase from colder to warmer forests. As such, meta-analyses show an increase in R s from boreal to temperate to tropical forests (Raich and Schlesinger, 1992;Subke et al., 2006). However, within FIGURE 4 | Mean soil respiration rates (R s ) predicted from the linear mixed effect models for each forest type using sampling date and soil temperature as the random variables. Predicted values were back-transformed as a square root transformation was used in both models to assure normality of the residuals. Standard errors are shown. Different letters indicate a significant difference at P < 0.01 between forest types using the Tukey's honest significant difference test. Contributions of both components to total soil respiration were also calculated. Estimates are means and standard errors (in parentheses or under S.E.) of the four sampling dates (June 15, June 30, July 21, and August 11). N.B. The autotrophic component does not account for rhizosphere respiration and thus, it is described as root respiration (R root ). This means that autotrophic respiration is underestimated and R h is overestimated. Different letters indicate a significant difference at P < 0.01 between forest types using one-way ANOVA.
a constrained geographical location, an increasing abundance in conifers can also modify the soil system to conditions (e.g., lower temperatures, pH and nutrient availability) that are less suitable for microbial degradation of soil organic matter and leaf litter (Binkley and Giardina, 1998;Binkley and Fisher, 2012;Prescott and Grayston, 2013;Joly et al., 2017). Such conditions would be expected to lower R h , i.e., a significant component of R s (25-35%) in temperate deciduous forests (Bowden et al., 1993;Cisneros-Dozal et al., 2007). Therefore, in a side-by-side study of deciduous and coniferous stands with similar photosynthetic rates, it would seem reasonable to find higher R s under deciduous stands because of the greater organic matter turnover (Raich and Tufekcioglu, 2000;Curiel Yuste et al., 2004;Fahey et al., 2005). At SBL, Collin et al. (2017Collin et al. ( , 2018 observed more acidic and nutrientpoor soils under mixedwoods and conifer dominated forests than hardwoods and hardwood-beech stands. Furthermore, in a study of four sites in the sugar maple-basswood (Tilia Americana) and sugar maple-yellow birch domains of Quebec, including SBL, Bélanger et al. (2019) found that decomposition rates of sugar maple leaf litter were slower under pure coniferous stands than under hardwood and hardwood-beech stands, and that decomposition rates under mixedwoods were intermediate. We thus hypothesized that the slower leaf litter decomposition in mixedwoods at SBL would lead to lower R s than hardwoods, but consistently higher R s under mixedwoods compared to the other forest types contradicts this hypothesis.
In a similar cool temperate climate in Korea, Lee et al. (2010) measured lower R s under deciduous forests than coniferous forests. They explained that there were greater photosynthetic constraints on leaves of deciduous trees. In turn, this led to limited photosynthate transport and lower root activity and R a . We do not have reliable dendrochronological data to infer differences in productivity among tree species and forest types at SBL. However, biomass data provided by Maliondo et al. (1990) suggest that balsam fir and white spruce stands produce equivalent or more biomass annually than sugar and red maple stands in northern temperate deciduous forests in eastern Canada. In this respect, mixedwoods at SBL, with balsam fir as the main coniferous species, are expected to perform very well. Conversely, growing conditions are sub-optimal for maple spp. and beech at SBL as they develop near their northern range limit and are supported by acidic and N-poor soils (Collin et al., 2017(Collin et al., , 2018. We thus argue that R s in mixedwood plots were balanced by a greater contribution of R a relative to hardwoods and hardwood-beech stands. This is well supported by our partitioning analysis of R root and R h . Despite colder and drier soils under mixedwoods, R root was substantially and significantly higher (about threefold) in mixedwoods than in the other forest types studied. We also estimated that R root in mixedwoods was about 50% of R s compared to only 20 and 25% in hardwoods and hardwood-beech stands, respectively. Raich and Schlesinger (1992) suggested that correlation between ecosystem productivity and R s was largely due to a regional/continental effect of climate, whereas Valentini et al. (2000) argued that ecosystem productivity can overshadow the influence of climatic variables on R s at large spatial scales. Similarly, Reichstein et al. (2003) were able to partition the contributions of climate and vegetation types on R s , but this was over an array of sites in Europe and America that encompassed a large gradient in soil microclimate and vegetation productivity. In our study, a dummy variable that characterized the mixedwood plots (0) (vs. hardwood-beech and hardwood plots 1) was the first predictor variable in a multiple linear regression of the residuals produced from the exponential model between R s and soil temperature. We were thus able to separate the large temperature effect on R s from the vegetation effect reflecting differences in productivity between mixedwoods and deciduous forests at a very small spatial scale. Raich and Potter (1995) and Fernández-Alonso et al. (2018) reported no difference in R s between deciduous and coniferous forests. In a Mediterranean ecotone forest of central Spain, Fernández-Alonso et al. (2018) partitioned the contribution of R a and R h to R s and found that deciduous stands had a lower R h and a greater R a than coniferous forests, but R s rates were similar between the two forests. This is an opposite scenario to the one at SBL where fir, not maple spp. or beech, produced higher R s due to higher R a . As discussed above, this finding is reasonable given that growing conditions at SBL are more suitable for fir than sugar maple and beech. It thus appear that patterns in R s between deciduous and coniferous forests under similar growth conditions cannot be easily generalized and that further sideby-side studies are needed to fully elucidate the factors, notably soils and climate, that govern differences in R s between trees with diverging phylogenetics and physiology.

Influence of Hardwood-Beech Stands on R s
Based on data extracted from Gosz et al. (1973); Melillo et al. (1982), Côté and Fyles (1994), and Moore et al. (1999), rates of American beech leaf litter decomposition are estimated to fall toward the lower end of the species present at the study sites: birch spp. ≥ red maple, aspen spp. ≥ sugar maple, spruce spp. > white pine > balsam fir, American beech > eastern white cedar. Our results are consistent with this sequence because the R h component is generally lower when leaf litter at SBL is enriched with beech and fir. However, due to the small sample and variability in the data, differences were not significant. In this respect, the poor quality and low decomposability of beech litter can be associated only in part to the significantly lower R s under hardwood-beech stands compared to hardwoods where there is no significant addition of beech litter. Collin et al. (2017Collin et al. ( , 2018) measured low canopy openness and light transmission under hardwood-beech stands compared to hardwoods and mixedwoods at SBL. For example, canopy openness in July was as low as 10.7% in hardwood-beech stands compared to 23.8 and 30.6% in hardwoods and mixedwoods, respectively. This resulted in light transmission of < 5 mol m −2 d −1 at the soil surface under hardwood-beech stands and > 10 mol m −2 d −1 under the other forest types. Castin-Buchet and Andre (1998) and Kunhamu et al. (2009) suggested that low light transmission can supress microfaunal and microbial activity in the leaf litter by negatively affecting immediate air temperature and relative humidity. At SBL, this is corroborated by significantly lower soil temperature and higher soil water potential in hardwood-beech stands than in other forest types. Differences in soil water potential also appeared to increase during drier periods. In this respect, the lower light transmission in these stands during the full leaf period likely limited overheating of the understory vegetation and soils and in turn, reduced water losses from evapotranspiration. In addition, two results at SBL suggest that differences in R s were due to the lower soil temperature and/or a different response to changes in soil temperature (i.e., greater sensitivity) under hardwood-beech stands compared to other forest types: (1) mixed model analysis with soil temperature as a random variable detected no difference in R s between hardwoods and hardwood-beech stands, and (2) the Q 10 value produced with data from hardwood-beech stands (2.97) was quite lower than Q 10 values produced for hardwoods (2.26) and mixedwoods (2.44). We conclude that low light transmission at the soil surface and low surface soil temperature during the full leaf period were acting in combination with the recalcitrant beech leaf litter in hardwood-beech stands to suppress R h and R s relative to hardwoods without beech or with only small proportions of the species in the canopy.
Annual soil respiration rates under four stands dominated by American beech at Hubbard Brook, New Hampshire, were estimated at 629 g C m −2 y −1 compared to 728 g C m −2 y −1 for three stands dominated by either sugar maple or white birch (Fahey et al., 2005). Although there are some differences in the approach to classify the stands between studies, forests at SBL and Hubbard Brook are quite similar in composition and grow on acidic soils. The main difference is that beech is at approximately 250 km south of its distribution at Hubbard Brook, whereas it is very near its northern limit at SBL. Thus, there are likely similar mechanisms associated with the presence of beech that are decreasing R s rates relative to the other forest types at these sites.

Influence of Soil Water and Magnesium on R s
Soil water potential was also an important effect on R s at SBL as it was the second predictor variable in the multiple regression model. The inclusion of soil water potential in the model is consistent with previous findings. Soil temperature, relative to soil moisture, generally exerts a more consistent control on R s during the growing season in temperate and boreal forests (Maier and Kress, 2000;Gough et al., 2005;Kumpu et al., 2018). However, soil moisture can become an important control of R s during dry periods, namely by affecting the decomposition of leaf litter and thus R h (Epron et al., 1999;Fang and Moncrieff, 2001;Cisneros-Dozal et al., 2007;Vogel et al., 2013;Santonja et al., 2015). In 2020, there were two heatwaves recorded before June 21 and the second one affected southern Quebec and SBL for more than 7 days. This is the first time that two heatwaves were recorded to hit Quebec before the summer solstice. An exceptional drought accompanied these heatwaves, with 70% less rain than normal in June. A shorter heatwave hit Quebec in mid-July and then conditions cooled off and precipitation was normal thereafter. This series of events led to unusually high soil temperature and water potential in June, only about 1 month after complete snow melt at the site, and likely explain the inclusion of an indicator of soil water availability in the multiple regression model.
Soil nutrient availability effects on R s also appeared important at SBL as Mg 2+ activity in the soil solution was selected as the third predictor variable in the multiple regression model. More specifically, Mg 2+ activity in the soil solution may reflect the benefits of an improved soil acid-base status on decomposer biota and activity (Blagodatskaya and Anderson, 1998;Bååth and Anderson, 2003;Han et al., 2008) under hardwoods. Indeed, 4 years of PRS data (2017-2020 in the plots at SBL suggest significantly higher Ca 2+ and Mg 2+ activities under hardwoods than under hardwood-beech stands and mixedwoods (Bélanger, unpublished data). Interestingly, soil solution Ca 2+ activity, another indicator of alkalinity, was selected as the fourth predictor variable in the model (results not shown). The higher activity of these base cations in the soil solution in hardwoods are likely associated with increased biocycling, especially by the abundant Betula spp. in these plots ( Table 1, Bélanger et al., 2004).

CONCLUSION
We assessed R s within a sugar maple forest near the northern limit of deciduous temperate forests and where the abundance of conifers and American beech varied. Seasonal variations in soil temperature exerted the largest influence on R s in these stands, but the results also illustrate how admixtures of balsam fir and beech are affecting R s in opposite directions. The admixture of fir consistently yielded the highest R s , whereas the admixture of beech yielded the lowest R s . The larger R s in the presence of fir is associated to increased photosynthate transport and root activity, whereas the lower R s in the presence of beech is due to the creation of a soil microclimate favored by the dense canopy as well as the poor litter quality (low decomposability). Other factors explaining the variation in R s include low soil moisture during subsequent heatwaves and an improved soil acid-base status due to biocycling by birch trees. Sugar maple was used in this study as the "control" species, but understanding the effects of a change in tree species composition on R s is relevant to all biogeographic contexts. Such studies will help to better predict future temporal and spatial changes in the biogeochemistry of forest ecosystems, including C pools and atmospheric CO 2 concentrations, and their changing trajectories in terms of species composition under climate change (Lafleur et al., 2010;Ettinger and HilleRisLambers, 2013).

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
NB conceived and designed the experiment and drafted the manuscript as the lead author. AC, RK, and SL-D acquired the data and revised and edited the manuscript. NB and AC analyzed the data. All authors contributed to the interpretation of the findings.

FUNDING
Financial support was provided through Natural Sciences and Engineering Research Council of Canada Discovery grants (RGPIN 2015-03699 and 2020-04931) and Canada Foundation for Innovation John R. Evans Leaders Fund (35370) and Innovation Fund (36014) grants to NB.

ACKNOWLEDGMENTS
We wish to acknowledge Charlène Mélançon, Theo Stathopoulos, Anne-Marie Riopel, Dominique Caron, Laurence Grimond, and Justin Bélanger for their help in the field and laboratory. We are also grateful to Mélanie Desrochers for preparing Figure 1. Finally, we thank the staff at SBL for providing access to the site and research facilities.