The studied regions, Kuujjuarapik and Radisson, are located within the discontinuous and sporadic permafrost zones in northwestern Quebéc, Canada (Figure 1). Kuujjuarapik (peatland sites K1P and K2) represents the subarctic, forested tundra ecoregion at the southernmost limit of the discontinuous permafrost zone. K1P peatland is characterized by palsa mounds, whereas K2 is a small fen with peat thickness of ca. 1–2 meters. Radisson (peatland sites LG2 and Rad) represents the northernmost boreal ecoregion. LG2 and Rad are ombrotrophic with peat thickness of ca. 3–4 meters. For more details of the studied peatlands, please refer to Piilo et al. (2019). In both study regions, mean annual temperatures are below 0°C (Table 1). In Kuujjuarapik, mean summer temperature has increased 1°C since 1961, while mean winter temperature has increased 0.5°C since 1971, and the increase in growing degree-days (GDD0) (between 1971–2000 and 1981–2010) is ca. 6%. In Radisson, since 1971, the biggest seasonal increase of 0.5°C is in mean autumn temperatures, and the increase in GDD0 is ca. 4%. In both regions, autumn rainfall especially has increased, and 35–40% of the precipitation falls as snow (Environment Canada, 2019).

The studied 12 peat cores, three from each four peatlands (K2, K1P, LG2, and Rad), were collected in early July 2017 (Figures 1, 2 and Table 1). The core lengths ranged between 32 and 39 cm. Peat sections were sampled by hand, sawing the layer overlying seasonal frost or permafrost from intermediate lawn microforms (water tables ranged from 13 to 23 cm) inhabited mainly by Sphagnum fuscum. Once in the laboratory, the peat cores were cut into contiguous 1 cm slices for further analyses. Plant macrofossil, peat property and carbon accumulation results of these cores were previously published in Piilo et al. (2019).

Testate amoeba analysis was performed at 2 cm resolution. Processing of testate amoeba samples followed a modified version of the standard method (Booth et al., 2010). Samples were simmered in distilled water for 15 min and stirred occasionally. The samples were then sieved with a 300-μm mesh and back sieved with a 15-μm mesh. Materials retained on the 15-μm sieve were centrifuged at 3000 rpm for 5 min. At least 100 individual shells for each sample were counted and identified to species or “type” level under a light microscope at the magnification of 200–400. Taxonomy followed Charman et al. (2000), supplemented with online sources (Siemensma, 2019).

We calculated the community weighted mean values of each standardized trait (Table 1; mixotrophy, biovolume, aperture size and aperture position), which is an index of functional composition expressed as the mean trait value of species present in the community, weighted by their relative abundances (Fournier et al., 2015). In addition, we calculated species richness and Shannon-Wiener diversity index (SWDI) for each sample. Values for test dimensions used to calculate all functional traits with the exception of mixotrophy were those published in the supplementary information of Fournier et al. (2015) except for species that were not included in that study. In those cases, we used average dimensions as described by Siemensma (2019) and in the case of one taxon (Pyxidicula type), an average of multiple replicates of our own light microscopic measurements. Mixotrophy was defined as a binary 0 (non-mixotroph) or 1 (mixotroph), also based on the definitions in Fournier et al. (2015).

Testate amoeba-based water-table depth (WTD) reconstructions were performed using the transfer function developed by Amesbury et al. (2018). The absolute WTD values were normalised to z scores over the length of cores from each region (Swindles et al., 2015b; Amesbury et al., 2016). Z > 0 indicates drier conditions than the region’s average, z < 0 indicates wetter than average. The analysis was conducted using “rioja” package in R version of 3.6.0 (R Core Team, 2019).

The chronologies for each peat record, developed using the “Bacon” package in R and based on ¹⁴C and ²¹⁰Pb dating, were previously published in Piilo et al. (2019). A locally estimated scatterplot smoothing (LOESS) function with a span-value (degree of smoothing) setting of 0.5 was applied to the supra-regional, regional, site and core level data of C accumulation rate, reconstructed WTD and testate amoeba functional trait to explore the temporal trends and spatial variations. The analysis was conducted using the function loess () in R. Linear regression analysis (95% confidence intervals displayed) was applied to C accumulation and testate amoeba functional trait, and reconstructed WTD data using the function lm() in R.

A non-metric multidimensional scaling (NMDS) ordination analysis was applied to testate amoeba functional traits and reconstructed WTD, vegetation types (Piilo et al., 2019) and testate amoeba richness and diversity data to investigate the spatial variations and drivers of the recorded functional trait features. The analysis was performed using “vegan” package in R.

In total, 64 taxa were found in the studied 12 peat records (Supplementary Figure S1). In general, the most commonly recorded taxa in all peat cores were Alabasta militaris, Archerella flavum, Hyalosphenia elegans, and H. papilio. Also, Difflugia pulex was dominant in Radisson cores. In addition, the Pyxidicula type (Supplementary Figure S2), currently unidentifiable to the species level, was found in all the peat records, and in some cases the proportion reached 50%. Site K2 showed higher species richness (average 21) than the other three sites (average around 15) (Figure 3). Sample-specific values for the Shannon-Wiener diversity index (SWDI) were higher at site K2 (average 2.5) than other sites (average 2.1). Temporally, both richness and diversity presented a remarkable decrease starting from ca. 1980 AD at 11 out of 12 records (Supplementary Figure S1). The compositional change of testate amoebae clearly indicated a consistent pattern in all the 12 records, i.e., mixotrophic taxa Hyalosphenia elegans, H. papilio and for some cores Heleopera sphagni became dominant (peak percentage around 50%) in the surface peat sections, even though the timing varied for different cores, ranging from early 1900 to late 1900 AD. Detailed testate amoeba community composition for each peat record can be found in Supplementary Figure S1.

Reconstructed water-table depth showed variations within each site and each region. The moisture conditions for the studied records mainly remained lawn features for the last few centuries. K1P site (especially core K1P3) recorded the largest variations and the wettest conditions (Figures 3, 4). Kuujjuarapik region was wetter than Radisson, especially for the period before 1950 AD (Figure 5), and Kuujjuarapik showed a clear drying starting from early 1900s, while the moisture conditions at Radisson remained relatively stable.

Mixotrophy did not show significant differences between different sites, but a consistent increasing trend toward recent decades (Figures 3–5). The lowest mixotrophy occurred around 1950–2000 and 1900 at regions Kuujjuarapik and Radisson, respectively. Biovolume and aperture size were significantly higher at sites LG2 and Rad than K2 and K1P. Biovolume showed a general decreasing trend toward the present, even though some records showed an increasing trend (LG2.1 and Rad1). Aperture position at site LG2 was lower than other three sites. Regionally, aperture size and position at Radisson showed a decreasing trend toward the present, while at Kuujjuarapik first increased until 1950 AD and then started to decrease. In addition, aperture size and position were much higher at Radisson than Kuujjuarapik starting from 1900 to 2000 AD.

Functional trait, C accumulation and reconstructed WTD results at different spatial scales showed differing degrees of variability (Figure 5). C accumulation and mixotrophy results were relatively consistent across all spatial scales, showing recent increases. Reconstructed WTD values tended to show a relatively consistent trend toward drier conditions, but the regional level data showed this was a more recent and sharper trend in the Kuujjuarapik region, compared to a more gradual trend in the Radisson region. A similar pattern is evident in the aperture size and aperture position data: both showed consistent decreasing trends using all data, but regional level sub-division showed that the trends were more severe and recent in the Kuujjuarapik region and more gradual in the Radisson region. Biovolume tended to decrease over the study period, but regional level data showed that biovolume was generally higher in the Kuujjuarapik compared to Radisson region and core level data showed that the trend toward lower biovolume was much stronger in some cores, whereas others displayed relatively little change.

There were significant correlations between C accumulation and testate amoeba functional trait data, but the patterns were not consistent (Figure 6). For example, despite qualitatively similar recent increases in both C accumulation and mixotrophy across many cores (Figure 5), correlations between these factors were significant in both positive (K2.1, Rad1) and negative (K2.3, LG2.2) directions, even at the same site (K2) (Figure 6). Correlations between aperture size and aperture position were generally more consistent, with 3 of 4 and 6 of 7 significant correlations, respectively, between C accumulation and these functional traits across all sites being negative (Figure 6).

The NMDS results, i.e., no clear clusters of samples from a certain site, showed that the testate amoeba functional traits in each peatland were relatively homogeneous (Figure 7). Mixotrophy was highly linked to the first axis and clearly preferred habitats with Sphagnum mosses. Aperture size located on the other side of the first axis showed a shorter length, i.e., smaller aperture in Sphagnum-dominated habitats. The second axis can be explained by water table conditions, with the upper part indicating dry conditions with abundant woody species and the lower part indicating wetter conditions with sedges and non-Sphagnum mosses and higher testate amoeba richness and diversity. Biovolume was closely linked to the second axis, with larger values occurring in wetter conditions. Aperture position with hidden apertures showed a preference to drier conditions with abundant woody plants.

Changes in response traits at the community scale are inferred as being driven by changes in environmental conditions, while changes in effect traits may indicate shifts in process at the ecosystem scale. Previous studies have shown that correlations between response traits and effect traits can be used to determine how environmental change will influence ecosystem functioning (Lavorel and Garnier, 2002; Suding et al., 2008; Hillebrand and Matthiessen, 2009).

Biovolume and aperture position are suggested to be correlated with moisture, with low water level promoting smaller biovolume and hidden and protected apertures (Fournier et al., 2015; Marcisz et al., 2016). Our data strongly supported such links, especially for biovolume; that higher biovolume highly correlated with wetter conditions. Unlike biovolume, in addition to water table, aperture position might be influenced by other factors like the abundance of woody plants; for example, Rad1 had protected aperture (higher values of aperture position) in response to drier conditions, while for more records, e.g., LG2.2, Rad2, Rad3, K1P1, K1P2, K2.2, and K2.3, less protected aperture increased under wetter conditions. This suggests that variations in water level may not be a decisive factor for changing aperture position. Also, it may be explained by our field sampling design, i.e., all the peat cores were collected from lawn microforms with relatively good water availability, so the detected variations in water table were not strong enough to drive aperture position composition changes. Previous studies have also shown that larger disturbances, like flooding, fire, very low water tables and deforestation might have stronger impacts on such trait compositions (Marcisz et al., 2016, 2019). While dry conditions that supported woody plant growth was not a driver for aperture position, the altered vegetation composition (i.e., more woody plants) might nevertheless be challenging for testate amoebae to adapt to (Marcisz et al., 2016 and references therein).

Mixotrophy, seem to show relatively consistent increasing trends over recent decades, which are in line with the increased carbon accumulation rates, even though these increases of C accumulation rates might partly due to less decomposition of recently accumulated peats (Piilo et al., 2019). Mixotrophic testate amoebae that acquire most of their C by photosynthesis of their symbionts, essentially functioning as autotrophs (Lara and Gomaa, 2017) are suggested to play an important role in driving peatland C accumulation, with higher abundances promoting C accumulation (Marcisz et al., 2016; Zhang et al., 2020), but the opposite pattern can also be expected depending on the specific condtions. For example, in this study both positive and negative linear links between mixotrophic testate amoebae and C accumulation were derived (Figure 6), which might depend on the prey availability, i.e., mixotrophic testate amoebae either favour photosynthesis and form a positive link between C accumulation and their abundance, or predation and constitute a negative link between C accumulation and their abundance (Jassey et al., 2012, 2013; Mieczan et al., 2015). Furthermore, in most of our cases the larger the aperture size, the lower the C accumulated, which supports the latter mechanism. Although it has been shown before that lowering of water table leads to a decrease in the abundance of mixotrophs (Marcisz et al., 2016), this is not visible in our records. This might also be due to the sample collecting habitats, i.e., lawns that supply sufficient water and did not reach the threshold for changes of mixotrophy (Jassey et al., 2018; Koenig et al., 2018; Lamentowicz et al., 2019; Zhang et al., 2020). Instead, we found that increases in mixotrophy were linked to increased Sphagnum (mainly S. fuscum) abundance (Van Bellen et al., 2018). This suggests that mixotrophs appear to thrive under increasingly acidic and ombrotrophic conditions (Fournier et al., 2015). In addition, higher light intensity at the surface could be another potential reason for increased near surface mixotrophy (Marcisz et al., 2014; Payne et al., 2016; Creevy et al., 2018), even though increase of mixotrophy might also occur under low light intensity conditions (Herbert et al., 2019). Aperture size that indicates the trophic position of testate amoebae in the microbial food web (Jassey et al., 2013) was found to mainly present a negative correlation to C accumulation, however, an adverse pattern was also captured (Figure 6). Similar to mixotrophy, aperture size also showed no links to water table, but a negative response to increase of Sphagnum, with smaller size observed in Sphagnum-dominated habitats.

As discussed above, mixotrophic testate amoebae might increase or decrease peatland C accumulation depending on different biotic conditions. Therefore, changes in their abundance will have impacts on peatland C dynamics. Our results show that the lowest abundance (close to zero) of mixotrophic testate amoebae occurred when Sphagnum dominance remarkably decreased. Post-LIA warming seems to play an important role in driving such vegetation change, followed by increasing of Sphagnum (Piilo et al., 2019) and correspondent mixotrophic testate amoebae, thus testate amoebae-peatland functioning. Particularly, the palsa site (K1P) showed synchronous vegetation changes, characterised by decreased Sphagnum abundance, for all records starting from ca. cal. AD 1950s and lasting until mid-1990s, which could be linked to changes in thawing permafrost dynamics in the site (Piilo et al., 2019). Increased snowfall since the 1950s and warm temperatures in the 1990s triggered rapid permafrost thawing on the eastern coast of Hudson Bay. Hence, permafrost is predicted to disappear from subarctic Québec in the coming decade (Payette et al., 2004), likely causing drying or wetting of the peat surface (Swindles et al., 2015a; Zhang et al., 2018b) and altering the vegetation (Zhang et al., 2018b; Piilo et al., 2019) and mixotrophic testate amoeba abundance, and ultimately changing peatland carbon dynamics.

In addition, boreal peatlands in, for example, western Canada are highly sensitive to fire because of the dry climate and high tree density, which sustains the spread of fire. Enhanced warming and drying during the 21st century will result in more intense regional fire regimes (De Groot et al., 2013). Fires were captured in 11 out of 12 studied records, especially in peatland Rad (Piilo et al., 2019). Atmospheric deposition of mineral dusts, whether they are supplied by natural or anthropogenic sources, may affect testate amoeba communities and their functional traits (Fialkiewicz-Koziel et al., 2015; Lamentowicz et al., 2015). For example, enhanced inputs of atmospheric dust are likely to limit mixotrophic testate amoeba abundance and therefore C uptake. However, the effects of drivers such as dust deposition are likely to have less effects on changes of testate amoeba taxa when compared with vegetation change, which partly depends on the volume of dust (Van Bellen et al., 2018), which seems to be supported also by our records that fires did not trigger clear shifts on testate amoeba composition.

Warming could enhance peatland C uptake due to increased growing season length (Charman et al., 2013; Gallego-Sala et al., 2018), more specifically increased plant photosynthesis. However, in terms of photosynthetic microorganisms, warming and subsequent changes in hydrology and vegetation all have impacts on mixotrophic testate amoebae. Studies have shown a reduction in mixotrophic testate amoebae with rising temperature (Wilken et al., 2013; Jassey et al., 2015) and lowering water tables (Marcisz et al., 2016; Basińska et al., 2020). However, one hypothesis is that global warming will lead into water level drawdown in peatlands, which favours new establishment of Sphagnum (Tahvanainen, 2011; Magnan et al., 2018). This Sphagnum-related change would likely to promote mixotrophic testate amoeba presence and mitigate the decrease of mixotrophic testate amoebae caused by warming and lowering of water table. On the other hand, drying might also cause shrub expansion (Gałka et al., 2017a; Zhang et al., 2018b) and decrease the openness of the peatland, thus hinder peatland C uptake through decreases in mixotrophic testate amoebae (Payne et al., 2016; Lamentowicz et al., 2020). Nevertheless, future climate change might influence peatland C sequestration through not only directly altering the hydrology, but also vegetation and consequent changes in mixotrophic testate amoebae, even though the magnitude and direction of such changes are not well constrained and might vary from one ecoclimatic region to another.

The spatial and temporal variations detected within both peatlands and regions highlight the importance of studying multiple peat sections from one study site and preferably the approach should be extended to regional-scales (Loisel and Garneau, 2010; Lamarre et al., 2012; Mathijssen et al., 2017; Zhang et al., 2018a, b; Piilo et al., 2019). To enable a spatial comparison of our records, we collected peat cores consistently from lawn habitats, however, even in a particular site large differences between the three replicated cores were captured. It has been previously reported that different peatland habitats may also undergo different peatland-climate feedback pathways (Zhang et al., 2020). Any one site-based discussion of regional feedbacks is likely to be biased because different peatland types can experience specific successions under similar climate conditions, such as the differences observed between palsa K1P and fen K2 in region Kuujjuarapik in this study, let alone different regions that have different climate (Sim et al., 2019). This problem is exacerbated when within-site data, even taken from a consistent microform, is also variable. As one main aim of small-scale studies is to provide information for global modeling of peatland-climate feedbacks, we need to keep it in mind that any single core is likely not represent the overall site or regional response. Consistent signals detected from replicated records are definitely more convincing in terms of interpretations on peatland-climate feedbacks. However, due to practical issues for particular studies, we do not provide any specific suggestions here about, for example, the number of records that are needed for a reliable signal, but emphasize the need for multiple cores in peatland-climate feedback studies.