Habitat severity characteristics structure soil communities at regional and local spatial scales along the Antarctica Peninsula

Abstract

Antarctic soils provide an excellent setting to test biogeographical patterns across spatial and environmental scales given their relatively simple communities and the dominance of physical factors that create strong environmental gradients. Additional urgency is given by the fact that their unique terrestrial communities are the subject of conservation efforts in a rapidly changing environment. We investigated relationships of soil community assembly and alpha and beta diversity with climatic and environmental parameters across regional and local scales in Maritime Antarctica. We sampled from a regional gradient of sites that differ in habitat severity, ranging from relatively favourable to harsher physicochemical conditions. At the regional scale, bacterial community characteristics and microarthropod abundance varied along this severity gradient, but most measures of fungal communities did not. Microarthropod and microbial communities differed in which soil and climate parameters were most influential, and the specific parameters that influenced each taxon differed across broad and fine spatial scales. This suggests that conservation efforts will need to focus on a large variety of habitat characteristics to successfully encompass diversity across taxa. Because beta diversity was the result of species turnover, conservation efforts also cannot focus on only the most biodiverse sites to effectively preserve all aspects of biodiversity.

Introduction

Scale is an important consideration in understanding the drivers of terrestrial biogeography. The factors that influence community composition, whether they are biotic interactions with other species or abiotic drivers of the physical environment, may differ in relative importance at local vs regional scales (Usher & Booth 1984, Nemergut et al. 2013, Danis et al. 2020, González-Caro et al. 2021). For example, climate (temperature and precipitation) tends to be an important driver of ecological processes at broad spatial scales, while other biotic and abiotic interactions become more important at fine spatial scales (e.g. Aerts 1997). To preserve biodiversity in a changing world, there is a need to understand the drivers of biological abundance and community composition across scales, including both alpha diversity (diversity within a location) and beta diversity (dissimilarity in species composition across locations).

An unanswered question in macroecology is the extent to which community assembly and diversity reflect either 1) the outcome of environmental filtering, wherein communities comprise a relatively consistent subset of the region's available species pool (nestedness) governed by local conditions (e.g. Lozupone & Knight 2007, Nemergut et al. 2011) or 2) neutral, stochastic drift and dispersal-limited, historically contingent assembly that result in species replacements (turnover) across localities (e.g. Chase 2007, Dickie et al. 2012, González-Caro et al. 2021). Dissimilarity measures such as the Sørenson and Simpson indices can demonstrate whether community dissimilarity increases over spatial distances or environmental gradients and partition that dissimilarity into nestedness and turnover (Baselga 2010). These macroecological concepts are difficult to test at large scales in the field, which limits our understanding of the spatial organization of soil community diversity at different scales, and of whether this organization differs across broad ranges of environmental conditions.

Antarctic soils provide an excellent setting to test biogeographical patterns given the relatively simple communities and dominance of physical factors that create strong environmental gradients (Convey et al. 2014). Low temperatures limit biological processes and water availability, but terrestrial life can be abundant in favourable locations (Wall & Virginia 1999, Convey et al. 2014). In particular, the less extreme climate of Maritime Antarctica, which includes much of the Antarctic Peninsula and the Scotia Arc archipelagos, allows more complex communities to develop compared to the continental regions of Antarctica (Convey & Biersma 2023). Aboveground plant cover is largely limited to mosses (~110 species) and lichens (~250 species; Ochyra et al. 2008, Øvstedal & Smith 2009), with just two indigenous vascular plants (Deschampsia antarctica and Colobanthus quitensis) present only in the Maritime Antarctic region. Megafauna consist, with the exception of one species of scavenging sheathbill, entirely of true marine vertebrates that come to land to breed and moult. Thus, Antarctic terrestrial communities are dominated by small invertebrates and microscopic soil organisms. These communities include a surprising diversity of microbes (bacteria, fungi and algae) and a less diverse community of invertebrates including nematodes, tardigrades, rotifers, springtails, mites and two native species of midge (Adams et al. 2006, Yergeau et al. 2007a, 2012, Cavicchioli 2015, Convey & Biersma 2023). Many are endemic and/or restricted to smaller regions within Antarctica due to dispersal limitation (Lawley et al. 2004, Nielsen et al. 2011, Convey et al. 2018b, 2020). Thus, conservation of Antarctic soil communities inherently conserves a unique and important element of global biodiversity.

For the Antarctic Peninsula, increasing knowledge is becoming available regarding the composition of soil communities. Recent studies in this region have focused on the relationship between soil properties and bacterial communities (e.g. Yergeau et al. 2007a, Ganzert et al. 2011, Gonzalez Garraza et al. 2011, Chong et al. 2012, Dennis et al. 2019), fungal communities (e.g. Lawley et al. 2004, Malosso et al. 2006, Arenz & Blanchette 2011, Dennis et al. 2012, Newsham et al. 2016, Misiak et al. 2021) and the influence of aboveground plants on belowground soil microbial communities (Yergeau et al. 2007a, Teixeira et al. 2010, Delgado-Baquerizo et al. 2018, Zhang et al. 2018, Wentzel et al. 2019, Ball et al. 2022b). Less is known about the relationship between invertebrate communities and soil properties on a broad spatial scale, as most work has been conducted at individual localities or pairings of distant sites. Such research suggests that the abiotic factors influencing these communities include a variety of physical climate factors such as temperature, precipitation and ultraviolet radiation, as well as organic matter content and water availability that are closely related to the aboveground plant community (Bölter 1997, Convey & Wynn-Williams 2002, Convey et al. 2002, Lawley et al. 2004, Yergeau et al. 2007b, Bokhorst et al. 2008, Nielsen et al. 2011, Convey et al. 2014).

The Antarctic Peninsula also provides distinct environmental gradients, enabling investigation of the biotic and abiotic influences on soil communities across scales. In addition to latitudinal and climatic gradients along the > 2000 km length of the Antarctic Peninsula, Ball et al. (2022b) recently described a gradient of soil habitat severity driven not by latitude but by a combination of climate and soil biogeochemical properties including soil nitrogen (N) and other nutrient and metal concentrations, pH, electrical conductivity (EC) and total annual precipitation, as well as measures of mean/maximum air temperature and variability. Habitat severity gradients exist both at the regional scale across the Antarctic Peninsula and at the local scale (within individual sites). Much of the within-site variability corresponds to differences in soil habitats that correlate with plant cover types. Thus, aboveground plant cover also has the potential to interact with habitat severity to structure belowground communities.

Given what is known from other ecosystems, how environmental gradients such as these influence soil communities in Antarctica may depend on the scale in question (e.g. Feeser et al. 2018, Lee et al. 2019). For many ecosystem properties and processes, the environmental factors that dominate at fine spatial scales are not necessarily the same as when considering patterns across broad scales (e.g. climate; Aerts 1997). Yet we lack an understanding of the larger-scale spatial patterning of Antarctic soil communities compared to the role of environmental factors at the local scale (Chong et al. 2015). In addition, studies of Antarctic soil communities typically focus on small numbers of taxonomic groups, mostly at individual localities and rarely at a large spatial extent with any fine-scale resolution. For example, invertebrate studies typically focus on individual taxonomic groups (e.g. mites or nematodes) in relation to particular soil parameters (e.g. moisture or temperature), and fewer studies link soil microbial and invertebrate communities together or address a suite of soil environmental parameters (e.g. Newsham et al. 2004). Working across large spatial scales in Antarctica is logistically challenging, and only a few projects have worked systematically across the Antarctic Peninsula. To date, these have focused on microbial communities (Yergeau et al. 2007a,b, Newsham et al. 2016, 2019, Dennis et al. 2019). Thus, it is not yet clear which environmental properties influence soil biodiversity at broad spatial scales across the Antarctic Peninsula. This strongly limits our ability to predict the consequences of global change for terrestrial communities and identify key areas for conservation.

This study set out to identify patterns in the relationships between soil community assembly and environmental parameters of climate and soil chemistry across scales along the Antarctic Peninsula. Given the rate of environmental change in this region, including increased temperatures (Turner et al. 2009), human presence (Chown et al. 2012) and associated changes in plant and soil communities (Convey & Peck 2019, Siegert et al. 2019), it is important to understand how these parameters structure terrestrial biodiversity to identify indicators of soil biodiversity, pristine reference sites and/or sites that may be particularly sensitive to change. We hypothesized that environmental structuring of the soil community, expressed using a variety of metrics such as community composition, alpha and beta diversity and abundance across multiple different taxa (bacteria, fungi and arthropods), will change over a broad environmental gradient, and that these connections will differ in strength when considered at fine vs broad spatial scales. Specifically, we hypothesize that sample-specific community metrics such as abundance and alpha diversity will decrease with increasing habitat severity at the regional scale, in relation to the parameters that pose biological constraints (temperature and resource limitations). We also expect a similar response for overall community composition, with altered community structure across the regional gradient of habitat severity along the entire Antarctic Peninsula. At the local scale (site level), we predict that habitat severity will similarly drive fine-scale community composition and alpha and beta diversity, but only where there are sufficient within-site environmental gradients. Specifically, our recent work has shown that both sites with very harsh or very favourable conditions lacked sufficient environmental variability to influence communities (Ball et al. 2022b). Furthermore, we hypothesize that, at the regional scale, communities will become increasingly dissimilar (beta diversity) as spatial distances and differences in habitat severity increase, with species turnover (dispersal limitation) explaining beta diversity across latitudes but nestedness (environmental filtering) across the habitat severity gradient. At the local scale, we expect increasing dissimilarity among communities across the habitat severity gradient only at the mid-severity sites, for the same reasons listed above.

Materials and methods

Study sites

Thirteen sites from 10 geographically independent locations covering ~15° of latitude were sampled along the western coast of the Antarctic Peninsula and the Scotia Arc archipelagos (Fig. 1 & Table S1). The transect spans nearly 2300 km from the northernmost point on Signy Island (60°42'S, 45°36'W) to the southernmost point at Sky Blu nunataks in Ellsworth Land (74°51'S, 71°36'W). Sites were selected based upon their location, accessibility and aboveground and/or belowground communities of interest. Specifically, the 13 sampling locations included logistically accessible sites where diverse and abundant plant or soil communities are known (e.g. Mars Oasis and Signy Island), as well as sites with very little knowledge of soil biology (e.g. Elephant and Jenny islands). Sites included two where the presence of non-native species is known (the grass Poa annua at Admiralty Bay and the chironomid midge Eretmoptera murphyi at Signy Island). The 10 latitudinal points include at least one site at every ~1° of latitude between 60°S and 67°S, plus more southerly sites at 71°S and almost 75°S. Further sampling of the southern regions between 67–71°S and 71–74°S was logistically unachievable, but replicate sites were sampled for the more southern locations (three in the vicinity of 67.5°S, two at 71°S and two separate nunataks at 74°S).

Fig. 1. Map showing the study sites along the west coast of the Antarctic Peninsula and in the Scotia Arc archipelagos. The colour of the marker corresponds to where the site falls along the gradient of habitat severity, approximately along the colour spectrum from red (least severe) to violet (most severe). The non-metric multidimensional scaling used to assign habitat severity is in Fig. S1.

Across the South Shetland Islands and South Orkney Islands (all sites north of 63˚S), soils are predominantly weakly developed Typic Haplorthels and Typic Gelorthents, with relatively shallow active layers recorded near the study sites (Beyer et al. 2004, Bockheim et al. 2015). To the south in the Palmer Archipelago, Graham Land and Palmer Land, soils are predominantly Typic Gelorthents followed by Typic and Lithic Humigelepts (Bockheim et al. 2015). The underlying bedrock is a mix of volcanic and granitic types throughout the region (Table S1) underlain by continuous or discontinuous permafrost.

As described by Ball et al. (2022b), these sites represent a gradient of overall habitat severity. Non-metric multidimensional scaling (NMDS) using a subset of non-multicollinear soil physicochemical and climate parameters across the 13 sampling locations summarizes site abiotic differences in one ordination of ‘habitat severity’ (Fig. S1). Of all the soil physicochemical properties measured (including a wide array of elemental composition, mineral nutrient compounds, pH, conductivity and moisture; see the ‘Soil physical and chemical properties’ section below) and the comprehensive collection of site climate data (including temperature and precipitation means and variability, both annually as well as during key seasonal periods; see ‘Climate data’ section below), a subset of 20 parameters were used in the NMDS. Only parameters identified as having variance inflation factors < 5 were included to avoid issues of collinearity that could skew the NMDS. The NMDS therefore depicts a composite ‘map’ of these site differences, demonstrating how sites varied along a gradient of overall habitat severity, ranging from more favourable conditions of warmer temperatures with moderate precipitation and high-nutrient but low-pH soils (upper left quadrant of Fig. S1) to sites that were either wetter or drier, often cooler and low in nutrients (lower right quadrant of Fig. S1). Sites were assigned to a qualitative habitat severity ranking according to the order in which they fall along this trajectory to enable interpretation of the biotic characteristics in the context of their relative severity. Notably, this measure of habitat severity did not change linearly with latitude. While the southernmost site (Sky Blu) was an outlier site on the extreme right of the ordination, the mid-latitude sites at ~67.5°S were the least severe, closely followed by Trinity (63.9°S) and Signy (60.7°S). Habitat severity also varies within sites, though environmental variability differs substantially among the sites. Because climate data describe the entire site (rather than differing for each individual soil sample), the variability within a site is the result of differences in soil physicochemical conditions only. Much of this variation results from the influence of different plant cover types sampled at each site, which are known to correlate with distinct soil characteristics (Yergeau et al. 2007b, Schmitz et al. 2020).

Sampling design

Over the 2014–2015 and 2015–2016 summer seasons, we collected samples at the 13 study sites from beneath different types of vegetation, thus encompassing local-scale variability in habitat severity through the span of microhabitats associated with each plant type. We prioritized locations at each site where plant types could be sampled from within a single area with limited variability in elevation and slope aspect, away from concentrated bird and mammal activity. At each site, we collected samples from beneath five vegetation types (if available): grass (D. antarctica), moss (the most abundant species present at the site, which was predominantly Sanionia uncinata, with Polytrichastrum alpinum or Syntrichia saxicola at some sites), algae (typically Prasiola spp. except at Mars and Ares oases), lichen mix (a mixture of lichens growing on dead moss or grass, given that pure lichen stands were largely found on rock not soil) and bare soils (control, with no visible vegetation cover). All plant cover types were collected from all sites with the following exceptions: Berthelot did not have large enough patches of algae without other plant types, Anchorage and Léonie did not have lichen mix growing on soil, Mars and Ares oases are too far south for grass and Sky Blu contains only bare soils with no macroscopic plant and minimal lichen growth. Because all plant cover types were sampled across most of the latitudinal gradient, regional-scale site differences do not appear to be the result of the presence/absence of a particular plant cover type.

Sites were sampled at one of two intensities (Table S1). At most sites, we collected five replicate samples from each of the focal cover types, referred to as ‘low-intensity sampling’ sites. At five of the study sites with more developed vegetation (Signy, Admiralty, Byers, Biscoe, Anchorage), we collected a greater number of replicate samples to allow more rigorous statistical analyses. Twenty replicate samples were taken from three of the aboveground cover types: grass, moss and bare soil, referred to as ‘high-intensity sampling’. These three cover types represent distinct aboveground influences, including one widespread across the Antarctic Peninsula (moss), the more prolific of the two vascular species (grass) and bare soil without plant cover. At the high-intensity sites, we also conducted low-intensity sampling of five replicate samples for algae and lichen mix.

Before taking each soil sample, we performed a vegetation survey to determine percentage cover of each cover type using a 25 cm × 25 cm quadrat to enumerate the aboveground cover related to the microhabitat differences. Cover type for this purpose included the five focal cover types, as well as rock, detritus and C. quitensis. An ethanol-wiped 10 cm-diameter metal corer was then used to remove a sample of the cover (and associated rhizosphere, if applicable) for later identification of the cover species and extraction for microarthropods. After removal of the cover species, the soil beneath was homogenized to ~10 cm depth, which represents the region of highest biological activity from across the sites that differ in their active layer depth. A subsample was taken for microbial community analysis, and the remaining homogenized soil was then collected for measurement of soil physicochemical properties.

Soil physical and chemical properties

Samples were analysed for a comprehensive array of physicochemical properties including water, nutrients, pH and organic content. The homogenized soil samples were either sieved to 2 mm (for the sites on southern Alexander Island and Sky Blu) or hand-sorted (all other sites) to remove pebbles and rocks. Gravimetric soil water content (SWC) was measured, and the remaining soil was frozen at -20°C for further analysis. For each sample, we followed standard soil methodologies to measure pH and EC on 2:1 and 5:1 di-H₂O:soil dilutions, respectively. Soluble ions were extracted in di-H₂O and analysed for anions (Cl, SO₄, Br, PO₄) on an ion chromatograph (Dionex Corporation model DX-120, Sunnyvale, CA, USA). A subsample of the extract was acidified to 5% HNO₃ and measured for cations (Al, Ba, Ca, Cu, Fe, K, Mg, Mn, Na, P) using inductively coupled plasma optical emission spectroscopy (Spectro Analytical Instruments GmbH & Co. KG, Kleve, Germany). Inorganic N (NO₃ + NO₂-N and NH₄-N) was extracted in 2 M KCl and analysed on a Lachat Autoanalyzer (QC8000; Lachat Instruments, Hach Company, Loveland, CO, USA). Total and organic C and N were analysed on a Perkin-Elmer Elemental Analyzer (PE2400; PerkinElmer, Santa Clara, CA, USA). Some sites had negligible inorganic C, and so only total C is reported. Loss on ignition was measured after 4 h at 550°C in a muffle furnace.

Soil microbial communities

The homogenized soil subsamples were filled with an equal volume of sterilized sucrose lysis buffer then frozen at -80°C for subsequent molecular analyses. DNA was extracted from either 0.7 g of soil for bacterial/archaeal community analysis or 5 g of soil for fungal communities using the cetyltrimethylammonium bromide (CTAB) method (Hall et al. 2008, Mitchell & Takacs-Vesbach 2008). Dual-index paired-end 2 × 301 bp amplicon sequencing on an Illumina (San Diego, CA, USA) MiSeq sequencer was used to characterize the microbial communities. Bacterial (and archaeal, henceforth referred to as ‘bacterial’ for simplicity) community composition was assessed by sequencing the 16S rRNA genes using V6 universal bacterial primers 939F (5′-TTG ACG GGG GCC CGC ACA AG-3′) and 1492R (5′-GTT TAC CTT GTT ACG ACT T-3′). Thermocycling conditions for the 16S amplification consisted of a denaturing step at 94°C for 5 min, 30 cycles of denaturation at 94°C for 30 s, 52°C for 30 s, 72°C for 1.5 min and a final extension at 72°C for 7 min. Fungal community composition was assessed by ITS-2 gene sequencing using the 5.8S-Fun forward primers (5′-GTC TCG TGG GCT CGG AGA TGT GTA TAA GAG ACA GAA CTT TYR RCA AYG GAT CWC T-3′) and ITS4-Fun reverse primers (5′-TCG GCA GCG TCA GAT GTG TAT AAG AGA CAG AGC CTC CGC TTA TTG ATA TGC TTA ART-3′). Thermocycling conditions for fungal amplification consisted of a denaturing step at 96°C for 2 min, 27 cycles of denaturation at 94°C for 30 s, 58°C for 40 s, 72°C for 2 min and a final extension at 72°C for 10 min (Taylor et al. 2016).

Bacterial and fungal datasets were processed separately. Next-generation sequencing of 16S rRNA genes resulted in 7 470 437 reads (average of 15 092 ± 22 331 reads per sample, n = 495 samples). Raw sequence data were demultiplexed and reads were trimmed using Sickle PE v1.33 to a minimum length of 200 bp (Joshi & Fass 2011). Paired-end sequences were merged using PANDAseq v2.11 (Masella et al. 2012). Mean sequence lengths after paired-end merging were 551 (± SD 47) bp. Gene sequences were processed as previously described (Van Horn et al. 2013, Feeser et al. 2018, Ball et al. 2022b) using the Quantitative Insights into Microbial Ecology (QIIME) pipeline (Caporaso et al. 2010b). Bacterial unique operational taxonomic units (OTUs) were clustered by the 97% DNA identity criterion using UCLUST (Edgar 2010) and aligned using the PyNAST aligner (Caporaso et al. 2010a). Taxonomy was assigned to representative OTUs using the Greengenes core set (version 13.8; DeSantis et al. 2006). The dataset was then rarefied to 1000 sequences per sample (random subsampling) to account for uneven sequencing depth and, from that, 119 170 bacterial and archaeal OTUs (97% sequence similarity) were identified from 452 samples. Samples that did not have a minimum of 1000 sequences were excluded from the dataset.

Due to time constraints, only a subset of 76 samples were analysed for fungal genes, including up to five randomly selected samples under three cover types (grass, moss and bare soil) from the five high-intensity sites in addition to Mars Oasis. Sequencing of ITS genes present in these 76 samples resulted in a total of 2 510 422 sequences (average: 33 032; SD: 32 623) and 378 514 OTUs. Sequences were trimmed with Sickle PE and paired-end reads were merged using PANDAseq as with the bacterial dataset. Mean sequence lengths after paired-end merging were 284.7 (± SD 63.1). Fungal sequences were clustered into OTUs using the pick_open_reference_otus.py command within QIIME. This script uses an iterative process that implements closed reference OTU picking against the UNITE reference database (the dynamic release version 8.0, dated 18 November 2018; Nilsson et al. 2018) in addition to de novo OTU picking. Relevant options for this step included --otu_picking_method = uclust anf --min_otu_size = 2 while all default parameters were used, with the exception of pick_otus:enable_rev_strand_match = True. Taxonomic assignments of fungal representative OTUs were made using UCLUST within QIIME (Edgar 2010). After rarefying to 5476 sequences per sample 69 samples remained, comprising 95 587 OTUs. Rarefaction depths were chosen for each dataset to include the largest number of sequences that eliminated no more than 10% of the total samples. Next, the fungal dataset was filtered to exclude all domains except those classified as ‘Fungi’, resulting in a final OTU count of 30 649.

Soil arthropod communities

The cover cores were placed on modified Tullgren funnels for heat extraction of microarthropods. The incandescent light source was gradually increased over the course of 5 days, at which point the vials containing the extracted microarthropods in ethanol were capped. Samples were then identified to the order level and enumerated, and the results were expressed as abundance per unit area based on the diameter of the core. Because at most three taxa of microarthropods are present (Acari, Collembola, Diptera), measures of alpha and beta diversity are not robust and analyses focus on their abundance. Instead, to explore community composition, we calculated the ratio of the two most abundant taxa (springtails and mites). The springtail:mite ratio was calculated by dividing the abundance of springtails by the abundance of mites to demonstrate the degree of springtail dominance in each sample.

Climate data

Climate data were extracted from WorldClim2 (http://worldclim.org/version2; Fick & Hijmans 2017), which provides global air temperature and precipitation data at a high spatial resolution. Data could not be extracted from precise locations for all sites. When necessary, the closest points with similar characteristics were identified and climate data for those points were used instead. This was the case for Signy, Elephant, Admiralty, Biscoe, Léonie, Jenny and Ares. Limited difference is expected for the sites where the neighbouring grid cell was used, given the small climatic differences over such a small distance near sea level. Abbreviations used in the analyses are: mean annual temperature (MAT), mean diurnal range (MDR), temperature seasonality (TSeasonality), maximum temperature of warmest month (TMaxWarm), minimum temperature of coldest month (TMinCold), temperature annual range (TAR), mean temperature of wettest quarter (MTWet), mean temperature of driest quarter (MTDry), mean temperature of warmest quarter (MTWarm), mean temperature of coldest quarter (MTCold), total annual precipitation (TAP), precipitation of wettest month (PWetMo), precipitation of driest month (PDryMo), precipitation seasonality (PSeasonality), precipitation of wettest quarter (PWet), precipitation of driest quarter (PDry), precipitation of warmest quarter (PWarm) and precipitation of coldest quarter (PCold).

Data analyses

Data are publicly available online (Ball et al. 2022a). All statistical analyses were conducted in R (version 3.4.4, The R Foundation, Vienna, Austria). Differences in individual soil chemical and biotic properties were assessed via analysis of variance on each parameter to determine how individual site location influenced soil chemistry and biology. Because habitat severity is defined here as a qualitative ranking according to sites‘ relative position on the ordination and not a numerical factor (which would not be possible from NMDS scores) that can be included as a categorical variable, we focused on analyses of site differences that are discussed in the context of their placement in the habitat severity gradient. Where the effect of site was significant, post hoc Tukey’s tests were conducted to determine which pairs of sites significantly differed from each other.

Microbial community analyses were conducted in the R packages phyloseq (v1.13.0; McMurdie & Holmes 2013) and vegan (v2.5.7; Oksanen et al. 2019). We used the ‘estimate_richness’ function in phyloseq to calculate alpha metrics including Chao1 (an estimate of total richness) and Shannon (effective richness that also factors in the evenness) in both datasets. Patterns in microbial (bacterial and fungal) community composition were analysed via NMDS using Bray-Curtis distance metrics and k = 5, yielding stresses of 0.110 for the bacterial community and of 0.102 for the fungal community. Statistical differences among sites were assessed using the ‘adonis2’ function (McArdle & Anderson 2001). Adonis2 is a permutational multivariate analysis of variance (PERMANOVA; n = 999) that partitions our Bray-Curtis distance matrices among sources of variation (Anderson 2001). When either main effect was significant, a pairwise PERMANOVA was conducted for that main effect using the ‘pairwise.perm.manova’ function of package RVAideMemoire (Herve 2020). NMDS was not, however, an effective method for depicting microarthropod community differences, given that only three taxa were collected (Acari, Collembola and Chironomidae), with not all being present at every site.

Additionally, linear mixed-effects models with site as a random effect were used to test for associations between the soil community (bacterial and fungal alpha diversity and microarthropod abundance) and environmental properties including latitude, climate metrics and soil physical and chemical properties. Models were run using R packages lme4 and lmerTest (Bates et al. 2015, Kuznetsova et al. 2017). An R ² value for the fixed effect of environmental properties was determined using R package MuMIn. R ² values from linear mixed-effects models should be interpreted with caution given the inclusion of the random effect of site that is not incorporated into the traditional R ² calculation of observed vs fitted values of the fixed effect from a traditional linear regression.

Finally, additional beta diversity patterns across the transect were explored by partitioning binary presence/absence OTU dissimilarity matrices into their nestedness and turnover components. Nestedness occurs if the community at a less diverse site is a subset of the diversity at richer sites resulting from loss of taxa from the more diverse community, while spatial turnover occurs by replacement of taxa by other taxa as a result of environmental sorting (Baselga 2010). Following the approach developed by Baselga (2010) and Gutiérrez-Cánovas et al. (2013) to leverage differences inherent in beta diversity indices, we calculated Sørensen‘s index (β_SOR) by measuring total/overall changes in beta diversity, Simpson’s index (β_SIM) by measuring turnover and the difference between Sørensen‘s and Simpson’s indices as the dissimilarity due to nestedness (β_NES = β_SOR - β_SIM). These calculations were performed on the bacterial and fungal communities at the five high-intensity sites using the betapart package (Baselga et al. 2022) at the regional scale. Following Gutiérrez-Cánovas et al. (2013), multiple regression models tested the relationships between these dissimilarity indices and the calculated matrix of Euclidean environmental distances of habitat severity (using the soil physicochemical and climate parameters used in the NMDS in Fig. S1) and of latitude (representing geographical differences). We also repeated these analyses within each of the high-intensity sites (local scale) according to differences in habitat severity. Because the habitat severity index includes climate factors that do not vary within sites, they were excluded from these local-scale analyses. Thus, variations in beta diversity according to habitat severity within sites are driven by only the physicochemical soil properties.

Results

Soil communities across the regional-scale habitat severity gradient

Biotic properties changed across the regional-scale habitat severity gradient (Fig. 2). Bacterial alpha diversity, in particular Chao1 richness, decreased by an order of magnitude with increasing habitat severity. Trinity was the only noticeable outlier in this pattern, with a Chao1 value similar to the more severe Maritime Antarctic sites rather than other similarly less severe sites. Fungal alpha diversity, on the other hand, did not correspond to habitat severity, but was the only biotic property that decreased with increasing latitude (Table I).

Fig. 2. Soil biotic parameters measured at each site across the latitudinal transect. Sites are listed approximately in order of habitat severity and colour-coded according to Fig. 1. Bars represent averages across the samples taken from beneath all plant cover types and error bars represent standard errors. Because properties at the neighbouring sites of Anchorage, Léonie and Jenny islands (A/L/J) were similar, they are presented as a latitudinal average. P-values of the analysis of variance testing for a significant influence of site are presented, and letters represent pairwise comparisons among sites determined by post hoc Tukey's tests.

Table I. R ² values from linear mixed-effects models testing the relationship of biotic metrics (microbial diversity and microarthropod abundance) with latitude and climate at the regional scale. The P-values for each significant relationship are denoted as P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***). See Table S2 for actual P-values. Cell colours denote positive (green) or negative (blue) relationships. Parameters include: mean annual temperature (MAT), mean diurnal range (MDR), temperature seasonality (TSeasonality), maximum temperature of warmest month (TMaxWarm), minimum temperature of coldest month (TMinCold), temperature annual range (TAR), mean temperature of wettest quarter (MTWet), mean temperature of driest quarter (MTDry), mean temperature of warmest quarter (MTWarm), mean temperature of coldest quarter (MTCold), total annual precipitation (TAP), precipitation of wettest month (PWetMo), precipitation of driest month (PDryMo), precipitation seasonality (PSeasonality), precipitation of wettest quarter (PWet), precipitation of driest quarter (PDry), precipitation of warmest quarter (PWarm) and precipitation of coldest quarter (PCold). Metrics with a preceding ⋅ are those included in the non-metric multidimensional scaling of habitat severity.

Bacterial community composition also significantly differed among the sites, though with considerable overlap in communities (Fig. 3a). As with diversity, community differences aligned with habitat severity. The more severe Continental Antarctic sites clustered to the right of the ordination followed by the more severe Maritime Antarctic sites (e.g. Biscoe, Admiralty) on the positive side of the x-axis, while the least severe Maritime Antarctic sites (Jenny, Léonie, Anchorage, Trinity, Signy) largely grouped on the negative side of the x-axis. Samples from only the three high-intensity aboveground types (grass, moss, bare soil) from a subset of the sites were sequenced for fungi, yielding fewer data for analyses of community composition. Fungal communities also significantly differed among sites, and the NMDS shows communities aligned more by latitude, with lower-latitude sites clustering to the left of the x-axis and the higher-latitude sites to the right (Fig. 3b).

Fig. 3. Non-metric multidimensional scaling (NMDS) ordination of a. soil bacterial and b. fungal operational taxonomic units, demonstrating differences in community composition across 13 sites along the latitudinal transect of the Antarctic Peninsula. Each point represents one soil sample, colour-coded according to site. Site names are listed in the legend from least severe (top) to most severe (bottom) along the colour spectrum. The P-values from the operational taxonomic units (‘adonis2’ in R) for site effects are provided. Sites significantly differed in community composition, and the letters next to the legend indicate the results of the post hoc pairwise comparisons, where sites with the same letter have similar communities. Note that fungal communities were only assessed from a subset of the sites.

Microarthropod abundance varied significantly across the transect, representing a ~30× range in abundance (Fig. 2). Notably, midges were only present in samples north of 65°S, though not at every site. Contrasting with the microbial findings, site differences did not correspond to habitat severity, other than very low abundance at the most severe sites. While many of the peaks in microarthropod abundance were at the least severe sites (e.g. Anchorage, Trinity, Signy), abundance was often also high at sites of moderate severity (e.g. Byers, Admiralty, Elephant). Among the three taxa, springtail abundance most closely aligned with habitat severity, and springtails dominated in most samples. Sites ranged by 100× in their degree of springtail dominance (springtail:mite ratio; Fig. 2), although this did not appear to be related to either latitude or habitat severity. Notably, all of these regional-scale patterns across taxa do not correspond to the variation in geological parent material.

Habitat influences on soil communities across scales

The different responses of microbial and arthropod communities to habitat severity reflect the different habitat characteristics to which they correlate at a large spatial scale. We explored these differences by evaluating relationships between biotic parameters and individual physicochemical characteristics. While bacterial alpha diversity was correlated with most of the climate factors identified by VIF for inclusion in the NMDS of habitat severity, fungal alpha diversity and arthropod abundance were not (Table I). At the regional scale, bacterial alpha diversity was promoted by warmer temperatures throughout the year (positively correlated with temperature means, maxima and minima), but it was robust against temperature changes (not correlated with measures of temperature variability) and lack of liquid water (not correlated with most precipitation metrics). Fungal alpha diversity, on the other hand, was sensitive to low precipitation and, overall, was less temperature sensitive, though Chao1 richness was somewhat sensitive to cold temperatures (positively related to MAT and MTCold) and variability. Mite abundance was less robust to temperature variability (therefore being negatively correlated with daily and annual ranges), sensitive to cold temperatures (positively correlated with temperature during cold periods) and associated with consistently wetter conditions (positively correlated with all metrics of precipitation except its seasonality). Springtail and midge abundances were unrelated to climate parameters.

Biotic parameters also differed in their relationships to individual soil chemical parameters, and these relationships differed at regional and local scales (Table II). Across all sites at the regional scale, bacterial alpha diversity was positively correlated with measures of organic matter, moisture and numerous nutrient elements. Fungal alpha diversity was overall unrelated to soil properties, except for the positive relationship between Shannon diversity and two soil micronutrients. Mite and springtail abundances were negatively correlated with pH, which decreases with greater organic matter content in less severe climates, and they were positively correlated with measures of organic matter, macronutrients and several micronutrients. Mite abundance was also strongly correlated with Al content but unaffected by changes in SWC, while springtails were positively correlated with SWC. Given the more limited occurrence of midges, no significant correlations were identified with any of the measured parameters.

Table II. R ² values from linear mixed-effects models testing the relationship between biotic metrics (microbial diversity and arthropod abundance) and soil properties at the regional scale (all sites) with site as a random effect, as well as the linear regression models testing the same relationships at each of the five ‘high-intensity’ sites. Data are only provided for significant relationships, with P-values denoted as P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***). Cell colours denote positive (green) or negative (blue) relationships. Metrics with a preceding ⋅ are those included in the non-metric multidimensional scaling of habitat severity. See Table S3 for full list of R ² values and Table S4 for associated P-values.

EC = electrical conductivity; LOI = loss on ignition; SWC = soil water content.

Within the individual high-intensity sites, however, these relationships with soil habitat parameters differed (Table II). The only persistent pattern at both regional and local scales was the lack of relationships between fungal diversity and soil chemistry, save for a couple of instances with fungal richness (Chao1) at one of the mid-severity sites. Similarities were otherwise very limited, and all biotic parameters were significantly correlated with some factors at individual sites that were not influential at other sites or at the regional scale across sites. At the least severe site (Anchorage), there were no significant correlations between microbial diversity and biotic and soil parameters, while at the more severe sites (Biscoe, Admiralty), there were relatively few significant correlations for springtail abundance.

Beta diversity across scales

Across the high-intensity sites, total bacterial beta diversity (β_SOR) significantly increased with distance along the habitat severity gradient and with increasing latitudinal difference (Fig. 4a). In other words, bacterial communities were increasingly distinct from each other as the physical distance among sites increased and as their difference in habitat severity increased regardless of physical distance. For both latitude and habitat severity, this pattern corresponded to increasing turnover (Simpson dissimilarity index) and decreasing nestedness in the community, suggesting that the change in total beta diversity was the result of species replacements. Although fungal alpha diversity was not related to habitat severity, community dissimilarity (β_SOR) significantly increased with distance along the habitat severity gradient, as well as with latitude. As with the bacterial community, this corresponded to increasing turnover and decreasing nestedness (Fig. 4b). When considered at the local scale, bacterial beta diversity (β_SOR) increased with increasing differences in habitat severity only at the two least severe sites, and fungal beta diversity was not related to habitat severity at any site (Table III). Significant changes in bacterial beta diversity corresponded to increasing turnover, not nestedness.

Fig. 4. Multiple regression models for distance matrices relating a. bacterial and b. fungal beta diversity to dissimilarity in habitat severity and latitude, with all variables being expressed by distance matrices instead of single values of raw data. Total beta diversity (Sørensen dissimilarity) was broken down into its turnover (Simpson dissimilarity index) and nestedness (Sørensen-Simpson dissimilarity) components. Linear relationships are depicted where statistically significant.

Table III. Results of the multiple regression models for distance matrices relating bacterial and fungal beta diversity with differences in latitude (at the regional scale across all sites) and habitat severity (at both the regional and local scales). Total beta diversity (Sørensen dissimilarity) was broken down into its turnover (Simpson dissimilarity index) and nestedness (Sørensen-Simpson dissimilarity) components. Significance levels of the intercept, linear coefficients and R ² values are denoted as: P < 0.05 (*), P < 0.01 (**) and P < 0.001 (***). Linear coefficients can be positive or negative, and, where values are small, they are represented as falling between 0.000 and 0.001 (as < 0.001) or between -0.001 and 0.000 (as > -0.001).

Discussion

At the regional scale, bacterial diversity (alpha and beta) and community composition were strongly influenced by habitat severity, and at the local scale the influence of soil parameters depended on the relative severity of the site. Thus, our hypotheses were supported in the bacterial community. Fungal communities, on the other hand, changed over the regional scale according to latitude, and at the local scale they were not strongly influenced by habitat characteristics. Invertebrate abundance and relative dominance were not directly related to habitat severity, spatial distances or geological differences at the regional or local scale but to various individual physicochemical parameters that differed in importance at the regional vs local scale.

Environmental severity structures belowground assemblages

Our data demonstrate that soil community composition on the Antarctic Peninsula correlates with a complex set of soil and climate parameters (summarized as ‘habitat severity’). While previous work has largely focused on specific taxonomic groups, our data demonstrate that relationships with environmental parameters vary across taxa. A key corollary of this is that biodiversity conservation cannot rely on a single or a small handful of habitat characteristics as a catchall for ‘soil communities’. For example, the characteristics correlated with bacterial diversity (higher temperatures, soil organic matter and micronutrients) were not the same as those for fungal diversity or mite abundance (low temperature variability, higher precipitation and macronutrients). The important factors also differed across sites, and those important at the regional scale were not necessarily important at the local scale, and vice versa.

Only a handful of studies have explored the drivers of soil communities across scales in Antarctica, and those have largely focused on microbial communities. Similar to our findings on the Antarctic Peninsula, bacterial communities in the Victoria Land Dry Valleys varied in the direction and magnitude of their relationships to edaphic characteristics when considered across fine vs broad scales (Van Horn et al. 2013, Feeser et al. 2018). Convey et al. (2014) concluded that water availability is a driving factor of terrestrial diversity at large scale while, at smaller spatial scales, other physical and chemical factors (such as salinity) become important. Chong et al. (2015) suggested that, at local scales, Antarctic communities are structured according to local physicochemical properties, whereas at larger regional scales (> 1000 km), historical influences of sympatric speciation and dispersal limitations explain differences among communities. This is a feature across many other taxonomic groups in Antarctica (Convey et al. 2020, Verleyen et al. 2021). At the broad scale of the entire Antarctic Peninsula, we found that habitat severity structured some aspects of the community but not others. Habitat severity influenced bacterial community composition and alpha diversity and bacterial and fungal beta diversity but not fungal alpha diversity or invertebrate abundance. Notably, available climate data included air temperature rather than soil temperature directly experienced by the soil fauna (see the discussion of the current lack of detailed knowledge of biologically relevant microclimates in Convey et al. 2018a). Had soil temperature data been available, it is possible that greater relationships would have been detected with the other taxa.

In addition, when viewed at a local scale within each site, these same ‘severity’ parameters did not have the same effects. Individual parameters that were important at the regional scale differed across sites at the local scale (Table II). Habitat severity also did not correlate with within-site beta diversity of bacterial communities (except at the two least severe sites) or fungal communities, supporting the prediction that community variation is driven by different factors at local (within-site) vs regional scales (across sites along the entire Antarctic Peninsula). A similar finding of context-dependent relationships between habitat variables and soil biota richness has been identified in other ecosystems, both across sites (Tedersoo et al. 2016, Hendershot et al. 2017, Alzarhani et al. 2019) and across scales (Alzarhani et al. 2019). Perhaps habitat severity was important at the least severe of the high-intensity sites (Anchorage, Signy) because these sites contain a greater span of habitat severity than the more severe sites (Fig. S1), creating enough of a gradient to have a detectable influence. It is also possible that different rates of demographic stochasticity inherent to each location contribute to the different responses to habitat severity within locations. Within the broad-scale response to habitat severity along the entire Antarctic Peninsula, the influence of physicochemical parameters at a local scale may vary due to differing trajectories of local adaptations, genetic drift, historical disturbances, etc. (sensu Alzate et al. 2019). Thus, any given physicochemical parameter may not promote or constrain species similarly across sites.

Our data support the hypothesis that, across the large scale of the entire Antarctic Peninsula (> 1000 km), historical influences, dispersal limitations or other related factors result in turnover of species across increasing spatial separation (Chong et al. 2015). However, we show here that, even at a local scale within a site, differences in microbial communities due to local physicochemical heterogeneity (at mid-severity sites) also result from species turnover rather than chemical and climate properties selecting for subsets of an overall larger population (nestedness) as we had hypothesized. The community differences here are more likely structured by the physicochemical environment as it varies across the local landscape rather than dispersal limitation (given the small distance that dispersal would require), but it does not appear that the environmental factors are selecting for or against members of a larger community (as would be the case with nestedness), instead creating a legacy of the community under each plant type (turnover).

To our knowledge, one other study exists of microbial diversity on the Antarctic Peninsula at a latitudinal scale similar to ours (Dennis et al. 2012, 2019, Newsham et al. 2016, 2018). Newsham et al. (2016) report that fungal Chao1 richness was strongly linked to air temperature. This is supported by our observation of a positive correlation between fungal Chao1 richness and mean temperature and of a negative correlation with temperature seasonality, leading to a significant decline with increasing latitude. However, our data suggest a stronger relationship with precipitation characteristics than temperature, a finding that is also supported at the global scale (Tedersoo et al. 2014). Other Antarctic studies have shown relationships between fungal communities and soil C and nutrient availability that were not strongly detected at either scale in our study (Arenz & Blanchette 2011, Dennis et al. 2012, Lee et al. 2019). Our data also support the observation of Dennis et al. (2019) of a relationship between bacterial alpha diversity and temperature, as well as with other characteristics of habitat severity such as pH (similar to other studies; e.g. Chong et al. 2009a,b).

Our study addressed multiple biological groups. Arthropods differed in their response to habitat severity as compared to microbial communities, while mite and springtail abundances also differed in their response to environmental parameters. In particular, mite abundance was influenced by both climate and soil chemistry parameters, while springtails were correlated with only soil chemistry. At the regional scale, both mites and springtails were positively correlated with similar soil parameters, while at the local scale, the similarities were fewer.

Conserving terrestrial biodiversity on the Antarctic Peninsula

The identification of environmental parameters that could be used as potential indicators of soil biodiversity hotspots would assist in not only biodiversity conservation, but also the identification of pristine reference sites (as discussed in Hughes et al. 2011, Bokhorst et al. 2019) and/or sites that may be particularly sensitive to climate change and/or provide optimal conditions for climate change manipulations. This proves to be a complicated task given the fact that no individual parameter or combination of parameters stood out as key drivers across taxa or locations.

In particular, the influences of habitat severity and latitude on microbial communities were the result of turnover, not nestedness. Beta diversity is therefore probably driven by dispersal limitation across latitudes and the historical influence of plant and soil habitat rather than environmentally deterministic patterns (Baselga 2010). Our data suggest that species exhibit specificity to particular habitats and are adapted to local conditions at any given site (Gutiérrez-Cánovas et al. 2013). It does not appear that stress-tolerant generalists occur everywhere while only less robust taxa are lost as habitat severity increases. Conservation efforts cannot, therefore, focus on just a few high-biodiversity sites. While beta diversity driven by nestedness would permit the prioritization of a small number of the richest sites, the observed dominance of species turnover requires conservation efforts to focus across sites along the entire Antarctic Peninsula, incorporating all aboveground types. Thus, our data strongly support recent calls for an expansion of Antarctica's protected area system to generate greater representativeness (Coetzee et al. 2017, Leihy et al. 2020).

Both microbial and invertebrate communities will be impacted by future environmental change. All taxa were correlated with some aspect of temperature, precipitation and soil chemistry that will probably change under a warming climate, but some taxa will benefit and others will suffer as a result of any given parameter (Convey 2011, Convey & Peck 2019). In polar ecosystems, microarthropod body size and integument characteristics influence sensitivity to environmental change, where small- and soft-bodied organisms tend to be more sensitive to drying and warming associated with environmental change (Convey et al. 2003, Makkonen et al. 2011, Bokhorst et al. 2012, Alatalo et al. 2017). In our results, mites (mostly dominated by small prostigmatid mites) preferred warmer, wetter environments with less variability, suggesting that they will be more sensitive to environmental changes (probably with a positive response). Springtails, on the other hand, are positively correlated with warmer maximum temperatures but are more robust against temperature variability and might suffer under scenarios of warming winters. Springtails are generally viewed as vulnerable to drying that could be associated with warming due to the lack of a waxy layer on their cuticle preventing them from controlling desiccation (Convey et al. 2002, 2003). Species expansions and invasions that alter vegetation will also result in altered soil communities. The ‘greening’ of the Antarctic Peninsula under warmer climates will probably reduce bare soil coverage, with increases in both mosses and flowering plants (Nielsen & Wall 2013, Cannone et al. 2016, 2017, 2022). Such changes may lead to more diverse and abundant soil communities, which may be particularly impactful at mid-severity sites where those plant types host unique microbial communities. The consequences for the biogeochemical functions performed under those communities will need to be explored in order to predict the wider ecosystem impacts of such shifts.