Effects of 5‐year experimental warming in the Alpine belt on soil Archaea: Multi‐omics approaches and prospects

Abstract We currently lack a predictive understanding of how soil archaeal communities may respond to climate change, particularly in Alpine areas where warming is far exceeding the global average. Here, we characterized the abundance, structure, and function of total (by metagenomics) and active soil archaea (by metatranscriptomics) after 5‐year experimental field warming (+1°C) in Italian Alpine grasslands and snowbeds. Our multi‐omics approach unveiled an increasing abundance of Archaea during warming in snowbeds, which was negatively correlated with the abundance of fungi (by qPCR) and micronutrients (Ca and Mg), but positively correlated with soil water content. In the snowbeds transcripts, warming resulted in the enrichment of abundances of transcription and nucleotide biosynthesis. Our study provides novel insights into possible changes in soil Archaea composition and function in the climate change scenario.


INTRODUCTION
Over the 21st century, temperature changes are affecting ecosystems, communities, and species all over the world, with severe consequences for biodiversity and related ecosystem processes. Particularly Alpine areas, characterized by high environmental harshness, such as low temperature, acidic pH, and limited water content (Donhauser & Frey, 2018), are experiencing warming well above the global average and are expected to be more sensitive to climate change (Diaz et al., 2003;Donhauser & Frey, 2018;Thompson, 2000).
Consequently, rising temperatures may result in changes of soil microbial communities composition and functioning related to biogeochemical processes such as carbon (C) and nutrients cycling Yu et al., 2021;Yuan et al., 2021). These patterns cannot be investigated solely with ecosystems observation studies; indeed, multiple environmental factors may interact with each other, resulting in complex direct and indirect effects on all pools and their fluxes in all ecosystems (Cavaleri et al., 2015). In this light, manipulative studies can distinguish between correlation and causation, which is why field-based climate warming manipulation experiments are frequently used (Cavaleri et al., 2015;Klein et al., 2005;Timonen & Bomberg, 2009). Such field experiments can isolate the effects of individual treatments to test mechanistic hypotheses, revealing important insights about ecosystem responses to extreme or abrupt climatic events (Klein et al., 2005;Timonen & Bomberg, 2009).
Noteworthy, most of the studies on the responses of microbial communities to experimental warming focus on bacteria and fungi, while knowledge on the abundance, diversity, and function of soil archaea, especially in high-altitude environments is still scarce (Hofmann et al., 2016). However, Archaea are omnipresent in terrestrial environments (DeLong & Pace, 2001;Timonen & Bomberg, 2009), with key roles in the global carbon and nitrogen cycles (Leininger et al., 2006). Metabarcoding screening targeting the 16S rRNA gene region showed a relative archaeal abundance of up to 10% of soil prokaryotes (Bates et al., 2011;DeLong & Pace, 2001) peaking in extreme habitats, such as high acidic and high altitude ones (Cavicchioli et al., 2000;Korzhenkov et al., 2019). Recently, it was reported that the relative abundance of archaeal proteins ranged between 0.6% and 6.8% of total soil proteins, with a significant positive correlation with aridity index, which was in turn correlated with a loss of bacterial and fungal diversity (Starke et al., 2021). Methyl-coenzyme M reductase (MCRA), tyrosine-tRNA ligase (SYY), and DNA protection during starvation (DPS) were identified among the most common functions across forest, grassland, and shrubland soils (Starke et al., 2021). MCRA is central to methanogenic pathways (Lyu et al., 2018) and is encoded by strictly anaerobic archaea converting a restricted number of substrates to methane (Hedderich & Whitman, 2006) and thereby potentially controlling methane production across climate (Starke et al., 2021). Not found in the proteins by Starke et al. (2021) but ubiquitous and abundant in the environment is the ammonia monooxygenase (amo) encoded by ammonia-oxidizing archaea (AOA) (Könneke et al., 2005;Treusch et al., 2005;Venter et al., 2004), which play a key role in nitrification (Stahl & de la Torre, 2012), carbon fixation (Hansman et al., 2009), and methane production in oceans (Metcalf et al., 2012). On a higher level of functionality, harsh soils are characterized by a large number of markers involved in the maintenance of cellular and high growth potential (e.g. higher metabolic quotients and genes related to central carbon metabolism, biosynthesis, DNA replication, transcription and translation) in addition to the traits involved in stress tolerance (Lyu et al., 2018).
We used metagenomics and metatranscriptomics to investigate the responses of the soil archaeal community to a 5-year in situ experimental warming (using open top chambers [OTCs]) in Alpine grassland and snowbed (D'Alò et al., 2021). We hypothesized (i) a higher abundance of archaea in snowbeds than in grasslands due to their ability to adapt to the lower temperatures and harsher environmental conditions characterizing the snowbed belt; (ii) an increased abundance of archaea upon warming compared to controls due to their ability to grow in more arid conditions; and (iii) functional differences between vegetation belts and upon warming.
The selected area was located around the Stelvio Pass in the Stelvio National Park, Italian central Alps (46 31 0 N, 10 25 0 E), due to evidence of recent climate change which included vegetational changes since 1953, such as shrub upwards migration and encroachment at the expense of Alpine grasslands and snowbeds (Cannone et al., 2007). The selected Alpine belt (2604-2624 m a.s.l.) comprised two different plant communities: a climax Alpine grassland (Caricetum curvulae) and an Alpine snowbed (Salicetum herbaceae). Both plant communities are highly vulnerable to the impacts of climate change and are at high risk of regression due to the ingression of species from neighbouring communities (Cannone et al., 2007;Cannone & Pignatti, 2014;Malfasi & Cannone, 2021).
To study the effect of short-term warming, a manipulative experiment was conducted in situ (D'Alò et al., 2021Sannino et al., 2022). In 2014, small hexagonal OTCs (2.08 m in diameter) were installed in both Alpine grassland and Alpine snowbed. OTCs passively increased the ambient summer temperature of 1 C (Malfasi & Cannone, 2021) to reach values predicted by the best future warming scenario (RCP2.6) by 2100 (IPCC, 2018). The paired plots outside the OTCs were used as controls, representing current conditions. Five soil samples were collected for each condition on the same day in July 2019 during the peak of the growing season. After removal of plant debris and roots, at least three soil scoops were passed through a 5-mm sterile mesh and mixed to obtain a composite sample, which was collected in Falcon tube. Data for edaphic and environmental variables, bacterial and fungal abundance using qPCR, and raw metagenomic and metatranscriptomic data were obtained from a previously published paper focusing on the abundance, diversity, and function of bacteria and fungi . Briefly, nucleic acid samples for the metagenomics (MG) and metatranscriptomics (MT) were pooled in equimolar volumes and sequenced on an Illumina NovaSeq6000 SP (DS-150) with a 2 Â 150 bp configuration. MG assembly and annotation were performed as previously described (Žifč akov a et al., 2017). Combined assembly was performed using MEGAHIT 1.1.3 (Li et al., 2015). Gene calling was performed using MG-RAST (Meyer et al., 2008). Taxonomic identification was performed both using MG-RAST identification based on Integrated Microbial Genomes (IMG) dataset as a reference and by BLAST against all published fungal genomes available in MycoCosm fungal genomics portal on January 2020 (Grigoriev et al., 2014). Of these two, that taxonomic identification with a higher bitscore was used as the best hit with minimal threshold of bitscore 54. The functions of the predicted genes were annotated with the hmmsearch function in HMMER 3.2.1 (Eddy, 2011) using the FOAM database as a source of HMMs for relevant genes (Prestat et al., 2014). To obtain the number of genes per sample, reads were mapped to the assembled contigs using the Bowtie 2 program (Langmead & Salzberg, 2012). Using the Kyoto Encyclopedia of Genes and Genomes (KEGG) Orthology database (Kanehisa & Goto, 2000), the dominant pathways were annotated at KEGG level 2. All analyses were performed in R (R Core Team, 2020). We extracted the archaeal genes and transcripts from these data as it is commonly done to describe prokaryotic communities (Clark et al., 2021;Kakuk et al., 2021;Souza et al., 2016). More information is available in the study by D' Alò et al. (2022) and the supplementary material S1.

RESULTS AND DISCUSSION
Experimental warming using OTCs successfully increased the ground surface temperatures of 1 C over 5 years ( Figure 1A). Particularly, temperature increased from 12.0 ± 2.1 C in control plots to 12.9 ± 2.0 C under warming in grasslands and from 11.3 ± 2.2 C in control plots to 12.0 ± 2.1 C in snowbeds, both with statistical evidence. The number of genes and transcripts assigned to Archaea in grasslands was not significantly different between control and warmed plots, whereas a weak statistical increase in archaeal abundance was observed in warmed snowbed plots, in both MG and MT ( Figure 1B). Control grassland samples showed a lower abundance of archaea genes (p = 0.2049) and transcripts ( p = 0.0753) than control snowbeds, and warming increased the magnitude of this trend in both genes ( p = 0.0004) and transcripts ( p = 0.0002) in both vegetations (data not shown). Among the physicochemical parameters (Table S1), only water content (WC) showed a slight significant positive Spearman correlation with the abundance of archaeal genes and transcripts, while calcium and magnesium, among the micronutrients, and fungal abundance, showed a significant negative Spearman correlation ( Figure 1C, Figure S1).
The functional profile of soil archaea showed the dominance of metabolic functions, as indicated by the yellow scale in Figure 2, regardless of sequencing method, vegetation belt, or treatment. Lower abundances were found for functions related to cellular processes, environmental information, and genetic information. However, only the genes in grasslands showed an almost identical functional composition in both control and warmed plots. Genes in samples subjected to warming in snowbeds comprised a higher number of functions related to metabolism, particularly of carbohydrates. Transcripts in grasslands showed a F I G U R E 1 Temperature during the growing season in grasslands and snowbeds in control plots and during warming (A), and the relative abundance of archaeal genes and transcripts (B). P values were estimated as paired t-tests. Significant (p < 0.05) spearman correlations of edaphic and environmental parameters as well as the abundance of bacteria and fungi based on qPCR against the relative abundance of archaea in genes (MG) and transcripts (MT) (C).
higher contribution of cellular processes and environmental information at the expense of genetic information and metabolism under warming compared to control plots. Otherwise, transcripts in warmed snowbeds plots showed a higher contribution of functions related to nucleotides and transcription associated to a reduction of cellular processes and environmental information, compared to the control plots. In line with our first hypothesis, we found a higher relative abundance of archaeal genes and transcripts in snowbeds than in grasslands. Indeed, snowbeds are characterized by harsher conditions favouring archaeal presence and activity, due to long-lasting snow cover, longer periods of low soil temperature, and a shorter growing season (D'Alò et al., 2021.
Archaeal sequences abundance further increased with warming only in snowbeds. Aside from the increased soil temperature for both vegetations inside the OTCs, no differences in water content or other physicochemical parameters were observed between warmed and control plots. Therefore, unlike our second hypothesis, the expected increase of archaea in warmer plots does not appear to be a result of increased aridity.
The total archaeal relative abundances in our samples ranged from 0.26% to 0.51% in MG and 0.14% to 0.61% in MT. These values are much lower than those recorded in prokaryotic communities around the world using 16S rDNA gene amplicon sequencing (Bates et al., 2011;DeLong & Pace, 2001). However, in our study, both MG and MT also included sequences from non-prokaryotic organisms and not just bacteria as in mentioned studies, which might explain the much lower relative contribution of archaea.
Archaea resulted to be negatively correlated with fungal biomass and micronutrients, implying that these organisms are more likely to be found in areas with low organic C, nutrient availability, and organic matter turnover (Burke et al., 2012). This appears to be consistent with previous findings in which an increase of fungi was recorded moving from Alpine to sub-Alpine belts, as a result of more favourable environmental conditions .
The positive correlation of the whole pool of archaeal genes and transcripts with water content, on the other hand, seems to contradict the previously observed positive correlation between archaeal proteins and aridity in several ecosystems (Starke et al., 2021). In any case, a direct comparison is questionable because the aridity index is defined as the ratio of mean annual precipitation to mean annual evapotranspiration rather than water content. Water is essential to all life forms, however, linear regression showed only weak coefficients of determination for both genes F I G U R E 2 Relative abundance of KEGG level 2 (L2) functions assigned to archaea in grasslands and snowbeds in control plots and during warming. Cellular processes are shown in red, environmental information in blue, genetic information in green, and metabolism in yellow.
(R 2 = 0.01) and transcripts (R 2 = 0.11), which is why manipulation experiments using different water contents would be necessary to test the validity of this correlation.
A dominance of Halobacteriota, Methanobacteriota, and Thermoproteota ( Figure S2) was recorded. These taxa did not match the previously reported dominance of Nitrososphaerales (Thermoproteota), and Thermoplasmata (Thermoplasmatota) across different terrestrial ecosystems (Bates et al., 2011;Isoda et al., 2017;Neilson et al., 2017;Starke et al., 2021). In this study, not only did the annotation produce fewer archaeal sequences than expected based on a previous research, but the quality of the annotation appeared to be lacking, in particular for non-cultivable taxa. Indeed, the community profile mirrored the profile from published genomes at JGI used for the annotation of the sequences, which is why we decided to remove community composition from this study. This seems to be a common issue of meta-omic studies. Similarly, in a previous metaproteomic analysis of archaea, only the functional part of the metaproteomic data was published (Starke et al., 2021), being the database (https:// www.uniprot.org) to which protein sequences are uploaded biased by the access to culturable strains. For archaeal taxonomy, to the current state, we suggest to use 16S rRNA gene amplicon sequencing instead.
Despite the lack of reliability in the taxonomic results, the general functionality should be equally well described in archaea as the general processes that we described are essential in any living organism. In line with our third hypothesis, we identified functional differences between vegetation belts and upon warming. A higher abundance of transcripts related to transcription and nucleotide biosynthesis during warming in snowbeds were observed, both considered related to strategies for maintaining high cellular growth potential under stress conditions (Diaz et al., 2003;Donhauser & Frey, 2018;Thompson, 2000). This difference was not visible in grassland transcripts, implying a greater impact of warming on archaea in more hostile environments (Diaz et al., 2003;Donhauser & Frey, 2018;Thompson, 2000). Despite the above results, no differences in the genes of both grasslands and snowbeds were observed, highlighting a difference between potential and real activity within the microbiome under different conditions. Apparently, only the active archaeal community changed its metabolic activity upon warming, which was consistent with the differential sensitivity of the soil microbiome in response to drought and forest management . Comparing between studies requires the use of similar functional classifications. In this study, we were only able to use KEGG for functional annotation, while in the work by Starke et al. (2021), the soil archaeal metaproteome was annotated by COG, making their comparison impossible.
Going deeper, investigating methane and nitrogen metabolisms or specific functions such as MCRA was impossible, due to the low number of archaeal sequences obtained. Certainly, these findings revealed the bottlenecks of investigating low abundant organisms in the soil food web, which could improve in the future when more archaeal genomes are sequenced.
To conclude, we investigated the impact of shortterm experimental warming by +1 C on the abundance, structure, and function of archaeal genes and transcripts in grasslands and snowbeds of the Italian Alpine belt. We found an increasing number of archaeal genes and transcripts in snowbeds but not in grasslands; moreover, an increased abundance was recorded upon warming in snowbeds, highlighting that warming is able to impact the soil microbiome, in particular archaeal community, starting from the more vulnerable, harsher, colder altitude site. A functional shift in the active but not the total snowbed community was unveiled during warming in line with the maintenance of high growth potential. However, due to the scarcity of archaea sequences, it was impossible to investigate specific functions, preventing us from going deeper into the ecological consequences of global ecosystem changes.

DATA AVAILABILITY STATEMENT
The datasets MG and MT have been deposited on MG-RAST under accession numbers mgm4901504.3 and mgm4901505.3, respectively.