Dispersal of bacteria and stimulation of permafrost decomposition by Collembola

Contrary to most soils, permafrost soils have the atypical feature of being almost entirely deprived of soil fauna. Abiotic constraints on the fate of permafrost carbon after thawing are increasingly understood, but biotic constraints remain scarcely investigated. Incubation studies, essential to estimate effects of permafrost thaw on carbon cycling, typically measure the consequences of permafrost thaw in isolation from the topsoil and thus do not account for the effects of altered biotic interactions because of e.g. colonization by soil fauna. Microarthropods facilitate the dispersal of microorganisms in soil, both on their cuticle (ectozoochory) and through their digestive tract (endozoochory), which may be particularly important in permafrost soils, considering that microbial community composition can strongly constrain permafrost biogeochemical processes.

Here we tested how a model species of microarthropod (the Collembola Folsomia candida) affected aerobic CO₂ production of permafrost soil over a 25 d incubation. By using Collembola stock cultures grown on permafrost soil or on an arctic topsoil, we aimed to assess the potential for endo- and ectozoochory of soil bacteria, while cultures grown on gypsum and sprayed with soil suspensions would allow the observation of only ectozoochory.

The presence of Collembola introduced bacterial amplicon sequence variants (ASVs) absent in the no-Collembola control, regardless of their microbiome manipulation, when considering presence–absence metrics (unweighted UniFrac metrics), which resulted in increased species richness. However, these introduced ASVs did not induce changes in bacterial community composition as a whole (accounting for relative abundances, weighted UniFrac), which might only become detectable in the longer term.

CO₂ production was increased by 25.85 % in the presence of Collembola, about half of which could be attributed to Collembola respiration based on respiration rates measured in the absence of soil. We argue that the rest of the CO₂ being respired can be considered a priming effect of the presence of Collembola, i.e. a stimulation of permafrost CO₂ production in the presence of active microarthropod decomposers. Overall, our findings underline the importance of biotic interactions in permafrost biogeochemical processes and the need to explore the additive or interactive effects of other soil food web groups of which permafrost soils are deprived.

1 Introduction

Carbon fluxes from soils are largely governed by the rate of decomposition of organic matter. Soil fauna is a crucial component in organic matter decomposition (García-Palacios et al., 2013; Griffiths et al., 2021), both directly through mechanical degradation of litter and corpses into smaller pieces and indirectly through feeding habits controlling the abundances of other decomposer groups such as fungi, microbial eukaryotes or bacteria (Hanlon and Anderson, 1979; Kaneda and Kaneko, 2008; Frouz et al., 2020; Potapov et al., 2020). In arctic soils, the scarcity of macrofaunal decomposers (e.g. earthworms – Blume-Werry et al., 2020) coupled with the high abundance of microbivorous microarthropods such as Collembola (Potapov et al., 2022) results in a particularly strong impact of Collembola on decomposition through microbial population control (Koltz et al., 2018; Crowther et al., 2012; Seastedt and Crossley, 1984).

Contrary to other important groups of soil fauna, such as earthworms, millipedes (Golovatch and Kime, 2009; Berman et al., 2015) or woodlice (Sfenthourakis and Hornung, 2018), Collembola are ubiquitous in arctic soils, where they can reach high densities – up to 130 000 individuals per square metre in high-Arctic Greenland (Sørensen et al., 2006). However, Collembola are mostly abundant in the topsoil and to our best knowledge have never been observed in the perennially frozen subsoil, the permafrost. Permafrost soils are a prominent feature of Arctic landscapes, and the huge carbon stock they represent and positive feedback to climate change that their thawing will likely induce cause concern (Masson-Delmotte et al., 2021). The frozen conditions over long periods of time have eliminated most fauna and non-microbial life from these environments, and despite the possibility of resuscitating organisms such as plants (Yashina et al., 2012), nematodes (Shatilovich et al., 2018) or rotifers (Shmakova et al., 2021), newly thawed permafrost soils generally harbour an extremely simplified soil food web entirely deprived of metazoans. Despite this particularity, the absence of Collembola and the consequence of their possible introduction into newly thawed permafrost on its biogeochemical cycling have been mostly overlooked.

In contrast to Collembola, some microorganisms survive and/or thrive in permafrost conditions, and microbial adaptations to frozen conditions have been studied (Mackelprang et al., 2011, 2017; Hultman et al., 2015; Bottos et al., 2018). However, not all microbes survive these conditions, and the combination of environmental constraints exerted over long periods of time (Mackelprang et al., 2017) and strong dispersal limitations (Bottos et al., 2018) results in microbial communities that can be deprived of some functions (Knoblauch et al., 2018; Monteux et al., 2020; Barbato et al., 2022). The re-introduction of such functions can result in drastic changes in permafrost processes, and sizable impacts on greenhouse gas production have been observed in vitro for CH₄ and CO₂ (Knoblauch et al., 2018; Monteux et al., 2020) and confirmed in situ for N₂O (Marushchak et al., 2021). Upon thawing, this re-introduction of missing functions or ecological rescue (Calderón et al., 2017) requires microorganisms to migrate into this newly available habitat, which could happen for instance laterally through airborne dispersal (Harding et al., 2011) – e.g. for permafrost exposed to the air in abrupt thaw processes (Inglese et al., 2017) – or vertically through percolation in the soil column when the active layer becomes deeper (Monteux et al., 2018; Johnston et al., 2019).

Collembola are another possible vector for the dispersal of microorganisms into newly thawed permafrost (Buse et al., 2014). With the deepening of the active layer, Collembola migration into newly thawed permafrost is unlikely since they mostly reside in topsoil layers, but it will likely occur in soil-mixing events (Väisänen et al., 2020) such as thermokarst, active-layer detachment or thaw slump processes where newly thawed permafrost is exposed to surface conditions. Like all large organisms, Collembola host a variety of microorganisms, their microbiome (Agamennone et al., 2015; Leo et al., 2021). Because microarthropods, such as Collembola, can move across large distances compared to fungi or bacteria, the microbial species in collembolan microbiomes might be among the first to colonize and establish in newly thawed permafrost when Collembola access it. This could occur through a combination of two main processes: ectozoochory, where the microbiome of the cuticle disperses into the new habitat, and endozoochory, where microorganisms disperse after transiting through the gut of the animals. To our best knowledge, whether Collembola affect the biogeochemical functioning of newly thawed permafrost and whether and how they can serve as a vector for microbial colonization have not been explored yet.

We incubated permafrost from the Yedoma domain, which represents a large carbon stock in parts of Siberia and North America (Strauss et al., 2017) and has previously been shown to lack certain microbial functions (Monteux et al., 2020), in the presence or absence of a model species of microarthropod (Collembola Folsomia candida), and assessed CO₂ production and bacterial community composition over a 25 d aerobic incubation. We used Folsomia candida Collembola from a stock culture, as well as Collembola subjected to manipulation of their cuticle microbiome or cuticle and gut microbiomes to test the following hypotheses:

1. Collembola presence in permafrost will alter bacterial community composition, through their grazing.

2. Collembola more closely exposed to topsoil bacteria will change permafrost bacterial community composition further than those exposed less intricately. In other words, we expect a gradual change between clean Collembola, Collembola with cuticle microbiome manipulated (ectozoochory), and Collembola with both gut and cuticle microbiomes manipulated (endozoochory).

3. Collembola presence in permafrost will increase CO₂ production, through both Collembola respiration and a stimulation of microbial activity (“priming effect”).

4. The gradual introduction of distinct bacterial communities hypothesized above (2) will result in increased CO₂ production, owing to the functional limitations of the Yedoma permafrost microbial community in terms of CO₂ production.

2 Materials and methods

2.1 Experimental design

2.1.1 Soils

The Yedoma sediment used to assess the impacts of Collembola on permafrost CO₂ production and bacterial community composition originated from the Cold Regions Research and Engineering Laboratory (CRREL) permafrost tunnel (Fox, Alaska, USA). The sediment was sampled from the upper silt unit and is an Upper Pleistocene silty deposit, previously described in detail (Shur et al., 2004; Mackelprang et al., 2011, 2017; Monteux et al., 2020). This sediment was chosen due to its microbial communities being vulnerable to invasions and exhibiting functional limitations, allowing the discernment of impacts of introduced microorganisms on broad proxies such as CO₂ production (Monteux et al., 2020). Approximately 35 g (fresh weight) of homogenized sediment was set in 200 mL glass jars, sealed with parafilm to allow for gas but not moisture or microorganism exchange, and pre-incubated at 10 ^∘C for 11 d before inoculation.

To manipulate the Collembola microbiome and make it more similar to that found in natural settings, we collected a topsoil (0–15 cm depth) from a subarctic meadow (Kärkevagge, 30 km west of Abisko, northern Sweden; 68^∘24^′23.8^′′ N, 18^∘18^′51.6^′′ E) in September 2019 and kept it frozen until the cultivation of Collembola. We could not obtain topsoil from the same location as the Yedoma permafrost due to practical constraints during the Covid-19 pandemic.

2.1.2 Collembola

A strain of the Collembola species Folsomia candida was obtained from Vrije Universiteit Amsterdam and cultured 6 months prior to the onset of the experiment. Folsomia candida is a parthenogenetic ground-dwelling Collembola, which has been routinely used as a model organism in soil ecology. Stock cultures were maintained on a gypsum and coal medium and fed baker's yeast; traces of mould were removed; fresh yeast and water were added once to twice a week, and fresh stock cultures were started monthly.

2.1.3 Collembola inoculation treatments

Two months prior to the onset of the experiment, separate stock cultures were established on gypsum and coal medium supplemented with a 2–3 cm layer of topsoil to obtain Collembola whose skin and gut microbiome were both colonized with topsoil microorganisms (“topsoil stock culture”). In parallel, similar stock cultures using permafrost sediment were established as an additional control (“permafrost stock culture”). These stock cultures with soil or sediment were supplemented with yeast and water like the gypsum stock culture to maintain high adult population densities prior to the experiment. One day before inoculating the incubation jars, soil suspensions were made from the topsoil and permafrost sediment (5 and 10 g, respectively, in 100 mL ddH₂O), shaken at 150 rpm for 1 h and filtered (Ahlstrom-Munksjö grade 006, 1.5 µm pore size) to manipulate the skin microbiome of Collembola in the ectozoochory treatments.

All incubation jars containing permafrost were randomly assigned a treatment on the day of inoculation from among the following, replicating each treatment across six jars (Fig. 1). The treatments consisted in adding to the jars with permafrost

• no Collembola (hereafter, “no-Collembola control”),

• Collembola grown on gypsum stock culture (“Collembola”),

• Collembola grown on gypsum stock culture and sprayed with topsoil suspension (“ectozoochory”),

• Collembola grown on gypsum stock culture and sprayed with permafrost suspension (“ectozoochory control”),

• Collembola grown on topsoil stock culture (“endozoochory”),

• Collembola grown on permafrost stock culture (“endozoochory control”).

To obtain a baseline allowing the quantification of the effect of Collembola on permafrost CO₂ production and bacterial community composition, Yedoma permafrost was incubated in the absence of Collembola (no-Collembola control). For all other treatments, the permafrost was supplemented with adult Collembola isolated from the stock cultures using a handheld vacuum cleaner mounted with a 10 mL pipet tip. The Collembola were transferred into a black plastic tray, allowing them to be spread and the adults to be picked out, excluding juveniles as much as possible. From this tray, a similar number of Collembola (30–80 individuals, in the range of values used in the literature, e.g. by Kaneda and Kaneko, 2008) was sampled with the vacuum cleaner and then inoculated into each jar by pouring them into a plastic funnel, with the Collembola provenance or manipulation depending on the treatment. Prior to transferring them into incubation jars, Collembola used in the ectozoochory treatment and its control were sprayed with the corresponding soil suspensions. To limit cross-contamination, separate funnels, as well as separate 10 mL pipet tips on the vacuum cleaner, were used for the different treatments. A picture of the inside of each jar was taken to count the exact number of Collembola, then the jars were closed with rubber septa and flushed with moisturized CO₂-free air and incubated. To estimate Collembola respiration per individual, Collembola were also incubated in six jars under identical conditions, except that the soil was replaced with a few drops of autoclaved ddH₂O to prevent dehydration (hereafter, no-soil calibration).

Figure 1 Summary of the experimental design. Folsomia candida Collembola were grown in a stock culture on gypsum and fed with yeast and then subjected to different treatments before being inoculated into jars filled with Yedoma permafrost sediment (30–80 individuals per jar). Ectozoochory was assessed by spraying the Collembola with soil suspensions, while endozoochory was assessed by growing Collembola in stock cultures supplemented with soil. The jars were incubated in the dark for 25 d at 10 ^∘C; CO₂ production was measured throughout the incubation period, and bacterial community composition was determined at the end of the incubation.

2.1.4 Incubation

We dark-incubated all flasks for 25 d under aerobic conditions at 10 ^∘C. This incubation temperature is similar to summer active-layer temperatures in permafrost-affected areas and within the thermal tolerance range of psychrophilic microorganisms (D'Amico et al., 2006). We used a short (25 d) incubation period to ensure a relatively stable Collembola population level, by limiting uncertain numbers of newly hatched Collembola individuals, since the eggs of Folsomia candida take 18–20 d to hatch at 16 ^∘C (Marshall and Kevan, 1962) and presumably longer than 25 d at 10 ^∘C.

2.2 Measurements

2.2.1 CO₂ production

Headspace air was sampled with a syringe to measure CO₂ concentrations (EGM-5 IRGA, PP Systems, Amesbury, Massachusetts, USA) at intervals ensuring CO₂ concentrations remained below 20 000 ppm to prevent a toxic CO₂ build-up (i.e. after 3, 7, 14 and 25 d). After each measurement, the jars were flushed with 0.45 µm filtered CO₂-free air moisturized by bubbling it through two 5 L bottles of ddH₂O, for 3 min at 1 to 2 L min⁻¹, i.e. with at least 15 times the volume of the jar. CO₂ concentrations were adjusted for changes in temperature and atmospheric pressure to calculate CO₂ production rates (τ) as follows:

$$\begin{array}{}\text{(1)}& {\mathit{\tau }}_{\left(i\right)}={\displaystyle \frac{{\left[{\mathrm{CO}}_{\mathrm{2}}\right]}_i\times \left(P_{i}V/R{T}_i\right)}{{\left({\mathrm{\Delta }}_t\right)}_i}},\end{array}$$

where (Δ_t)_i is the time interval between measurement (subscript i) and previous flushing, P_i atmospheric pressure at measurement time, V the headspace volume, R the ideal gas constant, and T_i the temperature. To calculate cumulative CO₂ production over the entire incubation, we summed up the quantity of CO₂ present in the headspace at each sampling, within each jar.

2.2.2 Bacterial community

Microcentrifuge tubes (1.5 mL) were filled with soil and snap-frozen in dry ice to analyse microbial communities from the jars harvested at the end of the incubation. The frozen tubes were kept at −20 ^∘C for up to 4 months before freeze-drying and then homogenized by bead beating (Precellys CK68 15 mL tubes, 2×30 s at 4500 rpm). DNA was extracted from 183 to 285 mg of homogenized freeze-dried soil using the DNeasy PowerSoil Pro Kit (Qiagen) according to the manufacturer's instructions, and DNA concentrations in the extracts were measured on a Qubit 1.0 fluorometer.

The V4–V5 region of the 16S ribosomal RNA gene was targeted in PCR amplification using primers 515F (5^′-GTGYCAGCMGCCGCGGTAA-3^′) and 926R (5^′-CCGYCAATTYMTTTRAGTTT-3^′) with Illumina sequencing adapters, using 12.5 µL Phusion Taq Green PCR Master Mix (Thermo Scientific), 0.25 µM of each primer, 2 µL of DNA extract diluted to 5 ng µL⁻¹ and nuclease-free water in 25 µL reaction volume. PCR conditions were as follows: initial denaturation (98 ^∘C, 3 min), 25 cycles of denaturation (98 ^∘C, 15 s), annealing (50 ^∘C, 30 s) and elongation (72 ^∘C, 40 s), and a final elongation (72 ^∘C, 10 min), after which PCR products were checked by electrophoresis on 1 % agarose SB gel. A total of 20 µL of each PCR product was cleaned, and the products' DNA concentrations were normalized using a SequalPrep Normalization Plate Kit (Invitrogen), according to manufacturer's instructions. Three DNA extraction blanks and two PCR blanks were included as negative controls, as well as a mock community as the positive control (ZymoBIOMICS, diluted to 5 and 0.5 ng µL⁻¹ in two replicates each).

A second PCR step was performed to add Nextera dual-indexing barcodes, using 30 µL reaction volume, 1 µM of each primer and 5 µL of cleaned PCR product. PCR conditions were as follows: initial denaturation (98 ^∘C, 3 min), eight cycles of denaturation (98 ^∘C, 30 s), annealing (55 ^∘C, 30 s) and elongation (72 ^∘C, 40 s), and a final elongation (72 ^∘C, 10 min), after which PCR products were checked by electrophoresis on 1 % agarose SB gel. A total of 25 µL of each PCR product was cleaned, and the products' DNA concentrations were normalized using a SequalPrep Normalization Plate Kit (Invitrogen). Serial elution across columns was used to increase concentration of the pooled products, i.e. using only 8×20 µL elution buffer instead of 96×20 µL. The eluted DNA was pooled; its concentration was measured on a Qubit fluorometer, and the size distribution of the amplicons was measured by automated electrophoresis (Agilent 2100 Bioanalyzer). The library was then sent for sequencing on an Illumina MiSeq with V3 chemistry (2×300 bp, 15 % PhiX spike-in) at the SNP&SEQ Technology Platform in Uppsala. Demultiplexing was performed by the sequencing facility, and data were deposited at ENA with accession number PRJEB51992.

2.3 Data analysis

2.3.1 Bioinformatics

All bioinformatics and statistics were performed in R v4.1.3 (R Core Team, 2022), unless specified otherwise. The whole analysis pipeline is found at the Bolin Centre Code Repository (Monteux, 2022) and the processed data and figure-generating script at Zenodo (Monteux et al., 2022). In short, amplicon sequence variants (ASVs) were created with DADA2 1.18.0 (pseudo-pooling; Callahan et al., 2016) after removing primers and adapters with Cutadapt (v3.10; Martin, 2011). Taxonomy was assigned to ASVs with the RDP naïve Bayesian classifier (v1.8; Wang et al., 2007), and ASVs resolved to the genus rank were further assigned a species rank by the exact string-matching algorithm implemented in DADA2 (assignSpecies), using SILVA v138.1 reference data (Quast et al., 2013). Putative contaminant ASVs were manually selected from those identified in silico using the decontam algorithm (Davis et al., 2018) with combined prevalence- and frequency-based methods using the default threshold of 0.1 and separate prevalence- and frequency-based methods with a threshold of 0.05. Eight contaminant ASVs amounting to up to 0.024 % of the total reads were removed. Appropriateness of the bioinformatics analysis parameters was judged by visually assessing the composition of mock communities (microbial community DNA standard, ZymoBIOMICS) at the genus level, leading to ASVs amounting to up to fewer than 10 reads and/or present in fewer than 3 samples being removed from the dataset. Sequencing depth was not associated with experimental treatments (ANOVA $F_{\mathrm{5},\mathrm{30}}=$ 0.701, P= 0.627); therefore read numbers were converted to proportional abundances within samples to normalize sample sizes.

2.3.2 Diversity analyses

The effect of the different treatments on bacterial communities was visualized using principal coordinates analysis (PCoA) and tested with permutational multivariate analyses of variance (PERMANOVAs, adonis function in the vegan package; Oksanen et al., 2022) after verifying homoscedasticity (betadisper). Two distance matrixes were computed for that purpose, using weighted and non-weighted UniFrac distances (Lozupone et al., 2011) to distinguish between compositional effects accounting for bacterial relative abundances and for only presence–absence, respectively. Pairwise contrasts were subsequently computed using the wrapper provided in the pairwiseAdonis R package (Arbizu, 2022).

In addition, the same effect was tested with the more robust and sensitive manyglm approach (Wang et al., 2012), which avoids certain pitfalls of distance-based methods (Warton et al., 2012). The manyglm models were fitted with negative binomial distribution, after visually checking that assumptions were met, on the non-normalized ASV count data using the default PIT-trap resampling with 1999 bootstrap permutations and a likelihood-ratio testing method. An analysis of deviance was carried out (anova.manyglm), subsequently using the provided wrapper for pairwise comparisons with free step-down P-value adjustment.

2.3.3 CO₂ production analysis

A repeated-measures ANOVA was carried out, using Greenhouse–Geisser ϵ correction to degrees of freedom to account for the violation of the sphericity assumption and thus to assess the interactive effect of the Collembola manipulation treatments over time. Since this interactive effect was not statistically significant (P= 0.56, Table A1), we removed the time dimension and used cumulative CO₂ production at the end of the incubation period in further analyses.

To assess our hypothesis that Collembola presence would increase CO₂ production we used two-sample t tests with unequal variances to compare the cumulative CO₂ production at the end of the incubation between the no-Collembola control and all other jars.

To explore differences between the Collembola microbiome manipulations, we used a one-way ANOVA followed by treatment contrasts to assess the difference from the control (i.e. all treatment and controls compared to the no-Collembola control) and selected orthogonal contrasts to assess the effects of each treatment (i.e. each treatment compared to its own control) using the “emmeans” package (Lenth, 2016).

2.3.4 Estimated Collembola respiration

Based on photographs taken upon inoculating the incubation jars with Collembola and 3 d after inoculating, we counted the exact number of live Collembola (i.e. ignoring Collembola which had apparently not moved between when the two pictures were taken). Only a few animals did not survive the transfer, and no dead animals were observed at the end of the incubation. The number of Collembola per jar varied across treatments ($F_{\mathrm{4},\mathrm{25}}=$ 3.19, P= 0.030), although no pairs of treated jars significantly differed from each other (when including the no-soil calibration jars, $F_{\mathrm{5},\mathrm{30}}=$ 4.048, P= 0.006; Fig. B1); we therefore needed to account for differing Collembola numbers. We averaged the respiration per individual in the no-soil calibration set for each of the CO₂ concentration measurement times and multiplied this amount by the number of Collembola present in each jar to estimate the amount of CO₂ produced by Collembola basal respiration in treatment jars.

2.3.5 Response-ratio calculations

A response ratio of CO₂ production in Collembola treatments was calculated by dividing cumulative CO₂ production by the average of the no-Collembola control. This was performed both for the net CO₂ production after subtracting the estimated Collembola basal respiration (RR_soil) and for the gross CO₂ production (RR_gross). This allows the partition of the difference in altered CO₂ production between what was respired by the Collembola and a putative priming effect on soil CO₂ production.

3 Results

3.1 Bacterial communities

Bacterial community composition was overall largely unaffected by the presence of Collembola. Taking relative abundances into account, no effect of the different treatments was identified with PERMANOVA (P > 0.05, Fig. 2a), but a significant treatment effect was observed using the presence–absence distance metrics (PERMANOVA $F_{\mathrm{5},\mathrm{30}}=$ 2.91, P= 0.001, Fig. 2b). More specifically, treatments containing Collembola differed from the no-Collembola control (Fig. 2b) and to some extent from each other (Table C1). Using manyglm, a significant treatment effect was observed (manyglm analysis of deviance P < 0.001; Table 1). The ectozoochory treatment and its control (Collembola sprayed with topsoil or permafrost soil suspension, respectively) differed significantly from the no-Collembola control in the post hoc test (0.01 < P < 0.05), while the other treatments including Collembola tended to differ from the no-Collembola control (0.05 < P < 0.1, Table 1). Similarly to the presence–absence analysis, treatments that included Collembola did not significantly differ from each other in the manyglm analysis. Overall, soil with Collembola exhibited higher alpha-diversity metrics (richness estimators Chao1, abundance-based coverage estimator – ACE, number of observed ASVs, as well as Shannon and Fisher diversity indexes; Fig. 3, Table C2) than the no-Collembola control, but those did not differ among the Collembola treatments.

Figure 2 Response of permafrost bacterial communities to the addition of Collembola undergoing different microbiome manipulation treatments. (a) Principal coordinates analysis (PCoA) of normalized abundance of bacterial ASVs (amplicon sequence variants) based on weighted UniFrac phylogenetic distance (i.e. accounting for relative abundances); (b) PCoA of normalized abundance of bacterial ASVs based on non-weighted UniFrac phylogenetic distance (i.e. presence–absence only). Ellipses indicate the 80 % confidence interval around treatment centroids, and axes are scaled to the percentage of explained variance. The asterisk indicates that the no-Collembola control (dark brown) in (b) significantly differs from all other treatments (pairwise PERMANOVA P < 0.05) and endo- and ectozoochory treatments did not differ from the Collembola treatment. Pairwise comparisons are summarized in Table C1.

Figure 3 Response of permafrost bacterial alpha diversity to the addition of Collembola undergoing different microbiome manipulation treatments. (a) Number of observed ASVs (amplicon sequence variants) per sample; (b) abundance-based coverage estimator (ACE). Coloured symbols and error bars are means and standard errors within a treatment (n=6); small grey symbols are individual samples. Error bars are only shown when exceeding the symbols' size. Asterisks denote that the no-Collembola control differs from all other treatments (estimated marginal means pairwise comparison, 0.001 < Holm-adjusted P < 0.01).

3.2 CO₂ production

Overall, the Collembola addition resulted in higher CO₂ production than the no-Collembola control, but the ectozoochory and endozoochory treatments did not result in higher CO₂ production than their respective controls (using gross CO₂ production, as in RR_gross; Table A2, Fig. A1). When accounting for Collembola basal respiration (as in RR_soil), only the ectozoochory treatment (P= 0.005) and to a lesser extent its control (P= 0.069) differed from the no-Collembola control (Table A2, Fig. A1).

We hypothesized increased CO₂ production in the zoochory treatments through an effect on bacterial communities; however that effect was absent. We therefore tested for the overall effect of Collembola presence across all treatments compared to the no-Collembola control (Fig. 4). Using a t test with unequal variances, we observed a 25.85 % increase in CO₂ production in the presence of Collembola (RR_gross; Fig. 4a). When subtracting estimated Collembola basal respiration, this increase was on average 13.22 % (RR_soil; Fig. 4b); thus roughly half of the observed increase in CO₂ production could be attributed to Collembola basal respiration.

Figure 4 Response ratios of permafrost cumulative CO₂ production rates to addition of Collembola. (a) RR_gross: measured CO₂ production rates. (b) RR_soil: measured CO₂ production rates with estimated Collembola basal respiration subtracted. All values are divided by the mean of the no-Collembola control; small grey circles are individual values; large symbols are means; error bars are standard errors of the mean (n= 6 and 30 for brown and black symbols, respectively); violin plots denote the shape of the distributions with the maximum width proportional to the number of observations. Statistics are derived from two-sample t tests with unequal variances.

4 Discussion

We showed that Collembola stimulated CO₂ production from permafrost sediment through “priming” while subtly altering the presence of certain bacterial ASVs, irrespectively of the microbiome manipulation imposed on the Collembola. Colonization of permafrost by Collembola therefore has the potential to result in bacterial community changes and to increase C turnover from tundra soils.

4.1 Introduction of bacteria by Collembola

Collembola introduced new bacteria into thawed permafrost soil, resulting in an increase in alpha diversity and significant differences in bacterial community composition, on a presence–absence basis. This pattern was consistent in all treatments where Collembola were added but did not depend on the microbiome manipulation treatment assigned to the Collembola (Figs. 2b, 3). Despite their effects on the presence of certain bacteria, Collembola did not impose broader community changes as bacterial communities were not affected when considering relative abundances as we had hypothesized (Fig. 2a, Table 1). However, given the subtle effect of Collembola on bacteria during this short incubation, this suggests that on longer timescales Collembola may modify bacterial communities more dramatically.

4.2 Differences between ecto- and endozoochory treatments

Although Collembola introduced new bacteria into permafrost soil, our results did not support that they did so differently through ecto- and/or endozoochory. This suggests that Collembola presence mostly introduced bacteria that were part of the core Collembola microbiome or that the way in which their microbiome changed was not well represented by our treatments. For example, spraying may not have altered the cuticle microbiome as Collembola cuticle is particularly hydrophobic (Hensel et al., 2013a, b; Nickerl et al., 2013). This should however not have been an issue for the endozoochory treatment, in which Collembola were in contact with the topsoil for a prolonged (> 1 month) period of time. Despite this prolonged exposure, no additional effect of the endozoochory treatment on bacterial community composition was observed, compared to ectozoochory or Collembola-only treatments. Collembola typically harbour a gut microbiome reflecting their diet and direct environment (Xiang et al., 2019; Leo et al., 2021), and observing soil-coloured faeces in the topsoil and permafrost stock cultures suggests that soil bacteria were present in the gut of Collembola used in the endozoochory treatments. In contrast with the permafrost stock cultures, where the yeast was necessary as an additional food source to maintain stable populations, Collembola populations remained stable in topsoil stock cultures similar to those used in the endozoochory treatment even in the absence of yeast, further indicating their consumption of soil particles and associated microorganisms. Therefore, the absence of the effect of the ectozoochory treatment on permafrost bacterial communities may indicate that Collembola cuticle does not serve as a vector for bacterial dispersal, possibly due to its omniphobic structure. Further, the lack of difference between the effects of Collembola and endozoochory treatments suggests that few to none of the bacteria present in this topsoil are able to establish in – or survive transit through – collembolan guts. This suggests against generalist “hitch-hiking” bacteria using microfauna guts as a means of dispersal, emphasizing instead that such hitch-hiking bacteria are similar regardless of the feeding context of the Collembola.

It is not unlikely that the observed limited response of the bacterial community to zoochory treatments was transient and that Collembola presence would have eventually affected community composition also in terms of relative abundances. A longer incubation period may have resulted in stronger effects on bacterial community composition, as observed by Coulibaly et al. (2019), but this would have been at the expense of controlling the number of Collembola in the jars and thus of being able to account for their basal respiration. Using an RNA-based approach to target the “active” bacterial community may be another way to detect such transient effects. It should also be noted that our semi-quantitative approach to bacterial community composition did not allow us to assess putative changes in absolute abundances, although such changes are not necessarily observed in similar studies (Kaneda and Kaneko, 2008).

Collembola effects on fungal communities may have indirectly affected bacterial communities and thus represent a potential additional mechanism by which Collembola have affected the presence of certain bacteria. It is not unlikely that fungi would have responded more strongly to the presence of mostly fungivorous Folsomia candida. Using fungi as the focus organism in this study, however, may have had its difficulties. For example, the use of Candida albicans as a food source for the Collembola stock cultures may have strongly distorted fungal community composition, thus requiring questionable bioinformatic workarounds. Further, assessing the effects on fungal communities may not allow for distinguishing between the consequences of Collembola as a vector for microbial migration and as a grazing consumer. Although assessing responses of the whole soil food web may have given better insights into potential mechanisms by which Collembola affect thawed permafrost soil, we here chose to focus on their effects on bacteria given their preponderant role on realizing soil functions in this system (Monteux et al., 2020).

4.3 Collembola effects on CO₂ production

Our hypothesis stating that soil CO₂ production would be increased in the presence of Collembola was supported by our results. However, since no differences in bacterial community composition were detected between the different Collembola treatments, the rationale for our follow-up hypothesis (attributing an increased influence on CO₂ production to the treatments exhibiting stronger changes in bacterial community) was not met. We therefore did not deem it relevant to formally evaluate this hypothesis. Our results suggest that it would be rejected (Fig. A1), but we could not say whether that is because the rationale is wrong or because the manipulation of bacterial communities did not yield the expected effect.

Overall, Collembola presence increased respiration by 25 %, half of which was attributed to respiration of Collembola themselves. We interpret the remaining half (13.22 %) as a priming effect, i.e. a stimulation of soil-organic-matter-derived CO₂ production induced by the presence of Collembola. The mechanisms behind this priming effect remain unclear and could for instance relate to stimulation of microbial activity by the input of more labile substrates by the Collembola (digestion by-products). Overall, our results add to a list of contrasting findings on the direction and magnitude of priming of soil-organic-matter decomposition by Collembola presence in soils in general. Earlier studies that compared CO₂ production in the presence or absence of Collembola or other faunal groups have mostly focused on litter (e.g. Hanlon and Anderson, 1979). The scarce studies on soil CO₂ production including Collembola reported contrasting findings, with some studies observing increases in soil CO₂ emissions in the presence of Collembola by up to +400 % (Addison and Parkinson, 1978; Bakonyi, 1989; Kaneda and Kaneko, 2008; Wang et al., 2017), while others found no significant changes (Theenhaus et al., 1999; Lubbers et al., 2020; Lucas et al., 2020). The Collembola effect on CO₂ production in litter can be species-specific (Hanlon and Anderson, 1979), and our results could support this as we found a similar increase to those in three studies where Folsomia candida was used (Bakonyi, 1989; Kaneda and Kaneko, 2008; Wang et al., 2017). Kaneda and Kaneko (2008) reported a density-dependent increase in CO₂ production with the addition of F. candida, which showed no effect at the density used in our study but an increase with 5–10 times higher density (i.e. 400 individuals in 30 g of soil compared to our 30–80). Accordingly, we did not find any trend supporting a density-dependent response in our data. Nevertheless, we suggest that permafrost soils colonized by Collembola in particular, as well as by soil fauna in general, will likely see an increase in their rates of decomposition, in CO₂ production in particular, due to both priming and the basal respiration of the newly established food web.

Our experimental setup did not allow us to distinguish further between soil and Collembola respiration. We estimated basal respiration, though this may not perfectly reflect reality, considering that Collembola in jars without soil may not have behaved similarly to those in jars with soil. Faced with a scarcity of food sources, they may have slowed their metabolism to endure that stress, which could result in decreased basal respiration and thus lead us to overestimate the magnitude of priming. We could not find reference values in the literature for Folsomia candida respiration rates, but Addison and Parkinson (1978) estimated the CO₂ production of two high-Arctic collembolan species to 182.6 and 250.6 µL CO₂ per gram live weight per hour at 10 ^∘C, for Hypogastrura tullbergi and Folsomia regularis, respectively. Assuming a CO₂ density of 1.977 kg m⁻³ and 170 µg per adult individual of Folsomia candida, our values range between 47 and 220 µL CO₂ per gram live weight per hour. Addison and Parkinson (1978) refer to Mitchell (1973) for a description of the method used for estimating respiration rates, where the measurements appear to be carried out over a period of 20 h. When considering only the first measurement date in our data (3 d after inoculation), the rates are closer to the values reported by Addison and Parkinson (1978) with 157.5±18.5 µL CO₂ per gram live weight per hour (mean ± SE, n= 6). In future studies of this mechanism, it would be interesting to use stable isotope methods to partition the Collembola-derived vs. soil-organic-matter-derived CO₂ production more accurately.

Beyond the CO₂ production increase due to priming by soil fauna, our Collembola addition treatments may have induced “microbial community priming”, in the form of the release of microbial functional limitations. The Yedoma sediment used in this study indeed lacks certain microbial taxa and functions, and the onset of an exotic microbial community can result in large increases in CO₂ production (Monteux et al., 2020). Considering the limited to non-significant (depending on the distance metrics) effect of Collembola addition on bacterial communities, we do not think this mechanism explains our observed results, although we do not rule out this explanation.

5 Conclusions

Collembola presence modified bacterial communities in newly thawed permafrost and primed its CO₂ production, highlighting the importance of non-microbial decomposers for the fate of soil organic matter in thawing permafrost. An emerging theme in permafrost research is the missing functions in permafrost microbial communities due to a lack of certain microbial groups, hampering the production of CO₂, CH₄ and N₂O (Knoblauch et al., 2018; Monteux et al., 2020; Marushchak et al., 2021). It is also becoming clearer that the functionality of permafrost microbes can vary across space (Barbato et al., 2022); therefore the modalities of microbial dispersal into newly thawed permafrost likely affect the fate of the permafrost's organic matter and the rate of release of greenhouse gases. Our findings suggest that Collembola may accelerate the release of greenhouse gases, possibly in part through the introduction of microorganisms, although those did not seem to depend strongly on the preceding Collembola environment. In nature, thawing permafrost mostly occurs at the bottom of the seasonally thawing active layer, often below the water table and at depths far from the topsoil where most Collembola reside. Different soil fauna, such as earthworms, may access such newly thawed permafrost, but their current geographical range does not overlap much with the permafrost region. However, rapid thaw events such as hillslope thermokarst, active-layer detachments or retrogressive thaw slumps expose newly thawed permafrost to surface conditions, thus making it susceptible to colonization by Collembola. Fauna-induced dispersal of microorganisms is most likely to take place in these environments and should be further investigated in realistic field conditions.

It remains unclear whether the increase in CO₂ production we attribute to priming stems from increased functionality of the microbial community due to bridging of functional limitations or from other mechanisms. Mechanical breakdown of organic material and digestion processes might provide microorganisms with more easily available substrates, while disruption of the soil pore structure may also result in increased microbial decomposition activity. Several studies indicate increases in CO₂ production with the presence of soil fauna; however to our best knowledge no studies have specifically partitioned to which extent such increases were due to faunal respiration or to priming effects. We argue that elucidating this question would be an important next step towards opening the “black box” that soil systems still often represent, thus helping to mechanistically address the effects of global changes.

Appendix A

Figure A1 Response ratios of permafrost cumulative CO₂ production rates to addition of Collembola. (a) RR_gross: measured CO₂ production rates. (b) RR_soil: measured CO₂ production rates with estimated Collembola basal respiration subtracted. All values are divided by the mean of the no-Collembola control; small grey circles are individual jar values; large symbols are means; error bars are standard errors of the mean (n= 6). Black symbols denote significant difference from the control; different symbols denote different statistical significance (estimated marginal means treatment contrasts with Holm adjustment for multiple comparisons, Appendix Table A2; dot (⋅) denotes 0.05 < P < 0.1; ^* 0.01 < P < 0.05; ${}^{**}$ 0.001 < P < 0.01; ${}^{***}$ P < 0.001).

Table A1 Effect of Collembola addition on daily CO₂ production over time, repeated-measures ANOVA. P[GG] indicates repeated-measures ANOVA P using Greenhouse–Geisser correction to degrees of freedom to account for violation of the assumption of sphericity (Mauchly's test). P values below 0.05 are denoted in bold; ${}^{***}$ P < 0.001; ^* 0.01 < P < 0.05. Num. DOFs: numerator degrees of freedom; Den. DOFs: denominator degrees of freedom; GES: generalized Eta-squared measure of effect size; GG ϵ: Greenhouse–Geisser ϵ.

Table A2 Effects of Collembola additions on response ratios of cumulative CO₂ production at the end of the incubation, excluding (RR_soil) or including (RR_gross) estimated Collembola respiration. Holm adjustment for P values of non-orthogonal contrasts. P values below 0.05 are denoted in bold; ${}^{***}$ P < 0.001; ${}^{**}$ 0.001 < P < 0.01; ^* 0.01 < P < 0.05; dot (⋅) denotes 0.05 < P < 0.1. DOFs: degrees of freedom; Res. DOFs: residual degrees of freedom.

Appendix B

Figure B1 Number of adult Collembola individuals per incubation jar. Small grey circles are individual jar values; large symbols are means; error bars are standard errors of the mean (n= 6); violin plots denote the shape of the distributions with maximum width proportional to the number of observations. Different black letters denote statistically significant pairwise differences (estimated marginal means pairwise contrasts with Holm adjustment for multiple comparisons, 0.05 < P < 0.1).

Appendix C

Table C1 Effects of Collembola additions on bacterial community composition (permutational multivariate ANOVAs and pairwise contrasts), using distance matrixes based on relative abundances (weighted UniFrac) or presence–absence (non-weighted UniFrac). Holm adjustment for P values of pairwise contrasts. P values below 0.05 are denoted in bold; ${}^{***}$ P < 0.001; ${}^{**}$ 0.001 < P < 0.01; ^* 0.01 < P < 0.05; dot (⋅) denotes 0.05 < P < 0.1.

Table C2 Effects of Collembola additions on bacterial alpha diversity (ANOVAs and pairwise contrasts), with total number of observed ASVs, abundance-based coverage estimator (ACE), and Shannon's and Fisher's diversity indexes. Holm adjustment for P values of pairwise contrasts. P values below 0.05 are denoted in bold; ${}^{***}$ P < 0.001; ${}^{**}$ 0.001 < P < 0.01; ^* 0.01 < P < 0.05; dot (⋅) denotes 0.05 < P < 0.1.