European Journal of Soil Biology 117 (2023) 103513 Contents lists available at ScienceDirect European Journal of Soil Biology journal homepage: www.elsevier.com/locate/ejsobi Original article Responses of soil hexapod communities to warming are mediated by microbial carbon and nitrogen in a subarctic grassland Miquel Ferrín a, b, *, Josep Peñuelas a, b, 1, Albert Gargallo-Garriga a, b, c, d, Amaia Iribar e, Ivan A. Janssens f, Sara Marañón-Jiménez b, d, Jérôme Murienne e, Andreas Richter g, Bjarni D. Sigurdsson h, Guille Peguero a, i, 1 a CSIC, Global Ecology Unit CREAF-CSIC-UAB, 08913, Bellaterra, Catalonia, Spain b CREAF, 08913, Cerdanyola del Vallès, Catalonia, Spain c Global Change Research Institute, The Czech Academy of Sciences, Belidla 986/4a, CZ-60300, Brno, Czech Republic d Universitat Autònoma de Barcelona, 08193, Bellaterra, Spain e Laboratoire Evolution et Diversité Biologique (UMR5174 EDB), Université Toulouse 3 Paul Sabatier, CNRS, IRD, 118 route de Narbonne, Toulouse, France f Department of Biology, Antwerp University, Universiteitsplein 1, B-2610, Wilrijk, Belgium g Centre for Microbiology and Environmental Systems Science, University of Vienna, Djarssiplatz 1, 1030, Wien, Austria h Agricultural University of Iceland, 112 Keldnaholt, Reykjavik, Iceland i Department of Evolutionary Biology, Ecology and Environmental Sciences, Universitat de Barcelona, 08028, Barcelona, Spain A R T I C L E I N F O A B S T R A C T Handling Editor: B. Griffiths 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 Keywords: been extensively studied, the responses of soil fauna have received much less attention, despite their important Climate change role in regulating key soil processes. We analyzed the response of soil hexapod communities in a subarctic Environmental DNA grassland exposed to a natural geothermal gradient in Iceland with increases of +3 and + 6 ◦ C above ambient Metabarcoding temperature. We characterized hexapod communities using environmental DNA (eDNA) metabarcoding. We Microbial communities Community composition analyzed the amounts of microbial carbon (Cmic), microbial N (Nmic), dissolved organic C (DOC) and dissolved DOC organic N (DON) and then assessed whether these variables could help to account for the compositional DON 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. 1. Introduction the impacts on nutrient dynamics are important for ecosystem func­ tioning in northern latitudes [4–6]. However, the impacts of C and N Northern regions are experiencing the fastest warming rate on Earth availability and warming on soil biodiversity of subarctic regions are [1,2]. Under these new temperature conditions soil microbial activity is mostly unknown, despite subarctic ecosystems accounting for the largest accelerated, increasing decomposition rates and reducing the stocks of pool of soil C on the surface of Earth [7]. labile carbon (C; Saad & Conrad, 1993). Despite the initial increase in Most information about the response of subarctic species to climate microbial activity, substrate depletion leads to the loss of microbial C change is limited to the aboveground ecosystems, with soil fauna being (Cmic), as well as microbial N (Nmic; Marañón-Jiménez et al., 2019; [3]. mostly unexplored [8]. Previous studies have not reached a consensus The ultimate depletion of labile C and N may have a larger impact on about how subarctic hexapod diversity responds to warming. While plant productivity than the direct effects of temperature alone, and so some studies have argued that the direct reaction should be strong, with * Corresponding author. CSIC, Global Ecology Unit CREAF-CSIC-UAB, 08913, Bellaterra, Catalonia, Spain. E-mail address: miquelferrin@gmail.com (M. Ferrín). 1 These authors share senior authorship. https://doi.org/10.1016/j.ejsobi.2023.103513 Received 13 April 2023; Received in revised form 29 May 2023; Accepted 5 June 2023 Available online 11 June 2023 1164-5563/© 2023 The Authors. Published by Elsevier Masson SAS. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 generalized reductions of hexapod species abundances and drastic al­ earthquake transferred geothermal energy from hot groundwater to terations in community composition [9], other studies have suggested previously unheated soils [24]. Belowground temperatures at 10 cm that these changes may only be transient due to their high resilience at a depth now display a permanent warming gradient reaching +10 ◦ C, community level [10,11]. Accordingly, some studies have reported that with a less severe increase in temperature at the soil surface of +0.2 C. soil warming can even improve the fitness of particular hexapod species, The warming has only been mildly disruptive respective to previous which could even expand their areas of distribution [12,13]. However, conditions, and the warmed area experiences a similar magnitude of global warming can also have multiple indirect effects over soil biodi­ warming and cooling over the seasons. Soil humidity was only versity [14]. Thus, downstream shifts in C and N availability as a result marginally affected, with volumetric water content changing from 40% of the effect of soil warming over soil microbial communities and to 38%, and water pH increased from 5.6 in unheated soil to up to 6.3 vegetation can strongly impact indirectly the communities of soil ar­ after warming. Geothermal groundwater has remained in the bedrock thropods [15]. Nonetheless, the impact of C and N throughout soil and has not reached the root zone, thus avoiding direct eco-toxicological food-webs remains poorly studied. Most soil ecosystems are relatively effects [24]. The resulting stable conditions and lack of artifacts provide poor in N, and many detritivore hexapods have low-N diets and Nmic is a realistic natural belowground experiment on soil warming under their main pathway for N assimilation [16]. Therefore, the abundances climate change. Even though, we acknowledge that aboveground effects of some hexapod species may decrease with lower C to N ratio (C:N), of the on-going warming should lead to an opposite warming gradient while other microvorous arthropods may exploit the increases in bac­ across the soil profile (i.e., warmer at the top) and it may also entail terial and fungal biomass (Peguero et al., 2022). This complexity and the complex plant-soil feedbacks that may not be entirely captured by the lack of a general pattern linking warming with C and N availability, geothermal warming present at our study site. Five transects were microbial decomposers and hexapod communities in soil hampers our established, each one consisting of three 2 × 2 m plots, and each plot at ability to predict their combined responses to climate change and the different temperature: an unheated control, a low warming level of ca. consequences for ecosystem functioning, which may be particularly +3 ◦ C and a higher warming level of ca. +6 ◦ C above the ambient sensitive at higher latitudes. reference in the control (henceforth referred as “+3 ◦ C” and “+6 ◦ C”). The role of soil hexapods in C and N cycling is still poorly known despite their prominent contribution to the dynamics of soil nutrients 2.2. Soil core sampling [17]. Hexapod detritivores break down soil organic matter into particles with higher surface to volume ratios, facilitating their further decom­ Soil cores were collected using an auger to a depth of ~10 cm, position and mineralization by the local microbial community [18,19]. excluding the O horizon. Soil cores were sampled seasonally four times: Microbial feeders, though, exert a strong top-down control on microbial August 2017, corresponding to late growing season; November 2017, at communities, thus regulating their abundance, diversity and activity start of winter and initial soil freezing; April 2018, with the first soil [19,20]. The typically cryptic morphological traits used in the identifi­ thaw in un-warmed soils, and June 2018, in the early part of the growing cation of soil hexapod species is likely one of the main reasons hindering season. We thus collected a total of 20 core samples for each warming our knowledge of these communities. Advanced molecular techniques treatment (5 replicates in 4 seasons for 3 temperature levels = 60 like metabarcoding have recently facilitated their identification, even samples). All samples were immediately sieved to remove roots and using DNA remnants in the soil [21,22]. This environmental DNA stones larger than 2 mm. Fifteen grams of each sample were then frozen (eDNA) can be easily extracted from soil samples and allows a detailed in plastic bags in liquid N in the field to immediately stop all biological community description and even the discovery of species that otherwise processes. All frozen samples were then freeze-dried in the laboratory could not be detected [23], thereby facilitating the study of soil faunal following the standard protocol of the commercial lyophilizer for this communities. type of sample (FreeZone 2.5 L -50C Benchtop Freeze Dryers, LabConco Here, we studied the impact of soil warming on soil hexapod com­ Corp., Kansas City, MO. USA). eDNA was extracted from 15 g soil munity composition in a subarctic grassland. We collected soil samples samples as previously described [22,28]. across a natural geothermal site in Iceland at ambient temperature, +3 and + 6 ◦ C above ambient. We characterized soil hexapod communities 2.3. . Metabarcoding analysis using eDNA metabarcoding and analyzed the amounts of Cmic, Nmic, soil dissolved organic C (DOC) and soil dissolved organic N (DON). We The soil hexapod communities were characterized based on Molec­ then assessed whether these variables could help to account for the ular Operational Taxonomic Units (MOTUs) using an eDNA meta­ compositional dissimilarity of hexapod communities across tempera­ barcoding approach. We amplified the 16S mitochondrial rDNA region tures. We hypothesized that soil warming would have an impact on the using the Ins16S_l primer pair (Ins16S_1-F: 5′ -TRRGACGAGAA­ structure and composition of soil hexapod communities likely through GACCCTATA-3′ ; Ins16_1-R: 5′ -TCTTAATCCAACATCGAGGTC-3′ ; Clarke indirect effects arising from alterations of microbial communities’ C and et al. 2014). This primer pair, specifically designed for hexapod meta­ N concentrations. barcoding, introduces a very limited taxonomic bias and performs very well for identifications at the species level throughout the Hexapoda 2. Material and methods subphylum (e.g. Refs. [29,30]. PCR amplification was performed in triplicate in 20-μL mixtures consisting of 10 μL of AmpliTaq Gold Master 2.1. Site description and experimental conditions Mix (Life Technologies, Carlsbad, USA), 5.84 μL of nuclease-free Ambion water (Thermo Fisher Scientific, Waltham, USA), 0.25 μM each primer, This study was conducted at the ForHot research site in Iceland [24] 3.2 μg of bovine serum albumin (Roche Diagnostic, Basel, Switzerland) between August 2017 and June 2018 (64◦ 0′ N, 21◦ 11′ W). Soil type was a and 2 μl of DNA template that was diluted 10-fold to reduce PCR inhi­ Brown Andosol [25]. Mean annual temperature at the site was 5.1 ◦ C. bition by humic substances. The thermal profile of the PCR amplification The coldest and warmest temperatures in the neighboring village of was 40 cycles of denaturation at 95 ◦ C (30 s), annealing at 49 ◦ C (30 s) Eyrarbakki in 2016 were − 12.3 ◦ C and 21.6 ◦ C, respectively. Average and elongation at 72 ◦ C (60 s), with a final elongation step at 72 ◦ C for 7 annual precipitation for the same year was 1153 mm [26]. The vege­ min. Tags had at least five differences between them to minimize am­ tation was an unmanaged grassland dominated by Agrostis capillaris L., biguities [31]. The sequenced multiplexes comprised extractions/PCR Galium boreale L. and Anthoxantum odoratum L.. Vascular plants cover blank controls, unused tag combinations and positive controls [29]. The 46% of the area over a moss mat which covers up to 88% of the ground. PCR products were then sequenced using the MiSeq platform (Illumina Natural N deposition in the area is 1.3 ± 0.1 kg N ha-1 y-1 [27]. This Inc., San Diego, USA), with the expected sequencing depth set at 400 grassland has been geothermally warmed since 29 May 2008, when an 000 reads per sample. 2 M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 The sequences were processed using ObiTools software [32]. while Nmic did so marginally (P < 0.05, P < 0.05, P < 0.05 & P = 0.08, Low-quality sequences (containing Ns, alignment scores <50, lengths respectively; Table 1). The confidence ellipses of the environmental <140 bp or >320 bp and singletons) were excluded. The remaining PCAs overlap between the control and the temperature levels (Fig. 1). sequences were clustered into MOTUs using SUMACLUST [33] at a threshold of sequence similarity of 97%. The final number of MOTUs 3.2. Insect community description after curation was 11785. The hexapod MOTUs were taxonomically assigned using Basic Local Alignment Search Tool, (BLAST), with a eDNA metabarcoding identified 40 species, corresponding to query coverage criterion of 95%. MOTUs showing <80% similarity with different lineages within the subphylum Hexapoda: Collembola (14 either the local or the EMBL reference databases were removed, leading species), Coleoptera (12 species), Diptera (six species), Hemiptera (five to 590 MOTUs [34]. These retained MOTUs included taxa from classes species), Hymenoptera (one species), Dermaptera (one species) and Insecta and Entognatha, which both belong to the subphylum Hexapoda. Lepidoptera (one species). Three springtail species were previously un­ All sequences with a frequency of occurrence below 0.05 per MOTU and recorded for Iceland: Sminthurinus bimaculatus (244 reads across nine per run were discarded. This threshold was empirically determined to samples), Pogonognathellus flavescens (110 reads across three samples) clear the negative sequencing controls (non-used tag combinations) and Megalothorax perspicillum (two reads in one sample). The most included in our global data production procedure. We then applied a abundant group was four species of springtails belonging to the order post-processing pipeline [35] to minimize PCR and sequencing errors, Neelipleona, which together accounted for 42% of all reads. The next contaminations and false-positive sequences, and a detailed curation of most abundant orders were Diptera (18.6%), Entomobryomorpha ecologically incongruent assignments also to deal with the ambiguous (9.1%), Hemiptera (9.1%), Symphypleona (8.4%) and Poduromorpha matches and provide more reliable information about species identifi­ (7.6%). cations (e.g., taxa with distributions outside the Palearctic and Nearctic zones, none barcoded MOTUs and redundant assignations). This con­ 3.3. Community response to soil warming servative approach retained a total of 40 identified species. We then used checklists of Icelandic hexapod species and information from pre­ The number of eDNA reads did not differ significantly across tem­ vious studies at the same study site [10,36] to assess the performance of perature levels (P = 0.54; Fig. 2a), nor did the richness of the soil our eDNA metabarcoding protocol to properly describe the hexapod hexapod communities (P = 0.29; Fig. 2b). The sPLS-DA with the communities in the soil. warming levels identified significant compositional dissimilarities be­ tween the control, the +3 and + 6 ◦ C treatments (P < 0.01 for both, 2.4. . Environmental variables 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 Dissolved organic C (DOC), dissolved organic N and total dissolved N including the control and the +6 ◦ C treatment accounted for 19%, and (TDN) in all soils were quantified in 1 M KCl extracts. Soil Cmic and the variate including both warming levels explained 16% of variance Nmic were determined by the chloroform-fumigation extraction method (Fig. 3). In all sPLS-DAs, higher scores in the x-variate implied increasing (48 h incubation period), followed by 1 M KCl extraction. All extracts compositional dissimilarity with higher temperatures (Fig. 3). Compo­ were analyzed for DOC and TDN with a TC/TN- Analyzer (Shimadzu, sitional changes with warming were driven mainly by springtails (e.g. TOC-VCPH/CPNTNM-1 analyzer). There were no correlations between Tomocerus ocreatus, Folsomia quadrioculata and Megalothorax svalbar­ environmental variables (Tables S1–S4). densis), beetles (e.g. Liogluta microptera and Hypnoidus sp. 1), and the dipteran Bradysia subvernalis (Fig. 4). 2.5. Data analysis None of the four environmental variables had a significant effect at +3 ◦ C when we assessed the capacity of the environmental variables to All data handling, visualization and statistical analyses were con­ account for such differences in the composition of the soil hexapod ducted using R v4.0.6 [37]. We first identified the relationships between community. Compositional similarities at +6 ◦ C, however, were nega­ all environmental variables (the amounts of DOC and DON, and Cmic tively correlated with the amounts of Cmic and Nmic and marginally and Nmic) across temperatures using principal component analyses correlated with the amount of DOC and DON found in the control (P < (PCAs) and simple general linear models (GLMs), with warming as 0.05, P < 0.05, P = 0.07 & P = 0.06, respectively; Table 2). The negative categorical fixed-effect terms. We then built simple GLMs, with tem­ correlation between the Control vs +6 ◦ C sPLS-DA x-variate and these perature as fixed-effects terms, to identify differences in rarefied mOTU environmental variables imply that Cmic, Nmic, soil DOC and DON abundance and richness due to warming. Changes in the composition of depletion lead to communities with species compositions similar to the hexapod communities for each pair of temperature levels were those found at the +6 ◦ C (Fig. 3). assessed using sparse partial least squares discriminant analysis (sPLS-DA) as implemented in the mixOmics package [38]. sPLS-DA se­ 4. Discussion lects variables and classifies the most discriminative taxa in the com­ munity matrices. The optimal number of components was chosen based Our results indicated that warming causes changes in the 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 Table 1 values based on the area under the receiver operating characteristic Variation of environmental variables across temperatures. curve to complement the sPLS-DA performance results. The main Environmental variable Estimate F P Adjusted R2 sPLS-DA variates (the x-variate, capturing the compositional similarity Microbial C 0.005 ± 0.002 6.26 0.02 * 0.27 of soil hexapod communities) were then subjected to simple GLMs − Microbial N − 0.01 ± 0.00 3.58 0.08 • 0.15 against the set of environmental variables: the amounts of soil DOC and Soil DOC − 0.06 ± 0.02 6.18 0.02 * 0.27 DON, Cmic and Nmic. Soil DON − 0.20 ± 0.07 7.70 0.01 * 0.32 Effect estimates followed by their standard errors are the output of general linear 3. Results models for each environmental variable. The intercept is always the control, and the explanatory variable is always the corresponding environmental variable: 3.1. Environmental variables microbial carbon (C), microbial nitrogen (N), soil dissolved organic carbon (DOC) or soil dissolved organic nitrogen (DON). •, P < 0.1; *, P < 0.05. All The amounts of Cmic, soil DOC and DON decreased with warming, models have 13◦ of freedom. 3 M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 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 exten­ sively distributed across the northernmost part of continental Europe (GBIF, 2021). 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 Fig. 1. Principal component analyses of the distribution of the environmental or abundance [18,39]. Shifts in community composition in our study variables, with 95% confidence ellipses denoting different treatments. 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 con­ trol 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 Hyp­ noidus 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 re­ sources 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 Fig. 2. Violin plots of the numbers of reads (a) and molecular operative taxo­ of Nmic might have followed this loss of C due to the relatively tight C:N nomic units –mOTU- (b) across sites and treatments. Horizontal lines inside stoichiometric balance of microbial metabolism at the site [46–48]. The each plot from the lowest to the highest denote the first, second and depletion of labile DOC and DON in the warmest soils in our study, third quartiles. however, only affected marginally the soil hexapod communities (Table 2, Fig. 3). The reduction of Cmic and DOC should have negatively compositions of the soil hexapod community in a subarctic grassland. In affected important soil hexapod trophic guilds such as microbial feeders this case of study, an increase in soil temperature of 6 ◦ C above ambient and detritivores relying on bacterial and fungal biomass as well as litter had an impact on the community composition, which correlated with input. In contrast, the lack of effect of the moderate +3 ◦ C warming may have been associated with the threshold dynamics of the soil 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 sig­ nificant and marginally significant slope parameters (P < 0.05, and P < 0.1, respectively). See Table 2 for further information on the linear models. 4 M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 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. the impact of warming correlated with lower DOC and DON availability Table 2 as well as with a loss of Cmic and Nmic. The impact of warming was Relationship between the compositional similarities of the soil hexapod com­ therefore quite different from the predicted short-term warming munity with the environmental variables across temperatures. response, which should be driven by the initial enhanced availability of Treatment Environmental Estimate F P Adjusted N. The requirements of hexapods for C and N, and the sources of C and N, variable R2 are of great interest for predicting the response of soil hexapods to Control vs Microbial C 0.000 ± 0.10 0.75 0 warming. In our study, we consider it likely that the loss of Cmic and þ3 ◦ C 0.001 Nmic due to warming depleted the food resources of hexapod bacter­ Microbial N 0.002 ± 0.23 0.64 0 0.004 ivores and fungivores [53], increasing litter concentration and exerting Soil DOC − 0.02 ± 2.42 0.15 0.13 a positive feedback on hexapod detritivores. These differences in C and 0.01 N availability, together with the species-specific temperature optima Soil DON − 0.02 ± 0.24 0.63 0 and sensitivity, may account for the lack of a common compositional 0.04 Control vs Microbial C − 0.002 ± 8.23 0.02 0.44 response across warming levels. The stoichiometric balance of micro­ þ6 ◦ C 0.000 * organisms and hexapods also highlight the importance of the C and N to Microbial N − 0.008 ± 7.38 0.02 0.41 understand the mechanistic response of soil hexapods to warming. The 0.003 * nutrient requirements and elemental stoichiometry of microorganisms Soil DOC − 0.03 ± 4.04 0.07 0.25 and plants have been well studied [54], but our knowledge of hexapod 0.01 • Soil DON − 0.07 ± 4.65 0.06 0.28 stoichiometry is a lot less advanced [55]. For instance, previous work 0.03 • suggests that arthropod stoichiometric response to temperature is not þ3 C vs ◦ Microbial C − 0.001 ± 0.65 0.44 0 necessarily monotonic in shape [56], and that resource quality (e.g. C:N þ6 ◦ C 0.001 and C:P ratios) plays a major role in ground arthropod feeding behavior Microbial N − 0.004 ± 0.70 0.42 0 0.005 [57,58]. Therefore, the complex behavior of arthropods hamper our Soil DOC − 0.02 ± 0.28 0.60 0 ability to assess their nutritional requirements [59,60]. 0.04 Furthermore, it is crucial to consider potential biases in these re­ Soil DON − 0.03 ± 0.12 0.72 0 sponses, as the study system differs from global warming in two signif­ 0.09 icant aspects. Firstly, geothermal warming exhibits a contrasting vertical Results are general linear models, with the first variate (X-variate 1) from an gradient compared to climatic warming, with warmer temperatures sPLS-DA of the hexapod communities as the response variable modeled against increasing from deep soil layers rather than the surface. Previous each environmental variable: microbial carbon (C), microbial nitrogen (N), soil research indicates that such warming influences the composition of soil dissolved organic carbon (DOC) or soil dissolved organic nitrogen (DON). Effect hexapod communities across the soil profile. For instance, in regions estimates are followed by their standard errors. •, P < 0.1; *, P < 0.05. All with a natural gradient of European climatic aridity, specific traits of models have 8◦ of freedom. Note that the higher the value of the sPLS-DA x- soil-dwelling springtail species are favored on the soil surface, whereas variate, the higher the similarity with the communities of the higher tempera­ short-term exposure to drought induces a different response, promoting ture level of the pair (see Fig. 3). traits associated with soil conditions (Ferrín et al., 2021). Secondly, global warming can lead to shifts in vegetation communities, which also microorganisms found in this subarctic grassland, where only increases have implications for hexapod populations. In subarctic regions, one in temperature equal or above 6 ◦ C have been shown to lead to abrupt common response to warming is the expansion of shrubs into grasslands shifts in microbial diversity and composition [49,50]. and tundra, a process known as shrubification. This vegetation change Global warming is generally expected to increase SOM decomposi­ alters nutrient dynamics, microbial activity, and hexapod communities tion leading to higher N availability [51,52]. However, at a decadal scale [61–63]. Therefore, it is essential to interpret these results with caution, 5 M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 considering the unique characteristics of the study system and the po­ [12] C.A. Deutsch, J.J. Tewksbury, R.B. Huey, K.S. Sheldon, C.K. Ghalambor, D.C. Haak, P.R. Martin, Impacts of climate warming on terrestrial ectotherms across latitude, tential influence of these factors. Nonetheless, experimental and mech­ Proc. Natl. Acad. Sci. USA 105 (18) (2008) 6668–6672, https://doi.org/10.1073/ anistic approaches allow the identification of the nutritional demands pnas.0709472105. affecting soil hexapods over a multiplicity of biomes, and by extension [13] F. Sánchez-Bayo, K.A.G. Wyckhuys, Worldwide decline of the entomofauna: a their community responses to warming, clarifying belowground review of its drivers, Biol. Conserv. 232 (2019) 8–27, https://doi.org/10.1016/j. biocon.2019.01.020. ecological dynamics. Future studies on this topic are therefore strongly [14] R.D. Bardgett, D.A. Wardle, Aboveground-Belowground Linkages: Biotic recommended. Interactions, Ecosystem Processes, and Global Change, Oxford University Press, 2010. [15] J. Filser, J.H. Faber, A.V. Tiunov, L. Brussaard, J. Frouz, G. De Deyn, A.V. Uvarov, Declaration of competing interest M.P. Berg, P. Lavelle, M. Loreau, D.H. Wall, P. Querner, H. Eijsackers, J.J. Jiménez, Soil fauna: key to new carbon models, SOIL 2 (4) (2016) 565–582, https://doi.org/ 10.5194/soil-2-565-2016. The authors declare that they have no known competing financial [16] A.E. Douglas, The microbial dimension in insect nutritional ecology, Funct. Ecol. interests or personal relationships that could have appeared to influence 23 (1) (2009) 38–47, https://doi.org/10.1111/j.1365-2435.2008.01442.x. the work reported in this paper. [17] J.F. David, The role of litter-feeding macroarthropods in decomposition processes: a reappraisal of common views, Soil Biol. Biochem. 76 (2014) 109–118, https:// doi.org/10.1016/j.soilbio.2014.05.009. Data availability [18] R.D. Bardgett, W.H. van der Putten, Belowground biodiversity and ecosystem functioning, Nature 515 (7528) (2014) 505–511, https://doi.org/10.1038/ nature13855. Data will be made available on request. [19] J. Filser, The role of Collembola in carbon and nitrogen cycling in soil, ferrinetal_c&nmediatewarmingimpactonsoilhexapods(Originaldata) Pedobiologia 46 (3–4) (2002) 234–245, https://doi.org/10.1078/0031-4056- (Figshare) 00130. [20] D.W. Freckman, Bacterivorous nematodes and organic-matter decomposition, Agric. Ecosyst. Environ. 24 (1–3) (1988) 195–217, https://doi.org/10.1016/0167- Acknowledgment 8809(88)90066-7. [21] P. Taberlet, A. Bonin, L. Zinger, E. Coissac, Environmental DNA: for Biodiversity Research and Monitoring, Oxford University Press, 2018, https://doi.org/10.1093/ This work was supported by the Catalan Government grant SGR oso/9780198767220.001.0001. 2017-10052021-1333, the Spanish Government grant TED2021- [22] L. Zinger, J. Chave, E. Coissac, A. Iribar, E. Louisanna, S. Manzi, V. Schilling, 132627B-I00, funded by Spanish MCIN and by the European Union H. Schimann, G. Sommeria-Klein, P. Taberlet, Extracellular DNA extraction is a fast, cheap and reliable alternative for multi-taxa surveys based on soil DNA, Soil NextGeneration EU/PRTR, and the Fundación Ramón Areces grant Biol. Biochem. 96 (2016) 16–19, https://doi.org/10.1016/j.soilbio.2016.01.008. CIVP20A6621. [23] M.E. Cristescu, From barcoding single individuals to metabarcoding biological communities: towards an integrative approach to the study of global biodiversity, Trends Ecol. Evol. 29 (10) (2014) 566–571, https://doi.org/10.1016/j. Appendix A. Supplementary data tree.2014.08.001. [24] B.D. Sigurdsson, N.I.W. Leblans, S. Dauwe, E. Guðmundsdóttir, P. Gundersen, G. Supplementary data to this article can be found online at https://doi. E. Gunnarsdóttir, M. Holmstrup, K. Ilieva-Makulec, T. Kätterer, B. Marteinsdóttir, M. Maljanen, E.S. Oddsdóttir, I. Ostonen, J. Peñuelas, C. Poeplau, A. Richter, org/10.1016/j.ejsobi.2023.103513. P. Sigurðsson, P. van Bodegom, H. Wallander, I. Janssens, Geothermal ecosystems as natural climate change experiments: the ForHot research site in Iceland as a case References study, Icel. Agric. Sci. 29 (2016) 53–71, https://doi.org/10.16886/IAS.2016.05. [25] O. Arnalds, The Soils Of Iceland, Springer, 2015, p. 180. https://www.springerpr ofessional.de/the-soils-of-iceland/4851244. [1] M. Scholze, W. Knorr, N.W. Arnell, I.C. Prentice, A climate-change risk analysis for [26] Icelandic Meteorological Office, Climate Archives for Weather Stations, 2016. world ecosystems, Proc. Natl. Acad. Sci. USA 103 (35) (2006) 13116–13120, Downloaded on May 2022 from, http://www.vedur.is/vedur/vedurfar/medaltalst https://doi.org/10.1073/pnas.0601816103. oflur/. [2] J.W. Williams, S.T. Jackson, J.E. Kutzbach, Projected distributions of novel and [27] N.I.W. Leblans, B.D. Sigurdsson, P. Roefs, R. Thuys, B. Magnússon, I.A. Janssens, disappearing climates by 2100 AD, Proc. Natl. Acad. Sci. USA 104 (14) (2007) Effects of seabird nitrogen input on biomass and carbon accumulation after 50 5738–5742, https://doi.org/10.1073/pnas.0606292104. years of primary succession on a young volcanic island, Surtsey, Biogeosciences 11 [3] S. Marañón-Jiménez, J.L. Soong, N.I.W. Leblans, B.D. Sigurdsson, J. Peñuelas, (22) (2014) 6237–6250, https://doi.org/10.5194/bg-11-6237-2014. A. Richter, D. Asensio, E. Fransen, I.A. Janssens, Geothermally warmed soils reveal [28] P. Taberlet, S.M. Prud’homme, E. Campione, J. Roy, C. Miquel, W. Shehzad, persistent increases in the respiratory costs of soil microbes contributing to C. Melodelima, Soil sampling and isolation of extracellular DNA from large amount substantial C losses, Biogeochemistry 138 (2018) 245–260, https://doi.org/ of starting material suitable for metabarcoding studies, Mol. Ecol. 21 (8) (2012) 10.1007/s10533-018-0443-0. 1816–1820, https://doi.org/10.1111/j.1365-294X.2011.05317.x. [4] L. Rustad, J. Campbell, G. Marion, R. Norby, M. Mitchell, A. Hartley, [29] A. Kocher, J.-C. Gantier, P. Gaborit, L. Zinger, H. Holota, S. Valiere, I. Dusfour, J. Cornelissen, J. Gurevitch, A meta-analysis of the response of soil respiration, net R. Girod, A.-L. Bañuls, J. Murienne, Vector soup: high-throughput identification of nitrogen mineralization, and aboveground plant growth to experimental ecosystem Neotropical phlebotomine sand flies using metabarcoding, Molecular Ecology warming, Oecologia 126 (4) (2001) 543–562, https://doi.org/10.1007/ Resources 17 (2) (2017) 172–182, https://doi.org/10.1111/1755-0998.12556. s004420000544. [30] S. Talaga, C. Leroy, A. Guidez, I. Dusfour, R. Girod, A. Dejean, J. Murienne, DNA [5] K.K. Treseder, Nitrogen additions and microbial biomass: a meta-analysis of reference libraries of French Guianese mosquitoes for barcoding and ecosystem studies, Ecol. Lett. 11 (10) (2008) 1111–1120, https://doi.org/10.1111/ metabarcoding, PLoS One 12 (6) (2017), e0176993, https://doi.org/10.1371/ j.1461-0248.2008.01230.x. journal.pone.0176993. [6] C. Zhang, Z. Song, D. Zhuang, J. Wang, S. Xie, G. Liu, Urea fertilization decreases [31] E. Coissac, T. Riaz, N. Puillandre, Bioinformatic challenges for DNA metabarcoding soil bacterial diversity, but improves microbial biomass, respiration, and N-cycling of plants and animals, Mol. Ecol. 21 (8) (2012) 1834–1847, https://doi.org/ potential in a semiarid grassland, Biol. Fertil. Soils 55 (3) (2019) 229–242, https:// 10.1111/j.1365-294X.2012.05550.x. doi.org/10.1007/s00374-019-01344-z. [32] F. Boyer, C. Mercier, A. Bonin, Y.L. Bras, P. Taberlet, E. Coissac, obitools: a unix- [7] J.P. Scharlemann, E.V. Tanner, R. Hiederer, V. Kapos, Global soil carbon: inspired software package for DNA metabarcoding, Molecular Ecology Resources understanding and managing the largest terrestrial carbon pool, Carbon Manag. 5 16 (1) (2016) 176–182, https://doi.org/10.1111/1755-0998.12428. (1) (2014) 81–91, https://doi.org/10.4155/cmt.13.77. [33] C. Mercier, F. Boyer, A. Bonin, E. Coissac, SUMATRA and SUMACLUST: fast and [8] N. Eisenhauer, A. Bonn, C.A. Guerra, Recognizing the quiet extinction of exact comparison and clustering of sequences, in: Programs and Abstracts of the invertebrates, Nat. Commun. 10 (1) (2019) 50, https://doi.org/10.1038/s41467- SeqBio 2013 Workshop, Abstract, 2013, pp. 27–29. 018-07916-1. [34] C. Kanz, The EMBL nucleotide sequence database, Nucleic Acids Res. 33 (2004) [9] J.C. Blankinship, P.A. Niklaus, B.A. Hungate, A meta-analysis of responses of soil D29–D33, https://doi.org/10.1093/nar/gki098. Database issue. biota to global change, Oecologia 165 (3) (2011) 553–565, https://doi.org/ [35] L. Zinger, C. Lionnet, A. Benoiston, J. Donald, C. Mercier, F. Boyer, metabaR: an r 10.1007/s00442-011-1909-0. package for the evaluation and improvement of DNA metabarcoding data quality, [10] M. Holmstrup, B.K. Ehlers, S. Slotsbo, K. Ilieva-Makulec, B.D. Sigurdsson, N.I. Methods Ecol. Evol. 12 (4) (2021) 586–592, https://doi.org/10.1111/2041- W. Leblans, J. Ellers, M.P. Berg, Functional diversity of Collembola is reduced in 210X.13552. soils subjected to short-term, but not long-term, geothermal warming, Funct. Ecol. [36] A. Fjellberg, Icelandic Collembola. Revised checklist and general comments, 32 (5) (2018) 1304–1316, https://doi.org/10.1111/1365-2435.13058. Entomol. Scand. Suppl. 64 (2007) 45–60. [11] Y. Peng, J. Peñuelas, L. Vesterdal, K. Yue, G. Peguero, D.A. Fornara, F. Wu, [37] R Core Team, R: A Language and Environment for Statistical Computing, R Responses of soil fauna communities to the individual and combined effects of Foundation for Statistical Computing, 2020, 4.0.0) [Computer software], https:// multiple global change factors, Ecol. Lett. 25 (9) (2022) 1961–1973, https://doi. www.R-project.org/. org/10.1111/ele.14068. 6 M. Ferrín et al. European Journal of Soil Biology 117 (2023) 103513 [38] F. Rohart, B. Gautier, A. Singh, K.-A. Lê Cao, mixOmics: an R package for ‘omics [51] E. Davidson, I. Janssens, Temperature sensitivity of soil carbon decomposition and feature selection and multiple data integration, PLoS Comput. Biol. 13 (11) (2017), feedbacks to climate change, Nature 440 (2006) 165–173, https://doi.org/ e1005752, https://doi.org/10.1371/journal.pcbi.1005752. 10.1038/nature04514. [39] M.J.I. Briones, Soil fauna and soil functions: a jigsaw puzzle, Front. Environ. Sci. 2 [52] D. Jenkinson, D. Adams, A. Wild, Model estimates of CO2 emissions from soil in (2014), https://doi.org/10.3389/fenvs.2014.00007. response to global warming, Nature 351 (1991) 304–306, https://doi.org/ [40] I. Drahun, K.F. Wiebe, P. Gohl, C.W. Koloski, A.J. Koiter, W.G. Van Herk, B. 10.1038/351304a0. J. Cassone, Three years of surveillance associates agro-environmental factors with [53] A. Gargallo-Garriga, J. Sardans, M. Ayala-Roque, B.D. Sigurdsson, N.I.W. Leblans, wireworm infestations in Manitoba, Canada, Pest Manag. Sci. 78 (1) (2022) M. Oravec, K. Klem, I.A. Janssens, O. Urban, J. Peñuelas, Warming affects soil 369–378, https://doi.org/10.1002/ps.6650. metabolome: the case study of Icelandic grasslands, Eur. J. Soil Biol. 105 (2021), [41] W. Li, Y. Yang, W. Xie, Q. Wu, B. Xu, S. Wang, Y. Zhang, Effects of temperature on 103317, https://doi.org/10.1016/j.ejsobi.2021.103317GBIF.org (16 December the age-stage, two-sex life table of Bradysia odoriphaga (Diptera: sciaridae), 2021) GBIF Occurrence Download 10.15468/dl.ukf963. J. Econ. Entomol. 108 (1) (2015) 126–134, https://doi.org/10.1093/jee/tou011. [54] S. Zechmeister-Boltenstern, K.M. Keiblinger, M. Mooshammer, J.P. Uelas, [42] J.M. Alatalo, A.K. Jägerbrand, P. Čuchta, Collembola at three alpine subarctic sites A. Richter, J. Sardans, W. Wanek, The application of ecological stoichiometry to resistant to twenty years of experimental warming, Sci. Rep. 5 (1) (2016), 18161, plant–microbial–soil organic matter transformations, Ecol. Monogr. 85 (2) (2015) https://doi.org/10.1038/srep18161. 23. [43] S.F. Bokhorst, J.W. Bjerke, H. Tømmervik, T.V. Callaghan, G.K. Phoenix, Winter [55] S.J. Moe, R.S. Stelzer, M.R. Forman, W.S. Harpole, T. Daufresne, T. Yoshida, Recent warming events damage sub-Arctic vegetation: consistent evidence from an advances in ecological stoichiometry: insights for population and community experimental manipulation and a natural event, J. Ecol. 97 (6) (2009) 1408–1415, ecology, Oikos 109 (1) (2005) 29–39, https://doi.org/10.1111/j.0030- https://doi.org/10.1111/j.1365-2745.2009.01554.x. 1299.2005.14056.x. [44] N. Verbrigghe, N.I.W. Leblans, B.D. Sigurdsson, S. Vicca, C. Fang, L. Fuchslueger, J. [56] T. Ruiz, A.M. Koussoroplis, M. Danger, J.P. Aguer, N. Morel-Desrosiers, A. Bec, U- L. Soong, J.T. Weedon, C. Poeplau, C. Ariza-Carricondo, M. Bahn, B. Guenet, shaped response Unifies views on temperature dependency of stoichiometric P. Gundersen, T. Kätterer, Z. Liu, M. Maljanen, S. Marañón-Jiménez, K. Meeran, E. requirements, Ecol. Lett. 23 (5) (2020) 860–869, https://doi.org/10.1111/ S. Oddsdóttir, I.A. Janssens, Soil Carbon Loss in Warmed Subarctic Grasslands Is ele.13493. Rapid and Restricted to Topsoil, vol. 25, 2022. [57] M. Jochum, A.D. Barnes, D. Ott, B. Lang, B. Klarner, A. Farajallah, U. Brose, [45] N. Verbrigghe, K. Meeran, M. Bahn, A. Canarini, E. Fransen, L. Fuchslueger, Decreasing stoichiometric resource quality drives compensatory feeding across J. Ingrisch, I.A. Janssens, A. Richter, B.D. Sigurdsson, J.L. Soong, S. Vicca, Long- trophic levels in tropical litter invertebrate communities, Am. Nat. 190 (1) (2017) term warming reduced microbial biomass but increased recent plant-derived C in 131–143, https://doi.org/10.1086/691790. microbes of a subarctic grassland, Soil Biol. Biochem. 167 (2022), 108590, https:// [58] D. Ott, C. Digel, B. Klarner, M. Maraun, M. Pollierer, B.C. Rall, U. Brose, Litter doi.org/10.1016/j.soilbio.2022.108590. elemental stoichiometry and biomass densities of forest soil invertebrates, Oikos [46] S. Marañón-Jiménez, J. Peñuelas, A. Richter, B.D. Sigurdsson, L. Fuchslueger, N.I. 123 (10) (2014) 1212–1223, https://doi.org/10.1111/oik.01670. W. Leblans, I.A. Janssens, Coupled carbon and nitrogen losses in response to seven [59] W.F. Fagan, E. Siemann, C. Mitter, R.F. Denno, A.F. Huberty, H.A. Woods, J. years of chronic warming in subarctic soils, Soil Biol. Biochem. 134 (2019) J. Elser, Nitrogen in insects: implications for trophic complexity and species 152–161, https://doi.org/10.1016/j.soilbio.2019.03.028. diversification, Am. Nat. 160 (6) (2002) 784–802, https://doi.org/10.1086/ [47] J. Séneca, A. Söllinger, C.W. Herbold, P. Pjevac, J. Prommer, E. Verbruggen, B. 343879. D. Sigurdsson, J. Peñuelas, I.A. Janssens, T. Urich, A.T. Tveit, A. Richter, Increased [60] J.A. Hódar, R. Zamora, J. Castro, Host utilisation by moth and larval survival of microbial expression of organic nitrogen cycling genes in long-term warmed pine processionary caterpillar Thaumetopoea pityocampa in relation to food quality grassland soils, ISME Communications 1 (1) (2021) 69, https://doi.org/10.1038/ in three Pinus species: food quality in a plant-insect interaction, Ecol. Entomol. 27 s43705-021-00073-5. (3) (2002) 292–301, https://doi.org/10.1046/j.1365-2311.2002.00415.x. [48] T.W.N. Walker, C. Kaiser, F. Strasser, C.W. Herbold, N.I.W. Leblans, D. Woebken, I. [61] I.C. Barrio, C.G. Bueno, D.S. Hik, Warming the tundra: reciprocal responses of A. Janssens, B.D. Sigurdsson, A. Richter, Microbial temperature sensitivity and invertebrate herbivores and plants, Oikos 125 (1) (2016) 20–28, https://doi.org/ biomass change explain soil carbon loss with warming, Nat. Clim. Change 8 (10) 10.1111/oik.02190. (2018) 885–889, https://doi.org/10.1038/s41558-018-0259-x. [62] B. Jia, L. Jia, X.M. Mou, J. Chen, F.C. Li, Q.J. Ma, X.G. Li, Shrubification decreases [49] D. Radujković, E. Verbruggen, B.D. Sigurdsson, N.I.W. Leblans, I.A. Janssens, soil organic carbon mineralization and its temperature sensitivity in alpine S. Vicca, J.T. Weedon, Prolonged exposure does not increase soil microbial meadow soils, Soil Biol. Biochem. 168 (2022), 108651, https://doi.org/10.1016/j. community compositional response to warming along geothermal gradients, FEMS soilbio.2022.108651. (Fed. Eur. Microbiol. Soc.) Microbiol. Ecol. 94 (2) (2017), https://doi.org/ [63] Z.A. Mekonnen, W.J. Riley, L.T. Berner, N.J. Bouskill, M.S. Torn, G. Iwahana, R. 10.1093/femsec/fix174. F. Grant, Arctic tundra shrubification: a review of mechanisms and impacts on [50] J.T. Weedon, E. Bååth, R. Rijkers, S. Reischke, E. Oddsdottir, J. van Hal, R. Aerts, I. ecosystem carbon balance, Environ. Res. Lett. 16 (5) (2021), 053001, https://doi. A. Janssens, M. van Bodegom, Community Adaptation to Temperature Explains org/10.1088/1748-9326/abf28b. Abrupt Soil Bacterial Community Shift along a Geothermal Gradient on Iceland, vol. 47, 2022, https://doi.org/10.1101/2022.04.04.486966. 7