Local-scale variability in snow chemistry drives distinct microbial communities in Alpine seasonal snowpack

Seasonal snowpack is a dynamic system that accumulates atmospherically transported particles, including dust and biological material such as microorganisms, from both local and long-distance sources. Here, we investigated the relationship between bacterial community composition and chemical properties of seasonal snowpack in mid-altitude valleys of the Swiss Alps. We characterised microbial community structure and diversity in surface snow and in 10–15 cm-deep snow across three adjacent valleys (elevations 1798–2578 m a.s.l.). Analysis of major ions and organic acids was used to identify key environmental factors influencing microbial community composition. While geographical location showed no clear influence on either the chemical composition or the bacterial community composition, bacterial community structure varied significantly with elevation. We identified significant differences between snow depths, with surface snow showing higher diversity than subsurface snow and exhibiting greater cross-site similarity. Surface snow contained twice as many bacterial genera in the core community as underlying snow, with nearly complete overlap. Total inorganic nitrogen and Ca²⁺ were key chemical variables associated with microbial community composition in the snowpack. Using Weighted Gene Co-expression Network Analysis (WGCNA), we identified modules of co-occurring bacterial taxa with distinct responses to these chemical gradients. Genera containing spore-forming taxa (Neobacillus, Niallia, Sporosarcina) were significantly associated with communities with higher nitrogen concentration, while Brevundimonas, Cryobacterium, and Polaromonas were associated with higher-calcium communities. The distinct patterns in chemistry and microbial community structure observed within just 10–15 cm reflect differences in atmospheric sources of deposited precipitation, with additional effects from post-depositional processes. This environmental filtering decreases local diversity while selecting different species based on initial community composition and local conditions. Understanding these bacterial-physicochemical relationships offers insights into how mountain ecosystems adapt to climate-driven changes in snow cover duration and atmospheric conditions.

1 Introduction

Mountain ecosystems comprise various cryospheric elements, including glaciers, permafrost, and snow. Among these, seasonal snow is the most extensive yet ephemeral component, covering the ground for periods ranging from several weeks to months during the cold season (Barry and Gan, 2022; Vaughan et al., 2013). In the European Alps, elevations above 2000 m a.s.l. experience snow cover for about 6.5 months each year, although this duration is decreasing due to climate change (Huss and Hock, 2018; Kosolapova and Altshuler, 2024; Marty et al., 2017). This seasonal snow cover performs multiple functions in alpine ecosystems: it lowers surface albedo (Barry and Gan, 2022; Wendler and Kelley, 1988; Zhang, 2005), acts as a ground insulator (Zhang, 2005), and serves as a water reservoir (Barnett et al., 2005; Fayad et al., 2017).

Furthermore, snow functions as an interface between the ground and atmosphere, collecting dust, microorganisms, and other biological particles on its surface (Xiang et al., 2009). Atmospheric deposition is the main route by which microorganisms colonise winter seasonal snowpack. In the European Alps, this relationship was demonstrated in a study comparing air and snow microbial communities during the accumulation period at Mount Sonnblick (European Alps, 3,106 m a.s.l.) (Els et al., 2019). The origin of atmospheric microorganisms varies considerably with geography (Margesin and Miteva, 2011). In the European Alps, local microbial sources such as vegetation, fauna, and anthropogenic activities contribute more substantially to snow microbiome composition at lower elevations, while high-altitude sites receive greater inputs from long-distance transport (Sanchez-Cid et al., 2022; Segawa et al., 2005; Wunderlin et al., 2016).

Once deposited, microorganisms face the challenging conditions of the snowpack environment (Maccario et al., 2015). Snow is a dynamic oligotrophic medium where atmospheric input constitutes the primary external source of nutrients (Kuhn, 2001). Beyond nutrient limitations, snow microorganisms face low temperatures, limited water availability, and frequent structural changes in the snowpack caused by freeze-thaw cycles. At high altitudes, increased ultraviolet radiation promotes photochemical reactions that produce reactive oxygen species (ROS), which can damage cellular components and impose additional oxidative stress on microbial communities (Ezraty et al., 2017; Grannas et al., 2007). Despite these challenging conditions, some microorganisms remain metabolically active within the snowpack (Amoroso et al., 2010; Carpenter et al., 2000; Holland et al., 2020; Lopatina et al., 2013; Price and Sowers, 2004; Rivkina et al., 2000; Zhu et al., 2020), indicating that these environmental stresses act as selective pressures favouring specific microbial taxa adapted to snow environments (Hell et al., 2013; Keuschnig et al., 2023; Larose et al., 2010; Maccario et al., 2014; Segawa et al., 2005). This post-depositional selection process contributes to the differences in the composition between snowpack microbial communities and their atmospheric sources (Els et al., 2020; Hell et al., 2013), as well as between different layers within the snowpack itself.

The larger scope of research in the field of snow microbiology focuses on polar snowpack (Hell et al., 2013; Keuschnig et al., 2023; Larose et al., 2010; Lopatina et al., 2013; Maccario et al., 2014; Malard et al., 2019). However, mid-latitude alpine ecosystems, including the European Alps, differ significantly from polar regions as air temperatures are higher and frequently fluctuate around the melting point, resulting in rapid changes in snowpack structure, depth, and liquid water availability (van Herwijnen et al., 2024). In recent years, several studies focusing on snow microbial communities in the European Alps have been conducted, primarily at high-altitude sites such as Jungfraujoch (3572 m a.s.l.) and Sonnblick (3106 m a.s.l.) mountains (Els et al., 2020; Fillinger et al., 2021; Wunderlin et al., 2016). These studies were mostly focused on temporal (seasonal and annual) dynamics in community composition and their relation to airborne microbial communities. However, it is unknown how mountain topography impacts snow bacterial communities at the local scale, and how these microbial patterns evolve within freshly deposited snow.

Mountain topography strongly influences the spatial distribution and properties of alpine snowpacks (Helbig et al., 2024; Mott et al., 2018). At large scales, orographic lifting of air masses enhances snowfall with increasing elevation, making altitude a major factor of snow accumulation in mountainous terrain (Helbig et al., 2024; Mott et al., 2018). At smaller ridge and slope scales, interactions between wind flow, terrain, and snow particles generate heterogeneous deposition patterns through processes such as preferential snowfall deposition, which enhances accumulation on leeward slopes and reduces it on windward slopes (Helbig et al., 2024; Lehning et al., 2008, 2011). These topographically driven processes also influence snow chemistry by shaping precipitation pathways, aerosol scavenging, and local dust inputs, leading to spatial variability in ion concentrations and deposition loads across altitudinal gradients (Lafrenière and Sinclair, 2011; Rogora et al., 2006). Such spatial variability in snow chemistry may in turn influence the composition of bacterial communities colonising the snowpack.

In this study, we investigated mid-altitude sites (up to 2500 m a.s.l.) in the Swiss Alps during spring. We sampled across 19 closely located sites (within 25 km) along an altitudinal gradient, spanning three adjacent valleys separated by mountain ridges, but receiving the same snowfall events. This design, combined with comparisons between surface snow (0–2 cm) and the underlying snow sampled at 10–15 cm depth, allows us to explore the roles of geography, elevation, and atmospheric deposition in shaping the chemical composition, bacterial community composition and diversity of the alpine snowpack. We hypothesise that bacterial community composition and chemical signatures would differ between valleys if local geographic factors play a significant role in shaping these variables. However, if atmospheric deposition plays a stronger role, we expect homogeneous bacterial communities and chemical signatures across valleys, particularly in recently deposited surface snow. Finally, by comparing communities separated by mere centimetres within the snowpack, we aim to assess whether post-depositional processes further filter bacterial communities over time, expecting reduced diversity and a shift in composition in older (subsurface) snow.

2 Methods

2.1 Study sites and sampling procedure

Snow samples were collected in late spring (18–19 April 2023) in the Western European Alps (Valais, Switzerland) at elevations between 1798 and 2578 m a.s.l. (Table 1, Fig. 1). At each sampling site, we measured air and snow temperatures, as well as snow depth.

Table 1 List of sampling sites.

Figure 1 Location of the sampling sites. Elevation data obtained from the USGS 3D Elevation Program (3DEP). USGS 3DEP elevation data courtesy of USGS.

A total of 114 samples were collected at 19 sampling sites. Sampling sites were distributed across three Alpine valleys and included several peripheral sites (Table 1), were selected to cover a compact geographic area (within ∼ 25 km) but spanning a substantial elevational gradient. Representative photographs of the sampling sites are provided in the Supplement (Fig. S1). Snow surface conditions were visually assessed at each site. No visible algal blooms or dust deposition were observed, except for site 17, which displayed a visibly darker and more granular snow surface. At each sampling site, snow samples were collected in triplicate from two fixed depth ranges using an ethanol-sterilised scoop: 0–2 cm (referred to as surface) and 10–15 cm (referred to as subsurface). The fixed-depth sampling design was chosen over the stratigraphic layer identification to ensure standardised and rapid sampling of all 19 sites within a 2 d sampling campaign. The depths were selected to capture the differences between snow samples of potentially different ages.

Samples were collected into 450 mL sterile whirl-pack bags, filled to near capacity. One sampling bag was damaged during transportation, resulting in a total of 113 samples. Snow samples remained frozen (with no freeze-thaw) before being transferred to the laboratory and stored at −20 °C until analysed.

2.2 Chemical analysis

Before chemical analysis, ∼ 150 mL of snowmelt water samples (N=113) were gradually thawed at 4 °C and filtered with Sterivex 0.22 µm filter units (EMD Millipore Corporation, USA), alongside a MQ water blank processed under identical conditions as a negative control. Filters were subsequently stored at −20 °C for DNA extraction. The anion (F⁻, Cl⁻, NO${}_{\mathrm{2}}^{-}$, NO${}_{\mathrm{3}}^{-}$, Br⁻, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$) and cation (Na⁺, NH${}_{\mathrm{4}}^{+}$, K⁺, Mg²⁺, Ca²⁺), as well as concentrations of several organic acids (formate, malate, lactate, butyrate, oxalate), were estimated with the use of ion chromatography (Dionex Integrion HPIC System, Thermo Scientific, USA) in the Central Environmental Laboratory (EPFL, Sion). Limits of quantification for each ion are provided in Tables 2 and S1 in the Supplement. Concentrations below the limit of quantification were set to zero for subsequent analyses.

2.3 Amplicon sequencing

Bacterial community composition was assessed with sequencing of the full-length 16S rRNA gene. DNA was extracted from filters with the DNeasy PowerWater Sterivex kit (Qiagen, Germany) according to the manufacturer's protocol. Extraction blanks were included in each batch of DNA extractions (negative controls), amplified and sequenced alongside samples to control for contamination. DNA concentrations were assessed using the Qubit HS DNA kit (Thermo Fisher Scientific, USA). Concentrations were generally low, with many samples below the detection limit (0.1 ng L⁻¹), and median of 0.175 ng L⁻¹ in samples with detectable DNA. For amplification of the full-length 16S rRNA gene, we used primer set 27f (AGAGTTTGATCMTGGCTCAG) and 1427r (CGGTTACCTTGTTACGACTT) (Zorz et al., 2023). The target sequence was amplified with KAPA HiFi HotStart ReadyMix polymerase (Roche, Switzerland). Successful amplification was confirmed by agarose gel electrophoresis, and extraction blanks showed no amplification. PCR products were cleaned up with AMPure beads (Beckman Coulter, USA) according to the manufacturer's protocol. DNA concentrations of PCR products after cleanup were estimated with the Qubit HS DNA kit. The amplicon DNA was barcoded using the Native Barcoding Kit 96 v14 (SQK-NBD114.96, Oxford Nanopore Technologies, UK). Samples were separated into two batches. Long-read sequencing was performed on a MinION Mk1C device (Oxford Nanopore Technologies, UK) using R10.4.1 flow cells (FLO-MIN114).

2.4 Bioinformatic analysis

Basecalling was performed using Dorado v0.7.2 in super accuracy (sup) mode. Adapter sequences were removed using Porechop v.0.2.4. Quality filtering (Q>12) and size selection (1300–1800 bp) were performed with Chopper v.0.9.0. Sequences were then dereplicated and clustered into OTUs using Vsearch v.2.22.1 with 97 % as the threshold for clustering (cluster_size function) (Rognes et al., 2016). Taxonomic classification of OTUs was carried out on the centroid sequences of the OTUs in QIIME 2 v.2024.10.1 (Bolyen et al., 2019) using the Greengenes2 database (McDonald et al., 2024). The phylogenetic tree was built with QIIME 2 using MAFFT alignment. The resulting OTU and taxonomy tables, along with the phylogenetic tree and metadata table, were then analysed with R (v.4.2.2) and were handled as a phyloseq object (package phyloseq v.1.42.0) (McMurdie and Holmes, 2013). Before statistical analysis, we removed all non-bacterial sequences and blank-related contaminants (package decontam v.1.18.0) (Davis et al., 2018). The OTU table was filtered so that each OTU must be present in at least two samples with a total count of at least 10. The rarefaction curve was created with the microeco package (v.1.14.0) (Fig. S2) (Liu et al., 2021).

2.5 Statistical analysis

For statistical analysis, we summed concentrations of nitrogen species (NO${}_{\mathrm{2}}^{-}$, NO${}_{\mathrm{3}}^{-}$, NH${}_{\mathrm{4}}^{+}$) as total inorganic nitrogen (TIN) and combined formate, malate, lactate, butyrate, and oxalate as organic acids. Paired comparisons of chemical composition between snow depths were performed using log-transformed chemical data with the Wilcoxon signed-rank test (package stats v 4.2.2), p-values were adjusted using the Benjamini–Hochberg procedure to control the false discovery rate (FDR). Spearman's correlation (rcorr function, package Hmisc v.5.2-3) and PCA (rda function, vegan package v.2.6-10) were calculated for log-transformed and scaled chemical variables and elevation. Correlation matrices were visualised with the corrplot function (corrplot v.0.95). PCA plot was visualised with ggplot2 package (v.3.5.2.).

Bacterial community composition was analysed with the phyloseq package. Before analysis, the OTU table was Hellinger transformed. Alpha-diversity metrics were calculated using the phyloseq function estimate_richness. Faith's phylogenetic diversity index was calculated with the calculatePD function (package biomeUtils v.0.022). Bray-Curtis dissimilarity matrix was calculated with vegdist function from vegan package (v.2.6-10), and weighted UniFrac distances were calculated with UniFrac function (phyloseq v.1.42.0). Hierarchical clustering was performed with the agnes function (method = “flexible”, par.method = 0.625; package cluster v.2.1.8.2). A core microbiome analysis was performed at the genus level; core genera were defined as those present in ≥ 50 % of samples with a minimum relative abundance of 0.5 %. This prevalence-based threshold was used to identify taxa consistently occurring across samples while accounting for expected spatial variability in snow microbial communities.

Terrain variables were derived from a digital elevation model retrieved using the elevatr package (v0.99.1) in R, projected to UTM zone 32N (EPSG:32632), and extracted at sampling sites using the terra package (v1.8.93). A slope-based resistance surface was constructed using the gdistance package (v1.6.5), in which each pixel was assigned a resistance value equal to its slope in degrees plus 1, and pairwise resistance distances between sites were computed as accumulated least-cost paths across this surface. Haversine geographic distances between sites were calculated using the geosphere package (v1.5.20). Mantel tests and partial Mantel tests were performed using Spearman rank correlations with 9999 permutations (vegan v2.7.2), using Bray-Curtis dissimilarity and weighted UniFrac distances as community dissimilarity measures and Euclidean differences in elevation, slope, Beers aspect index (cos((45° − aspect) $\times \mathit{\pi }/\mathrm{180}\mathit{°}$) +1; ranging from 0 for southwest-facing to 2 for northeast-facing slopes), geographic distance, and resistance distance as predictors; p-values were adjusted for multiple comparisons using the Benjamini–Hochberg false discovery rate (FDR) procedure. Spatial autocorrelation of alpha-diversity metrics was assessed using Moran's I with an inverse geographic distance weight matrix (ape v5.8.1), computed on site-level means. Spearman rank correlations between alpha-diversity metrics and elevation, slope, and aspect were computed separately for each snow depth.

Distance-based redundancy analysis was performed with the capscale function (package vegan v.2.6-10). The best models were selected with an automatic stepwise model selection procedure with the ordistep function (vegan v.2.6-10). Differential abundance analysis at the genus level was performed in relation to TIN and Ca²⁺ concentrations with ANCOM-BC2 (package ANCOMBC v.2.0.3) with p-values adjusted for multiple comparisons using the Holm method. For differential analysis, we used the following parameters: a prevalence threshold of 0.1, a library size cutoff of 1000 reads per sample, and a smoothing parameter of 0.05. Weighted gene co-expression network analysis (WGCNA) was performed with the package WGCNA (v.1.73) (Langfelder and Horvath, 2008). Before analysis, the OTU table was additionally filtered by removing OTUs whose total abundance across all samples was lower than 0.05. The soft thresholding power was set to 9 based on scale-free topology criteria (R²=0.907), and the minimal module size was set to 100 OTUs. To identify genera enriched within each WGCNA module, we performed a hypergeometric test (phyper function, package stats v. 4.2.2). P-values were adjusted for multiple comparisons using the Benjamini-Hochberg method.

3 Results

3.1 Chemical composition of the alpine seasonal snowpack

To assess the spatial and vertical variability in snowpack chemistry, we measured the concentrations of ions, organic acids, and pH in two depths of seasonal snow samples collected from 19 alpine sampling sites in the late spring: the surface snow (top 0–2 cm; n=57) and the subsurface snow (10–15 cm deep; n=56). The last precipitation events in the area occurred a day before (with air masses transported from the east, over Central Europe) and three days before sampling (with the air masses originating from the Atlantic Ocean and transported from the west) (Figs. S3, S4).

Paired comparisons of samples from the same sites showed significant differences in chemical composition between surface and underlying 10–15 cm deep snow (Wilcoxon test, p<0.05; Fig. 2a). Measured ion and organic acid concentrations are summarised in Table 2, with detailed per-sample data provided in Table S1. Total inorganic nitrogen (TIN) was higher in the surface snow, while organic acids, Ca²⁺, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, Mg²⁺, Na⁺ and K⁺ were higher in the subsurface snow. No significant differences were found for pH and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ (Fig. 2a).

Table 2 Chemical composition of surface (0–2 cm) and subsurface (10–15 cm) snow.

Figure 2 Chemical composition comparison between 0–2 and 10–15 cm snow depths. (a) Paired differences in chemical concentrations between 0–2 and 10–15 cm snow depth. Positive values indicate higher concentrations in surface snow (blue dots), and negative values indicate higher concentrations in subsurface snow (orange dots). Statistical significance determined by Wilcoxon signed-rank test with FDR correction. (b) Principal Component Analysis biplot of chemical composition data. Filled circles represent surface (0–2 cm) snow samples, empty circles represent subsurface (10–15 cm) snow samples, with colours indicating the valley of origin. The ${}^{***}$ p<0.001, ${}^{**}p<\mathrm{0.01}$, ${}^{*}p<\mathrm{0.05}$ indicate significance; n.s. indicates not significant. (c, d) Spearman correlations among chemical compounds and elevation for 0–2 cm and 10–15 cm snow depths, respectively. Colour intensity and size indicate correlation strength.

To assess broader patterns in chemical composition, we performed a Principal Component Analysis (PCA) (Fig. 2b). The first two principal components explained 56.78 % of the total variance. Samples were separated by snow depth, with differences primarily driven by total inorganic nitrogen (TIN) content. In contrast, no distinct clustering was observed based on the valley of origin. Notably, six samples from a single site, characterised by high calcium concentrations, formed a separate cluster.

To further explore the relationships among chemical compounds in the seasonal snowpack, we performed Spearman rank correlation analysis separately for each snow depth (Fig. 2c and d). In the subsurface snow, most chemical compounds showed positive correlations with one another. Total inorganic nitrogen (TIN) and organic acids were positively correlated with each other, as well as with SO${}_{\mathrm{4}}^{\mathrm{2}-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, F⁻, while Ca²⁺ exhibited a positive correlation exclusively with Mg²⁺ at both depths. Elevation was positively correlated with TIN, organic acids, and SO${}_{\mathrm{4}}^{\mathrm{2}-}$ at both sampling depths. In contrast, Ca²⁺showed a negative correlation with elevation in subsurfacesnow, indicating compound-specific response to the elevation.

3.2 Bacterial community composition of the seasonal snowpack

To investigate microbial community composition and diversity patterns in seasonal alpine snowpack, we performed full-length 16S rRNA sequencing. After quality filtering (see Methods), which removed low-abundance and potentially erroneous sequences, we detected 16 050 OTUs across all samples.

To characterise the overall microbial community composition, we analysed the relative abundances of taxonomic groups across both snow depths (Fig. 3). Five phyla comprised over 90 % of the total community at both 0–2 and 10–15 cm snow depths, with nearly identical relative abundances: Pseudomonadota (mean ± SE: 45.5 ± 1.8 %, n=57 in surface; 45.3 ± 2.68 %, n=56 in subsurface), Bacillota I (20.4 ± 2.57 % in surface, 24.0 ± 3.3 % in subsurface), Bacteroidota (11.5 ± 1.36 % in surface, 10.4 ± 1.69 % in subsurface), Acidobacteriota (11.3 ± 1.11 % in surface, 10.2 ± 1.34 % in subsurface), Actinomycetota (5.06 ± 0.66 % in surface, 6.43 ± 1.28 % in subsurface).

Figure 3 (a) Microbial community composition across alpine snow samples. Top – Hierarchical cluster analysis using weighted UniFrac distances with flexible linkage method; bottom – Relative abundance of major bacterial phyla in surface and subsurface snow samples. (b, c) NMDS ordination of bacterial communities based on Bray-Curtis dissimilarity for surface (b) and subsurface snow (c), coloured by elevation (m a.s.l.).

To characterise bacterial genera consistently present across sampling sites, we performed core microbiome analysis at the genus level. Core genera were defined as those present in at least 50 % of samples with a minimum relative abundance of 0.5 %. This revealed 17 genera in surface snow and 9 genera in subsurface snow (Fig. S5). The core microbiomes of surface and subsurface snow largely overlapped, with Domibacillus emerging as the only genus exclusive to the subsurface snow core microbiome. These core genera included taxa with known ecological adaptations consistent with survival strategies expected in snow, including several cold-adapted (e.g., Hymenobacter, Sphingomonas) (Tighe et al., 2025), lichen-associated (e.g., Lichenibacterium, Lichenicola) (Touchette et al., 2023), and spore-forming taxa (e.g., Sporosarcina, Domibacillus) (Bauer et al., 2002).

Within-sample diversity analysis revealed systematic differences between snow depths (Fig. 4). Surface snow consistently exhibited higher microbial diversity than subsurface snow across multiple metrics. Shannon diversity index was 13 % higher in the surface snow (mean ± SE: 4.90 ± 0.12) than in the subsurface snow (mean ± SE: 4.34 ± 0.12; Wilcoxon test, p<0.01). Similarly, inverse Simpson diversity showed a significant difference (surface: 47.04 ± 4.53; subsurface: 30.16 ± 3.09; Wilcoxon test, p=0.003), indicating that surface snow not only contains more species but also has more even abundance distributions. At the same time, the differences in Faith's phylogenetic diversity were not statistically significant (p=0.06), suggesting that while surface snow has more species, the evolutionary diversity of communities is similar.

Figure 4 Alpha diversity metrics in surface (0–2 cm) and subsurface (10–15 cm) snow samples. Boxplots comparing Shannon diversity, inverse Simpson diversity, and Faith's phylogenetic diversity between surface and subsurface snow samples. Statistical significance determined by the Wilcoxon signed-rank test (${}^{***}p<\mathrm{0.001}$, ${}^{**}p<\mathrm{0.01}$).

To compare beta-diversity patterns between snow depths, we analysed Bray-Curtis dissimilarity distributions within surface and subsurface snow samples. Microbial community dissimilarity was significantly higher between the deeper samples compared to the surface samples (Wilcoxon test, p<0.001), indicating that while surface samples are similar to each other, subsurface samples are more dissimilar from one another. Analogous comparison of weighted UniFrac distances confirmed this pattern (p<0.001), indicating greater homogeneity in surface snow communities. Despite these structural differences, cluster analysis using both Bray-Curtis and weighted UniFrac distances did not reveal distinct clusters associated with the snow depth or valley of origin (weighted UniFrac, Fig. 3a).

Geographic distance exhibited only a weak correlation with community dissimilarity in surface snow (Bray-Curtis distance: r=0.14, p=0.009), and showed no significant link with subsurface snow phylogenetic composition (weighted UniFrac distance: r=0.05, p=0.26), although a signal persisted for subsurface Bray-Curtis distances (r=0.16, p=0.005) (Table S2 in the Supplement). Resistance distance (distance accounting for terrain barriers) yielded similar results in surface snow (r=0.14, p=0.006) but demonstrated a stronger association in deeper snow (r=0.21, p=0.002), remaining significant after controlling for elevation (partial Mantel test, r=0.15, p=0.002) (Table S2). Slope difference between sites showed an independent association with surface snow community dissimilarity after controlling for both elevation and geographic distance (partial Mantel tests: r=0.08, p=0.046 and r=0.09, p=0.037) but showed no significant association in deeper snow (all p>0.05) (Table S2). Aspect showed no consistent association with community dissimilarity across metrics or sample depths (Tables S2).

Elevation difference between sites was the strongest of investigated geographical predictors of community dissimilarity at both depths, maintaining its signal after controlling for geographic distance (partial Mantel test; surface: r=0.23, p<0.001; subsurface: r=0.20, p=0.002) (Table S2). In contrast, geographic distance signals in surface snow largely disappeared after controlling for elevation (r=0.07, p=0.069) (Table S2). NMDS ordination coloured by elevation further indicated an elevational gradient in community composition, particularly in surface snow (Fig. 3bc).

At the within-sample level, community evenness increased with elevation in both surface (r=0.35, p=0.008) and subsurface snow (r=0.34, p=0.010), whereas Shannon diversity and Faith's phylogenetic diversity showed no significant trends with elevation. Slope and aspect showed no significant correlation with alpha-diversity metrics at either sampling depth (all p-values >0.05). Moran's I revealed no significant spatial autocorrelation in alpha-diversity across sites (all p-values >0.1).

3.3 Environmental drivers of bacterial communities in the seasonal snowpack

To identify environmental predictors of microbial community composition, we performed distance-based redundancy analysis (dbRDA) using both weighted UniFrac and Bray-Curtis dissimilarities. For weighted UniFrac distances, stepwise model selection identified TIN, Ca²⁺, and Mg²⁺ as significant chemical predictors, explaining 17.5 % of phylogenetic variation, with TIN showing the strongest effect (pseudo-F =9.71, p=0.001). Including elevation as a candidate predictor significantly enhanced the model, with only TIN, Ca²⁺, elevation, and organic acids explaining 24.0 % of the variation (Fig. 5a). For Bray-Curtis distances, models were more complex, with the chemical-only model including six predictors (Ca²⁺, SO${}_{\mathrm{4}}^{\mathrm{2}-}$, PO${}_{\mathrm{4}}^{\mathrm{3}-}$, TIN, pH, organic acids) explaining 12.5 % of variation, while adding elevation improved model performance to 14.5 % (Fig. 5b). Marginal effects analysis of the Bray-Curtis model revealed Ca²⁺, elevation, and TIN as the strongest predictors (pseudo-F =2.99, 2.79, and 2.65, respectively, all p=0.001). Together, these results highlight Ca²⁺ and TIN as the chemical variables most consistently explaining variation in bacterial community composition in seasonal snow, showing strong influence across both phylogenetic and taxonomic analyses. Elevation also plays a notable role, particularly in shaping the phylogenetic structure of microbial communities (Fig. 5a).

Figure 5 Environmental drivers of microbial community variation. (a, b) Distance-based redundancy analysis (dbRDA) ordination plots using weighted UniFrac (a) and Bray-Curtis (b) distances. (c, d) Fold-change in genus-level abundance associated with total inorganic nitrogen (TIN) and calcium concentrations (ANCOM-BC2 analysis). (e) Weighted gene co-expression network analysis (WGCNA) heatmap showing correlations between microbial modules (y-axis) and environmental factors (x-axis), values show Pearson correlation coefficients with p-values in parentheses. Colour intensity reflects the direction and strength of the correlation, with positive correlations shown in red and negative correlations shown in blue. Module names are arbitrary colour-based identifiers automatically assigned by the WGCNA package.

To identify specific taxa driving the observed associations with TIN and Ca²⁺ concentrations, we performed differential abundance analysis using ANCOM-BC2. For TIN concentration, we identified 13 significantly differentially abundant genera (Fig. 5c), with 12 genera showing positive associations and only 1 genus showing a negative association. The positively associated taxa included several genera previously identified as core taxa in surface snow (Neobacillus, Sporosarcina, Thermoactinomyces_A) or both snow depths (Niallia_299899 and Pelomonas). At the same time, Telluria_571537 was the only genus negatively associated with nitrogen concentration. Similarly, seven bacterial genera showed increased abundance in calcium-rich snow samples (Fig. 5d).

While differential abundance analysis identified individual taxa responding to variation in snowpack chemistry, we next explored how these taxa are organised into co-occurring groups within the community. To capture these complex environmental interactions, we applied Weighted Gene Co-expression Network Analysis (WGCNA) to identify modules of co-occurring bacterial taxa and their collective responses to chemical variation. The analysis identified nine co-occurrence modules, ranging in size from 105 to 3952 OTUs, with 56 OTUs remaining unassigned to a module (Table 3). Three modules showed strong and significant correlations with chemical compounds ($\left|r\right|>\mathrm{0.5}$, p<0.05; Fig. 5e). Notably, the darkgreen (#5) module (n=3626 OTUs) was positively correlated with Ca²⁺ (r=0.79, p= 3e−25) and was significantly enriched for 19 genera, including three (Brevundimonas, Cryobacterium_381841, and Polaromonas) previously identified in the differential abundance analysis (Table 3). The lightcyan (#1) module (n=964 OTUs) exhibited the strongest association with TIN (r=0.52, p=4e−09) and was enriched in five genera, including three (Neobacillus, Niallia_299899, and Sporosarcina), also detected through the differential abundance analysis (Table 3). Interestingly, the largest module, darkred (#5), showed a significant but weaker negative correlation with elevation ($r=-\mathrm{0.44}$, p=1e−06) and several chemical compounds, including organic acids ($r=-\mathrm{0.26}$, p=0.005) (Fig. 5e). This module was enriched in 16 genera, including seven members of the core microbiome, such as Granulicella_C_415415, Terriglobus_A, as well as lichen-associated genera Lichenibacterium_504423 and Lichenicola (Table 3).

Table 3 Co-occurrence modules identified by WGCNA with associated environmental factors and enriched genera. ns: non-significant.

4 Discussion

In this study, we investigated the relationship between bacterial community composition and chemical content of the seasonal snowpack in mid-altitude valleys in the European Alps. We compared the surface (0–2 cm) to the underlying snow (10–15 cm) across 19 sites spanning an elevation gradient. We demonstrated that both chemical content and bacterial community composition differed between depths, with total inorganic nitrogen (TIN) and calcium concentrations being the variables most strongly associated with both chemical and microbial community differentiation. At the geographical scale, elevation was the dominant predictor of bacterial community dissimilarity. Together, these findings highlight the complex interplay between mountain topography, snow chemistry, and bacterial communities in seasonal alpine snowpack.

Analysis of the major ions and several organic acids concentrations in spring snowpack revealed no clear association between snow chemistry and the valley of origin. This suggests that atmospheric transport of air masses from different sources over medium to long distances plays the predominant role in determining snow chemistry at our study sites. However, one site formed a distinct cluster in the PCA analysis due to elevated calcium concentrations, likely reflecting short-range transport of dust eroded from surrounding exposed slopes (Nickus et al., 1998). In contrast, elevation emerged as a significant predictor of chemical composition, with higher concentrations of organic acids, sulphate, and TIN at higher elevation sites. This pattern may result from the warmer daytime temperatures at lower elevations that may promote melting that redistributes solutes downward through the snowpack, resulting in lower concentrations compared to high-elevation sites. For comparison, a study analysing complete snowpack profiles in the French Alps (1100–3300 m a.s.l.) found the opposite relationship between elevation and concentrations of nitrate and ammonium; however, this trend was not consistent as the Southern Alps did not follow that pattern (Dambrine et al., 2018).

In the broader context, our measurements of the major ion concentrations fell into the expected ranges defined by previous studies on precipitation and snowpack chemical composition in the European Alps, although enrichment was particular for the certain ions. While nitrate concentrations fell into the expected range, ammonium concentrations were elevated compared to the previously reported full-column snowpack measurements across Alps, including multi-year means from the Austrian Alps (Greilinger et al., 2016) and 2005–2006 season snowpack measurements from the Aosta Valley, Italy (Filippa et al., 2010). As our measurements reflect recently deposited surface and subsurface snow rather than whole-snowpack composition, the elevated ammonium concentrations may indicate enhanced ammonium deposition at our sites during spring, likely driven by the fertilizer application and agricultural ammonium emissions from the adjacent Rhône valley (Ammann and Valach, 2024; Filippa et al., 2010; Greilinger et al., 2016; Novak et al., 2025). In contrast, both nitrogen species concentrations far exceeded values observed in remote Arctic snowpacks, highlighting the anthropogenic influence at our sampling sites (Krnavek et al., 2012; Nickus, 2003).

Comparison between surface and underlying 10–15 cm deep snow revealed significant dissimilarities in chemical composition, which likely originated from the differences in chemical composition of various snowfall events characterised by distinct air mass sources. Surface snow, with mean nitrate and ammonium concentrations of 0.58 and 0.24 mg L⁻¹ respectively, is consistent with air masses coming from Central Europe's industrialised and densely populated areas (Fig. S4), explaining the higher concentrations of nitrogen species, which correlated positively with sulphate, commonly associated with anthropogenic pollution (Filippa et al., 2010; Greilinger et al., 2016; Novak et al., 2025). Subsurface snow showed higher levels of chloride and sodium (0.17 and 0.089 mg L⁻¹, respectively), potentially indicative of a preceding precipitation event with stronger Atlantic influence, though no firm attribution to a specific snowfall can be made. These concentrations nonetheless fall within the range previously reported for other Alpine sites, and are substantially lower than those observed at Arctic sites dominated by marine aerosol input. This variability at both sampling depths may originate not only from differences in the initial chemical composition of each snowfall event, but also from processes occurring after deposition. For example, the lower TIN concentrations in the 10–15 cm deep snow may reflect losses of nitrogen species through photolysis and volatilisation caused by prolonged UV exposure (Jacobi and Hilker, 2007; Trachsel et al., 2019), as well as biological uptake of TIN. The strong positive correlation between major ions and organic acids in the deeper snow may reflect post-depositional redistribution processes occurring in the snowpack (Raben and Theakstone, 1994).

Despite these chemical differences, likely due to deposition from alternative air masses, the core bacterial community showed a substantial overlap between snow depths. Nearly all subsurface core genera were also present in surface snow, although surface snow contained twice as many core genera as subsurface snow. This overlap may reflect atmospheric biological filtration, where taxa resilient to freezing temperatures, UV radiation, and desiccation are more likely to survive atmospheric transport (Lappan et al., 2024), resulting in a similar core community across deposition events regardless of atmospheric source (Maccario et al., 2015). Indeed, common bacterial taxa have been identified across geographically distant snow environments, suggesting that snow-adapted ecotypes are globally distributed (Brown and Jumpponen, 2019; Wunderlin et al., 2016), and that different deposition events may draw from a shared pool of airborne bacteria capable of surviving atmospheric transport.

Beyond the shared core, surface and subsurface snow differed in the overall diversity and community structure, with surface snow having higher diversity and greater cross-site similarity compared to the deeper snow. The observed depth-related differences of taxa (particularly low-abundance ones) possibly reflect differences in the initially deposited communities and/or the subsequent post-depositional processes. Post-depositional selection may differentiate these communities over time, reducing local diversity while increasing between-site variability in subsurface snow through selection for taxa adapted to local conditions. Similar selection toward snow-specific communities has been documented in Greenland snow-ice communities and was accompanied by changes in functional profiles (Maccario et al., 2019). Physical redistribution within the snowpack may also have occurred under late-spring Alpine conditions, leading to vertical movement of solutes and possibly microbial cells and environmental DNA (Lazzaro et al., 2015). However, with the available data, it is difficult to disentangle the relative importance of these processes or to determine the extent to which they occurred in the sampled snowpack.

Similar to the chemical composition patterns, the valley of origin showed no clear association with the bacterial community composition. Consistent with the lack of valley effect, neither geographical nor terrain-corrected distance remained significantly associated with community composition after controlling for elevation, highlighting that spatial separation and terrain barriers among sites had little influence on community composition at the local scale (< 25 km). As for the snowpack chemistry, these findings suggest the predominant influence of the large-scale atmospheric transport, with bacterial communities across sites reflecting a shared atmospheric microbial pool rather than local or valley-specific inputs.

Topographic variables such as slope and aspect showed little to no association with the community composition. These results are likely connected, as the ecological relevance of aspect depends strongly on slope angle. On near-flat terrain, north- and south-facing surfaces receive similar solar irradiance regardless of their orientation, and strong microclimatic contrasts between aspects emerge only as slopes steepen (Barry, 1992). With slopes ranging from gentle to moderate (0.8–20°), aspect-driven radiation contrasts across our sites were likely too subtle to drive microbial differentiation. At steeper slopes, where snow redistribution and differential melt between aspects are more pronounced, such contrasts could become sufficient to structure microbial communities (Pomeroy et al., 2004). Aspect-driven effects on snowpack communities may also require sustained differential solar exposure to accumulate over time, through processes such as UV-mediated selection, melt-freeze cycling, and differential nutrient dynamics (Sanchez-Cid et al., 2023). The relatively short snow residence time of studied snow depth may therefore have been insufficient for such signals to emerge.

Among the topographic variables examined, elevation most strongly structured bacterial community composition, with its signal persisting at both sampling depths. Community evenness increased with elevation at both depths while Shannon diversity and Faith's phylogenetic diversity showed no such trend, indicating that elevation reshapes dominance patterns rather than overall diversity. Consistent with this pattern, the largest cluster of co-occurring OTUs was negatively correlated with elevation and enriched in lichen-associated genera, suggesting stronger terrestrial inputs at lower elevation sites closer to vegetated terrain. More broadly, higher evenness at greater elevations may reflect limited local terrestrial inputs and reduced bacterial growth, with communities more closely mirroring the regional atmospheric pool, while at lower elevations lichen-associated deposition and warmer conditions likely favour dominance of fewer taxa. Without direct activity measurements, however, we cannot distinguish whether this pattern results from active ecological processes or differences in depositional sources.

Beyond topographic factors, snowpack chemistry also served as a significant structuring force, with calcium and TIN consistently identified as key chemical predictors of microbial community composition. We identified 964 OTUs associated with elevated TIN levels, with the genera Neobacillus, Niallia, and Sporosarcina enriched in these nitrogen-enriched bacterial communities. All three genera are spore-forming bacteria, an adaptation that enables them to survive in the harsh environmental conditions associated with aerial transport and subsequent deposition in the snowpack (Maccario et al., 2015). Additionally, the genus Sporosarcina includes psychrophilic species that remain active in cold environments (e.g., ice), which may explain their positive correlation with TIN, necessary for metabolic activity under such challenging conditions (Bakermans and Skidmore, 2011; Yadav et al., 2016). In Arctic snowpack, nitrogen was shown to be the limiting factor for bacterial growth and activity, with assimilation serving as the primary N cycling pathway (Larose et al., 2013).

A much larger set of OTUs (3626 OTUs) was associated with high calcium levels, as well as magnesium. These ions are associated with the mineral particles, the concentration of which is positively correlated with the bacterial deposition from the atmosphere to the snowpack (Keuschnig et al., 2023). The high number of OTUs positively associated with calcium may reflect the diversity of bacteria that were deposited together with the dust particles, though the particulate fraction was not assessed and calcium concentration measurements reflect dissolved fraction only. A positive correlation between calcium and bacterial abundance was demonstrated previously in melted snow (Margesin and Miteva, 2011). We identified three genera that were consistently associated with calcium-enriched microbial communities: Brevundimonas, Cryobacterium, and Polaromonas – all previously identified in snow communities (Harding et al., 2011; Hell et al., 2013; Segawa et al., 2005; Yan et al., 2012; Zhang et al., 2010, 2012). Notably, Brevundimonas species have been identified as most abundant in microbial communities from alpine ice caves, a habitat characterised by constant low temperatures and high calcium concentrations, where these bacteria remain metabolically active and participate in the precipitation of calcium carbonate (Lange-Enyedi et al., 2024).

5 Conclusions

Seasonal snowpack serves as a dynamic medium that accumulates atmospheric particles and microorganisms, creating unique microbial habitats in alpine environments. In this study, we examined bacterial community composition and chemical properties in seasonal snowpack across 19 sites in three adjacent Swiss Alpine valleys. Our findings demonstrate that temporal differences in atmospheric deposition and post-depositional processes, more so than geography, drive the observed patterns in chemical composition and community structure within alpine snowpack. Elevation emerged as the main topographical predictor of bacterial community dissimilarity across sites. Surface snow exhibited higher bacterial diversity and greater cross-site similarity compared to subsurface snow at 10–15 cm depth, reflecting distinct atmospheric sources and post-depositional selection processes that operate at centimetre scales. Total inorganic nitrogen and calcium emerged as the chemical variables most strongly associated with community structure, with spore-forming genera (Neobacillus, Niallia, Sporosarcina) associated with nitrogen-enriched surface communities and cold-adapted genera (Brevundimonas, Cryobacterium, Polaromonas) linked to calcium-rich environments.

Future research incorporating functional annotation and assessment of the community's activity may improve our understanding of post-depositional selection and reveal differences and similarities in bacterial activity in alpine snowpack compared to the more extensively studied snowpack of polar regions. Additionally, inclusion of fungi, non-fungal eukaryotes such as protists and algae, and viruses could provide a more comprehensive understanding of microbial diversity and ecological interactions within alpine snowpack ecosystems.