Responses of soil hexapod communities to warming are mediated by microbial carbon and nitrogen in a subarctic grassland

Warming in subarctic ecosystems will be two-fold higher compared to lower latitudes under current climate change projections. While the effects of warming in northern ecosystems on plants and microorganisms have been extensively studied, the responses of soil fauna have received much less attention, despite their important role in regulating key soil processes. We analyzed the response of soil hexapod communities in a subarctic grassland exposed to a natural geothermal gradient in Iceland with increases of + 3 and + 6 ◦ C above ambient temperature. We characterized hexapod communities using environmental DNA (eDNA) metabarcoding. We analyzed the amounts of microbial carbon (Cmic), microbial N (Nmic), dissolved organic C (DOC) and dissolved organic N (DON) and then assessed whether these variables could help to account for the compositional dissimilarity of ground hexapod communities across temperatures. The increases in soil temperature did lead to changes in the composition of hexapod communities. The compositional differences caused by + 6 ◦ C plots were correlated with a decrease in Cmic and Nmic, soil DOC and DON. Our results highlight the response of soil hexapods to warming, and their interaction with microbial biomass ultimately correlated with changes in the availabilities of soil C and N.


Introduction
Northern regions are experiencing the fastest warming rate on Earth [1,2].Under these new temperature conditions soil microbial activity is accelerated, increasing decomposition rates and reducing the stocks of labile carbon (C; Saad & Conrad, 1993).Despite the initial increase in microbial activity, substrate depletion leads to the loss of microbial C (Cmic), as well as microbial N (Nmic; Marañón-Jiménez et al., 2019; [3].The ultimate depletion of labile C and N may have a larger impact on plant productivity than the direct effects of temperature alone, and so the impacts on nutrient dynamics are important for ecosystem functioning in northern latitudes [4][5][6].However, the impacts of C and N availability and warming on soil biodiversity of subarctic regions are mostly unknown, despite subarctic ecosystems accounting for the largest pool of soil C on the surface of Earth [7]. Most information about the response of subarctic species to climate change is limited to the aboveground ecosystems, with soil fauna being mostly unexplored [8].Previous studies have not reached a consensus about how subarctic hexapod diversity responds to warming.While some studies have argued that the direct reaction should be strong, with generalized reductions of hexapod species abundances and drastic alterations in community composition [9], other studies have suggested that these changes may only be transient due to their high resilience at a community level [10,11].Accordingly, some studies have reported that soil warming can even improve the fitness of particular hexapod species, which could even expand their areas of distribution [12,13].However, global warming can also have multiple indirect effects over soil biodiversity [14].Thus, downstream shifts in C and N availability as a result of the effect of soil warming over soil microbial communities and vegetation can strongly impact indirectly the communities of soil arthropods [15].Nonetheless, the impact of C and N throughout soil food-webs remains poorly studied.Most soil ecosystems are relatively poor in N, and many detritivore hexapods have low-N diets and Nmic is their main pathway for N assimilation [16].Therefore, the abundances of some hexapod species may decrease with lower C to N ratio (C:N), while other microvorous arthropods may exploit the increases in bacterial and fungal biomass (Peguero et al., 2022).This complexity and the lack of a general pattern linking warming with C and N availability, microbial decomposers and hexapod communities in soil hampers our ability to predict their combined responses to climate change and the consequences for ecosystem functioning, which may be particularly sensitive at higher latitudes.
The role of soil hexapods in C and N cycling is still poorly known despite their prominent contribution to the dynamics of soil nutrients [17].Hexapod detritivores break down soil organic matter into particles with higher surface to volume ratios, facilitating their further decomposition and mineralization by the local microbial community [18,19].Microbial feeders, though, exert a strong top-down control on microbial communities, thus regulating their abundance, diversity and activity [19,20].The typically cryptic morphological traits used in the identification of soil hexapod species is likely one of the main reasons hindering our knowledge of these communities.Advanced molecular techniques like metabarcoding have recently facilitated their identification, even using DNA remnants in the soil [21,22].This environmental DNA (eDNA) can be easily extracted from soil samples and allows a detailed community description and even the discovery of species that otherwise could not be detected [23], thereby facilitating the study of soil faunal communities.
Here, we studied the impact of soil warming on soil hexapod community composition in a subarctic grassland.We collected soil samples across a natural geothermal site in Iceland at ambient temperature, +3 and + 6 • C above ambient.We characterized soil hexapod communities using eDNA metabarcoding and analyzed the amounts of Cmic, Nmic, soil dissolved organic C (DOC) and soil dissolved organic N (DON).We then assessed whether these variables could help to account for the compositional dissimilarity of hexapod communities across temperatures.We hypothesized that soil warming would have an impact on the structure and composition of soil hexapod communities likely through indirect effects arising from alterations of microbial communities' C and N concentrations.

Site description and experimental conditions
This study was conducted at the ForHot research site in Iceland [24] between August 2017 and June 2018 (64 • 0 ′ N, 21 • 11 ′ W).Soil type was a Brown Andosol [25].Mean annual temperature at the site was 5.1 • C. The coldest and warmest temperatures in the neighboring village of Eyrarbakki in 2016 were − 12.3 • C and 21.6 • C, respectively.Average annual precipitation for the same year was 1153 mm [26].The vegetation was an unmanaged grassland dominated by Agrostis capillaris L., Galium boreale L. and Anthoxantum odoratum L.. Vascular plants cover 46% of the area over a moss mat which covers up to 88% of the ground.Natural N deposition in the area is 1.3 ± 0.1 kg N ha-1 y-1 [27].This grassland has been geothermally warmed since 29 May 2008, when an earthquake transferred geothermal energy from hot groundwater to previously unheated soils [24].Belowground temperatures at 10 cm depth now display a permanent warming gradient reaching +10 • C, with a less severe increase in temperature at the soil surface of +0.2 C. The warming has only been mildly disruptive respective to previous conditions, and the warmed area experiences a similar magnitude of warming and cooling over the seasons.Soil humidity was only marginally affected, with volumetric water content changing from 40% to 38%, and water pH increased from 5.6 in unheated soil to up to 6.3 after warming.Geothermal groundwater has remained in the bedrock and has not reached the root zone, thus avoiding direct eco-toxicological effects [24].The resulting stable conditions and lack of artifacts provide a realistic natural belowground experiment on soil warming under climate change.Even though, we acknowledge that aboveground effects of the on-going warming should lead to an opposite warming gradient across the soil profile (i.e., warmer at the top) and it may also entail complex plant-soil feedbacks that may not be entirely captured by the geothermal warming present at our study site.Five transects were established, each one consisting of three 2 × 2 m plots, and each plot at different temperature: an unheated control, a low warming level of ca.+3 • C and a higher warming level of ca.+6 • C above the ambient reference in the control (henceforth referred as "+3 • C" and "+6 • C").

Soil core sampling
Soil cores were collected using an auger to a depth of ~10 cm, excluding the O horizon.Soil cores were sampled seasonally four times: August 2017, corresponding to late growing season; November 2017, at start of winter and initial soil freezing; April 2018, with the first soil thaw in un-warmed soils, and June 2018, in the early part of the growing season.We thus collected a total of 20 core samples for each warming treatment (5 replicates in 4 seasons for 3 temperature levels = samples).All samples were immediately sieved to remove roots and stones larger than 2 mm.Fifteen grams of each sample were then frozen in plastic bags in liquid N in the field to immediately stop all biological processes.All frozen samples were then freeze-dried in the laboratory following the standard protocol of the commercial lyophilizer for this type of sample (FreeZone 2.5 L -50C Benchtop Freeze Dryers, LabConco Corp., Kansas City, MO.USA).eDNA was extracted from 15 g soil samples as previously described [22,28].

. Metabarcoding analysis
The soil hexapod communities were characterized based on Molecular Operational Taxonomic Units (MOTUs) using an eDNA metabarcoding approach.We amplified the 16S mitochondrial rDNA region using the Ins16S_l primer pair (Ins16S_1-F: 5 ′ -TRRGACGAGAA-GACCCTATA-3 ′ ; Ins16_1-R: 5 ′ -TCTTAATCCAACATCGAGGTC-3 ′ ; Clarke et al. 2014).This primer pair, specifically designed for hexapod metabarcoding, introduces a very limited taxonomic bias and performs very well for identifications at the species level throughout the Hexapoda subphylum (e.g.Refs.[29,30].PCR amplification was performed in triplicate in 20-μL mixtures consisting of 10 μL of AmpliTaq Gold Master Mix (Life Technologies, Carlsbad, USA), 5.84 μL of nuclease-free Ambion water (Thermo Fisher Scientific, Waltham, USA), 0.25 μM each primer, 3.2 μg of bovine serum albumin (Roche Diagnostic, Basel, Switzerland) and 2 μl of DNA template that was diluted 10-fold to reduce PCR inhibition by humic substances.The thermal profile of the PCR amplification was 40 cycles of denaturation at 95 • C (30 s), annealing at 49 • C (30 s) and elongation at 72 • C (60 s), with a final elongation step at 72 • C for min.Tags had at least five differences between them to minimize ambiguities [31].The sequenced multiplexes comprised extractions/PCR blank controls, unused tag combinations and positive controls [29].The PCR products were then sequenced using the MiSeq platform (Illumina Inc., San Diego, USA), with the expected sequencing depth set at 000 reads per sample.
The sequences were processed using ObiTools software [32].Low-quality sequences (containing Ns, alignment scores <50, lengths <140 bp or >320 bp and singletons) were excluded.The remaining sequences were clustered into MOTUs using SUMACLUST [33] at a threshold of sequence similarity of 97%.The final number of MOTUs after curation was 11785.The hexapod MOTUs were taxonomically assigned using Basic Local Alignment Search Tool, (BLAST), with a query coverage criterion of 95%.MOTUs showing <80% similarity with either the local or the EMBL reference databases were removed, leading to 590 MOTUs [34].These retained MOTUs included taxa from classes Insecta and Entognatha, which both belong to the subphylum Hexapoda.All sequences with a frequency of occurrence below 0.05 per MOTU and per run were discarded.This threshold was empirically determined to clear the negative sequencing controls (non-used tag combinations) included in our global data production procedure.We then applied a post-processing pipeline [35] to minimize PCR and sequencing errors, contaminations and false-positive sequences, and a detailed curation of ecologically incongruent assignments also to deal with the ambiguous matches and provide more reliable information about species identifications (e.g., taxa with distributions outside the Palearctic and Nearctic zones, none barcoded MOTUs and redundant assignations).This conservative approach retained a total of 40 identified species.We then used checklists of Icelandic hexapod species and information from previous studies at the same study site [10,36] to assess the performance of our eDNA metabarcoding protocol to properly describe the hexapod communities in the soil.

. Environmental variables
Dissolved organic C (DOC), dissolved organic N and total dissolved N (TDN) in all soils were quantified in 1 M KCl extracts.Soil Cmic and Nmic were determined by the chloroform-fumigation extraction method (48 h incubation period), followed by 1 M KCl extraction.All extracts were analyzed for DOC and TDN with a TC/TN-Analyzer (Shimadzu, TOC-VCPH/CPNTNM-1 analyzer).There were no correlations between environmental variables (Tables S1-S4).

Data analysis
All data handling, visualization and statistical analyses were conducted using R v4.0.6 [37].We first identified the relationships between all environmental variables (the amounts of DOC and DON, and Cmic and Nmic) across temperatures using principal component analyses (PCAs) and simple general linear models (GLMs), with warming as categorical fixed-effect terms.We then built simple GLMs, with temperature as fixed-effects terms, to identify differences in rarefied mOTU abundance and richness due to warming.Changes in the composition of the hexapod communities for each pair of temperature levels were assessed using sparse partial least squares discriminant analysis (sPLS-DA) as implemented in the mixOmics package [38].sPLS-DA selects variables and classifies the most discriminative taxa in the community matrices.The optimal number of components was chosen based on the error rate of t-tests, which suggested the use of single-component sPLS-DAs in both warming levels.The mixOmics package also delivers P values based on the area under the receiver operating characteristic curve to complement the sPLS-DA performance results.The main sPLS-DA variates (the x-variate, capturing the compositional similarity of soil hexapod communities) were then subjected to simple GLMs against the set of environmental variables: the amounts of soil DOC and DON, Cmic and Nmic.

Environmental variables
The amounts of Cmic, soil DOC and DON decreased with warming, while Nmic did so marginally (P < 0.05, P < 0.05, P < 0.05 & P = 0.08, respectively; Table 1).The confidence ellipses of the environmental PCAs overlap between the control and the temperature levels (Fig. 1).

Community response to soil warming
The number of eDNA reads did not differ significantly across temperature levels (P = 0.54; Fig. 2a), nor did the richness of the soil hexapod communities (P = 0.29; Fig. 2b).The sPLS-DA with the warming levels identified significant compositional dissimilarities between the control, the +3 and + 6 • C treatments (P < 0.01 for both, Fig. 3a and b).The first sPLS-DA variate (x-variate) from the control and the +3 • C treatment accounted for 11% of the variance, the variate including the control and the +6 • C treatment accounted for 19%, and the variate including both warming levels explained 16% of variance (Fig. 3).In all sPLS-DAs, higher scores in the x-variate implied increasing compositional dissimilarity with higher temperatures (Fig. 3).Compositional changes with warming were driven mainly by springtails (e.g.Tomocerus ocreatus, Folsomia quadrioculata and Megalothorax svalbardensis), beetles (e.g.Liogluta microptera and Hypnoidus sp.1), and the dipteran Bradysia subvernalis (Fig. 4).
None of the four environmental variables had a significant effect at +3 • C when we assessed the capacity of the environmental variables to account for such differences in the composition of the soil hexapod community.Compositional similarities at +6 • C, however, were negatively correlated with the amounts of Cmic and Nmic and marginally correlated with the amount of DOC and DON found in the control (P < 0.05, P < 0.05, P = 0.07 & P = 0.06, respectively; Table 2).The negative correlation between the Control vs +6 • C sPLS-DA x-variate and these environmental variables imply that Cmic, Nmic, soil DOC and DON depletion lead to communities with species compositions similar to those found at the +6 • C (Fig. 3).compositions of the soil hexapod community in a subarctic grassland.In this case of study, an increase in soil temperature of 6 • C above ambient had an impact on the community composition, which correlated with decreases in the amounts of Cmic, Nmic, DOC and DON.These results highlight how the responses of soil hexapod community assemblages to warming can be associated with changes in microbial biomass and shifts in C and N availability.Moreover, eDNA metabarcoding was an efficient methodological approach to obtain a description of the soil hexapod communities.Additionally, eDNA metabarcoding allowed us to detect three springtail species not previously reported for Iceland but extensively distributed across the northernmost part of continental Europe (GBIF, 2021).

Our results indicated that warming causes changes in the
Our results support previous findings that exposure to soil warming for over a decade in this subarctic grassland may not result in a shift of either soil hexapod abundance or richness [10].We found, however, compositional dissimilarities that may threaten the functionality of soil hexapod communities more than can simple changes in species richness or abundance [18,39].Shifts in community composition in our study were mainly driven by springtail, beetle and fly species, which were also the richest lineages with the most DNA reads.The species driving such changes in community composition differed not only between the control and higher temperature levels, but also between the +3 and + 6 • C treatments.These differences may have been associated with different temperature optima, sensitivity to temperature maxima or resilience to warming by the populations of hexapod species [12,13].For instance, the rove beetle Liogluta microptera and wireworms of the genus Hypnoidus are known to prefer warmer and drier soils (Ottesen 1996; [40], while some fungus gnats of the genus Bradysia are also known to have a temperature-sensitive phenology [41] (Fig. 4).Changes in food resources caused by warming are also an important factor capable of altering soil hexapod communities [10,42], both by affecting the amount of litter input or the microbial biomass (S.F [43]. The impacts on soil hexapod communities under the most extreme conditions of soil warming correlated with the depletion of Cmic and Nmic in the topsoil.Previous studies at the same site have found that warming leads to the loss of soil organic C (SOC) stocks in the topsoil ultimately reducing the amount of Cmic [3,44,45].A proportional loss of Nmic might have followed this loss of C due to the relatively tight C:N stoichiometric balance of microbial metabolism at the site [46][47][48].The depletion of labile DOC and DON in the warmest soils in our study, however, only affected marginally the soil hexapod communities (Table 2, Fig. 3).The reduction of Cmic and DOC should have negatively affected important soil hexapod trophic guilds such as microbial feeders and detritivores relying on bacterial and fungal biomass as well as litter input.In contrast, the lack of effect of the moderate +3 • C warming may have been associated with the threshold dynamics of the soil   Confidence ellipses are set at 95%.The linear relationships between the x-variate scores with the corresponding soil environmental variables (microbial carbon, microbial nitrogen, soil dissolved organic carbon and nitrogen -DOC and DOC, respectively-) are shown on top of each sPLS-DA.Solid and dashed lines show significant and marginally significant slope parameters (P < 0.05, and P < 0.1, respectively).See Table 2 for further information on the linear models.microorganisms found in this subarctic grassland, where only increases in temperature equal or above 6 • C have been shown to lead to abrupt shifts in microbial diversity and composition [49,50].
Global warming is generally expected to increase SOM decomposition leading to higher N availability [51,52].However, at a decadal scale the impact of warming correlated with lower DOC and DON availability as well as with a loss of Cmic and Nmic.The impact of warming was therefore quite different from the predicted short-term warming response, which should be driven by the initial enhanced availability of N. The requirements of hexapods for C and N, and the sources of C and N, are of great interest for predicting the response of soil hexapods to warming.In our study, we consider it likely that the loss of Cmic and Nmic due to warming depleted the food resources of hexapod bacterivores and fungivores [53], increasing litter concentration and exerting a positive feedback on hexapod detritivores.These differences in C and N availability, together with the species-specific temperature optima and sensitivity, may account for the lack of a common compositional response across warming levels.The stoichiometric balance of microorganisms and hexapods also highlight the importance of the C and N to understand the mechanistic response of soil hexapods to warming.The nutrient requirements and elemental stoichiometry of microorganisms and plants have been well studied [54], but our knowledge of hexapod stoichiometry is a lot less advanced [55].For instance, previous work suggests that arthropod stoichiometric response to temperature is not necessarily monotonic in shape [56], and that resource quality (e.g.C:N and C:P ratios) plays a major role in ground arthropod feeding behavior [57,58].Therefore, the complex behavior of arthropods hamper our ability to assess their nutritional requirements [59,60].
Furthermore, it is crucial to consider potential biases in these responses, as the study system differs from global warming in two significant aspects.Firstly, geothermal warming exhibits a contrasting vertical gradient compared to climatic warming, with warmer temperatures increasing from deep soil layers rather than the surface.Previous research indicates that such warming influences the composition of soil hexapod communities across the soil profile.For instance, in regions with a natural gradient of European climatic aridity, specific traits of soil-dwelling springtail species are favored on the soil surface, whereas short-term exposure to drought induces a different response, promoting traits associated with soil conditions (Ferrín et al., 2021).Secondly, global warming can lead to shifts in vegetation communities, which also have implications for hexapod populations.In subarctic regions, one common response to warming is the expansion of shrubs into grasslands and tundra, a process known as shrubification.This vegetation change alters nutrient dynamics, microbial activity, and hexapod communities [61][62][63].Therefore, it is essential to interpret these results with caution, Fig. 4. Identity of the main hexapod species driving the compositional differences between the controls and the warming treatments.Only species with a loading score >0.1 in the corresponding discriminant analysis are shown.

Table 2
Relationship between the compositional similarities of the soil hexapod community with the environmental variables across temperatures.Results are general linear models, with the first variate (X-variate 1) from an sPLS-DA of the hexapod communities as the response variable modeled against each environmental variable: microbial carbon (C), microbial nitrogen (N), soil dissolved organic carbon (DOC) or soil dissolved organic nitrogen (DON).Effect estimates are followed by their standard errors.•, P < 0.1; *, P < 0.05.All models have 8 • of freedom.Note that the higher the value of the sPLS-DA xvariate, the higher the similarity with the communities of the higher temperature level of the pair (see Fig. 3).
considering the unique characteristics of the study system and the potential influence of these factors.Nonetheless, experimental and mechanistic approaches allow the identification of the nutritional demands affecting soil hexapods over a multiplicity of biomes, and by extension their community responses to warming, clarifying belowground ecological dynamics.Future studies on this topic are therefore strongly recommended.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Fig. 1 .
Fig. 1.Principal component analyses of the distribution of the environmental variables, with 95% confidence ellipses denoting different treatments.

Fig. 2 .
Fig. 2. Violin plots of the numbers of reads (a) and molecular operative taxonomic units -mOTU-(b) across sites and treatments.Horizontal lines inside each plot from the lowest to the highest denote the first, second and third quartiles.

Fig. 3 .
Fig. 3. Compositional variation of the soil hexapod communities based on a single-component sparse partial least squares discriminant analysis (sPLS-DA) for control against +3 • C warming (a), control against +6 • C warming (b) and +3 • C against +6 • C (c).The amount of explained variance of the x-variate is shown in each figure.Confidence ellipses are set at 95%.The linear relationships between the x-variate scores with the corresponding soil environmental variables (microbial carbon, microbial nitrogen, soil dissolved organic carbon and nitrogen -DOC and DOC, respectively-) are shown on top of each sPLS-DA.Solid and dashed lines show significant and marginally significant slope parameters (P < 0.05, and P < 0.1, respectively).See Table2for further information on the linear models.

Table 1
Variation of environmental variables across temperatures.
Effect estimates followed by their standard errors are the output of general linear models for each environmental variable.The intercept is always the control, and the explanatory variable is always the corresponding environmental variable: microbial carbon (C), microbial nitrogen (N), soil dissolved organic carbon (DOC) or soil dissolved organic nitrogen (DON).•,P < 0.1; *, P < 0.05.All models have 13 • of freedom.M. Ferrín et al.