Microorganisms play critical roles in sea ice biogeochemical processes. However, microbes living within sea ice can be challenging to sample for scientific study. Because most techniques for microbial analysis are optimized for liquid samples, sea ice samples are typically melted first, often applying a buffering method to mitigate osmotic lysis. Here, we tested commonly used melting procedures on three different ice horizons of springtime, first year, land-fast Arctic sea ice to investigate potential methodological impacts on resulting measurements of cell abundance, photophysiology, and microbial community structure as determined by 16S and 18S rRNA gene amplicon sequencing. Specifically, we compared two buffering methods using NaCl solutions (“seawater,” melting the ice in an equal volume of 35-ppt solution, and “isohaline,” melting with a small volume of 250-ppt solution calculated to yield meltwater at estimated in situ brine salinity) to direct ice melting (no buffer addition) on both mechanically “shaved” and “non-shaved” samples. Shaving the ice shortened the melting process, with no significant impacts on the resulting measurements. The seawater buffer was best at minimizing cell lysis for this ice type, retaining the highest number of cells and chlorophyll a concentration. Comparative measurements of bacterial (16S) community structure highlighted ecologically relevant subsets of the community that were significantly more abundant in the buffered samples. The results for eukaryotic (18S) community structure were less conclusive. Taken together, our results suggest that an equivalent-volume seawater-salinity buffered melt is best at minimizing cell loss due to osmotic stress for springtime Arctic sea ice, but that either buffer will reduce bias in community composition when compared to direct melting. Overall, these findings indicate potential methodological biases that should be considered before developing a sea ice melting protocol for microbiological studies and afterwards, when interpreting biogeochemical or ecological meaning of the results.
Introduction
Sea ice is a critical component of polar marine ecosystems. Prokaryotic and eukaryotic microorganisms growing in sea ice, particularly at the ice-seawater interface, are significant contributors to polar biogeochemical cycles (Vancoppenolle et al., 2013). Primary production from ice algae is responsible for an estimated 2%–24% of the total primary production in sea ice covered regions, where bottom ice can contain up to 10 times more biomass carbon than the underlying seawater (van Leeuwe et al., 2018). The unique ecology of these microbial communities derives in part from the nature of sea ice. When ice forms from seawater, impurities such as salt and particles are excluded from the growing ice crystals, creating a brine that increases in salinity as temperature decreases and the ice continues to grow. This brine accumulates between growing ice lamellae at the ice-seawater interface, forming a habitat for the sea ice microbial community within the ice matrix (Grossi and Sullivan, 1985; Worster and Wettlaufer, 1997; Junge et al., 2004). As air temperature drops and sea ice grows, the ice develops strong vertical gradients in temperature and brine salinity, with pore spaces in the upper ice becoming colder and more saline (Grossi and Sullivan, 1985; Worster and Wettlaufer, 1997). These conditions impose considerable stress on microbial communities. However, sea ice microbial communities appear to be well adapted to such conditions, with a percentage of the community remaining metabolically active even through the polar night (Junge et al., 2004).
The inherently stratified and heterogeneous nature of the sea ice matrix presents a challenge for sampling internal brine communities. Many methods commonly applied to sea ice microbial communities are designed for liquid samples and require melting the ice before analysis (Miller et al., 2015). Because the bulk salinity of sea ice (melted ice + brine) is lower than the brine salinity of the pore spaces inhabited by the sea ice microbial community, melting can have undesired effects on cell physiology and morphology, including cell lysis due to osmotic stress. Simply melting the ice (“direct” melt) has been shown to result in 13%–97% eukaryotic cell loss, with the highest losses in ciliates and flagellates (Garrison and Buck, 1986; Mikkelsen and Witkowski, 2010), and up to 55% bacterial cell loss in winter ice (Ewert et al., 2013). These losses can alter measurements of dissolved inorganic nutrients and other parameters. Recent work identified a significant difference in dissolved ammonium concentration between buffered and direct melting, which the authors hypothesized to be the result of lysis of ammonium-rich flagellates and ciliates (Roukaerts et al., 2019). An alternative explanation to increased ammonium from direct melting is stimulated bacterial respiration of nitrogen-containing compatible solutes, abundantly released from marine microorganisms during sea ice melt (Firth et al., 2016; Torstensson et al., 2019). The relevant explanation may be influenced by ice-melting time, ice community composition, and sampling location (depth) within the sea ice matrix.
To combat these biases when melting sea ice, several techniques have been implemented by the scientific community to control the thermohaline properties of ice during the melting process and minimize osmotic stress. In addition to direct melt methods, common protocols include melting with “seawater” or “isohaline” buffers. A seawater-buffered melt involves melting an ice section in an equivalent (or larger) volume of either filtered seawater or a NaCl solution prepared to near-oceanic salinity (i.e., 34–35 ppt). An isohaline-buffered melt involves melting an ice section in a specific volume of a more concentrated or hypersaline NaCl solution, so that the resulting meltwater has a salinity similar to the estimated internal brine salinity for that ice horizon (Miller et al., 2015). Internal brine salinities for sea ice are calculated using the measured parameter of ice temperature (Cox and Weeks, 1983).
Mechanical ice crushing, for example, using a hammer on ice core sections in a plastic bag covered by a towel, is sometimes implemented in sea ice melt protocols to reduce melting time (Junge et al., 2004). A shorter melting time is thought to help minimize shifts in community structure during the melting period; rapid melting also appears to reduce impacts on phytoplankton photophysiology (Campbell et al., 2019). We explored an alternative option for reducing the melting time, using a commercial ice shaver, commonly used in the food and beverage industry for producing ice-cold drinks, to shave the sea ice core sections.
The general recommendation is that ice samples should be melted at low temperature using a buffer compatible with the intended biological analysis. Yet, there is little consensus in the literature on which exact melting procedures to employ, making comparing measurements across studies difficult (Miller et al., 2015) and increasing the complexity of designing field campaigns. Recently, Campbell et al. (2019) found that variations in melt protocols can account for over 60% of the variability within Arctic primary production estimates. Additionally, while biases in phytoplankton community structure due to melt procedures have been well documented (Garrison and Buck, 1986; Mikkelsen and Witkowski, 2010), there has been less effort to quantify the effects on bacterial community structure. By measuring both prokaryotic and eukaryotic communities simultaneously, our data provide a broader perspective to fully understand the impacts and inherent biases of melt procedures on analyses of microbial abundance and community structure in sea ice.
We tested the direct, seawater, and isohaline melt procedures on shaved and non-shaved land-fast spring sea ice collected from the Chukchi Sea to quantify their impacts on microbial (eukaryotic and bacterial) abundance, photophysiology, and community structure. Specifically, we sought to identify which melting method was best at minimizing cell loss and whether specific bacterial (16S) or eukaryotic (18S) rRNA gene phylotypes were differentially impacted by the melt procedures. The results of this experiment provide a basis for recommending optimal melt approaches in sea ice research and offer an example of potential biases in such studies.
Methods
Site description and ice core processing
A total of 18 first-year sea ice cores were collected from Chukchi Sea land-fast ice approximately 400 m northwest of the Iḷisaġik College (Naval Arctic Research Laboratory campus) in Utqiaġvik, AK (Figure 1) on April 10, 2019. The cores were collected approximately 15 cm apart with a 9-cm diameter Kovacs ice corer and sectioned in the field using an ethanol-sterilized hand saw. Snow depth was consistent across the sampling site (5 cm) and total ice thickness was measured at 68 cm. Air temperature at the time of sampling was −11°C, and ice temperature was measured immediately after collection of the first core at the top and bottom of the core using a thermistor probe. A linear model fit was then used to estimate temperatures throughout the core (ice-seawater interface: −1.8°C, top: −8.2°C). For all cores, the bottom 50-cm portion of the ice was sectioned into three distinct depth horizons: 0–10 cm (midpoint temperature: −2.3°C), 10–30 cm (midpoint temperature: −3.7°C), and 30–50 cm (midpoint temperature: −5.6°C), where the ice-seawater interface is 0 cm. Sections were sealed in sterile bags (Whirl-Pak™) before transport in a light-shielded cooler to the Barrow Arctic Research Center/Environmental Observatory, with approximately 2.5 hours between collection and the completion of shaving and buffer addition.
Map and images of sampling site. Ice cores were collected from springtime land-fast sea ice near Utqiaġvik, AK, at a coring location (red circle, all panels) located approximately 400 m northwest of Iḷisaġik College, Naval Arctic Research Laboratory campus (blue triangle, panel A). (A) Map created in Ocean Data View version 5.6.2 (Schlitzer, 2022), where dark grey indicates land and shades of brown indicate ocean bathymetry (IBCAO; Jakobsson et al., 2012); inset created using ESRI satellite imagery in QGIS version 3.18.3 (QGIS.org, 2021). (B) Sampling site determined by field notes and images, one example of which is provided, from the date of sampling. Photo credit: JS Bowman.
Map and images of sampling site. Ice cores were collected from springtime land-fast sea ice near Utqiaġvik, AK, at a coring location (red circle, all panels) located approximately 400 m northwest of Iḷisaġik College, Naval Arctic Research Laboratory campus (blue triangle, panel A). (A) Map created in Ocean Data View version 5.6.2 (Schlitzer, 2022), where dark grey indicates land and shades of brown indicate ocean bathymetry (IBCAO; Jakobsson et al., 2012); inset created using ESRI satellite imagery in QGIS version 3.18.3 (QGIS.org, 2021). (B) Sampling site determined by field notes and images, one example of which is provided, from the date of sampling. Photo credit: JS Bowman.
Procedural variations in ice-melting method (direct melting and two buffering methods) were assessed for shaved and non-shaved ice sections across all ice horizons (Figure 2). First, half of the cores (9 cores, 27 horizons) were shaved using an industrial ice shaver (Swan SI-100E) rinsed with Milli-Q water. Then, shaved and non-shaved cores were split evenly among the three melting methods (6 cores, 18 total horizons per method; Figure 2). Samples in the “direct melt” (DM) method received no addition of buffer for an approximate final melt volume of 700 mL in the 0–10 cm ice horizon and 1300 mL in the 10–30 cm and 30–50 cm horizons. In the “isohaline melt” (IM) method, a volume of autoclaved 250 ppt reagent grade NaCl (Fisher Scientific™) solution (stored at 5°C) was added to each sample to reach a final melt salinity matching the estimated brine salinity for each ice horizon: 43 ppt, 67 ppt, and 95 ppt, respectively. Brine salinities were estimated as a function of ice temperature. Volumes for brine addition were calculated using the equations outlined in Cox and Weeks (1983) as a function of ice temperature and density (calculated from bulk ice salinity and temperature; Table S1). Bulk ice salinity was measured using a refractometer on a directly melted subsample of the midpoint of each ice section prior to further treatment. Specifically, a 5 cm, 20 cm, and 40 cm subsample were collected into individual 50 mL falcon tubes from a representative ice core at the start of coring and the small volumes allowed for rapid melt at room temperature. Samples in the “seawater melt” (SM) method were weighed on a balance and an equivalent mass of autoclaved 35-ppt NaCl solution (stored at 5°C) was added, reaching a volumetric dilution factor of 2.06. Approximate final melt volumes were 1300 mL for 0–10 cm ice sections and 2700 mL for 10–30 cm and 30–50 cm sections. Following processing at room temperature, laboratory lighting was turned off and samples were left to melt in their original sample bags for approximately 5 hours before bags were moved into dark refrigerators (1–5°C) to finish melting overnight. Samples were observed over the course of the next morning (April 11, 2019) and subsampled for the following measurements as quickly as possible following the complete melt of each section.
Schematic of the experimental set-up. The bottom 50 cm of 18 ice cores were sectioned into 3 horizons (0–10 cm, 10–30 cm, and 30–50 cm, labeled as 5, 20 and 40, respectively), which were distributed evenly among treatments. Half of the ice sections (27 of 54, indicated by the dashed box) were shaved using an industrial ice shaver. Eighteen sections (9 shaved and 9 non-shaved) were used in each of the 3 melting methods. Direct melt (DM) sections were melted with no osmotic buffer added; isohaline melt (IM) sections were melted in an individualized volume of 250 ppt NaCl solution to reach a final salinity matching estimated internal brine salinities; and seawater melt (SM) sections were melted in an equivalent mass of 35 ppt NaCl solution. All cores were processed further in the same manner.
Schematic of the experimental set-up. The bottom 50 cm of 18 ice cores were sectioned into 3 horizons (0–10 cm, 10–30 cm, and 30–50 cm, labeled as 5, 20 and 40, respectively), which were distributed evenly among treatments. Half of the ice sections (27 of 54, indicated by the dashed box) were shaved using an industrial ice shaver. Eighteen sections (9 shaved and 9 non-shaved) were used in each of the 3 melting methods. Direct melt (DM) sections were melted with no osmotic buffer added; isohaline melt (IM) sections were melted in an individualized volume of 250 ppt NaCl solution to reach a final salinity matching estimated internal brine salinities; and seawater melt (SM) sections were melted in an equivalent mass of 35 ppt NaCl solution. All cores were processed further in the same manner.
Cell abundance and photophysiology
Samples for cell abundance (flow cytometry) and photophysiology (active fluorescence) measurements were collected at room temperature in ambient laboratory lighting. After gentle mixing, an aliquot of 4 mL of meltwater was transferred to a polycarbonate cuvette and run immediately on an AquaFlash Handheld Active Fluorometer (Turner Designs™) to measure chlorophyll a concentration (µg L−1) and effective quantum yield according to the following equation, where Fm represents the fluorescence maximum and F0 represents the fluorescence minimum:
We aliquoted 2 mL of meltwater for measurements of cell abundance using a flow cytometer. These samples were stored in plastic cryovials, fixed to a final concentration of 0.25% glutaraldehyde, and shipped frozen on liquid nitrogen to Scripps Institution of Oceanography (SIO) at the University of California San Diego. At SIO, they were stored at −20°C until processed (April 23, 2019) on a Cyflow Space flow cytometer (Sysmex America Inc.™). Total cell abundance was measured by staining 1 mL of sample with the nucleic acid stain SYBR-green (Molecular Probes™) which binds to all DNA. Additionally, autofluorescent cell abundance was measured by running the remaining 1 mL of unstained sample to count only the photosynthetic, chlorophyll-containing microbial community. Autofluorescent cell abundance was measured in addition to total cell abundance to assess whether the photosynthetic population was any more or less sensitive to osmotic stress and buffer type than the rest of the microbial community. Forward scatter, side scatter, and fluorescence channels FL1 (ex/em: 488/536 nm), FL2 (ex/em: 488/590 nm), FL3 (ex/em: 488/675 nm), FL4 (ex/em: 488/748 nm), FL5 (ex/em: 405/455 nm), and FL6 (ex/em: 638/748 nm) were used to characterize and count all (DNA-containing) and autofluorescent (chlorophyll-containing) cells. A standard volume of 1:10 diluted 123-eCount beads (Fisher Scientific™) were used to calculate absolute cell counts. Both absolute cell counts (autofluorescent and total cells) were then normalized to their respective melt-method dilution factors: 1 for the direct-melt method, 2.06 for the seawater-melt method, and between 1.2 and 1.6 for the isohaline-melt method (Table S1).
DNA collection, extraction, and sequencing
Immediately following subsampling for cell abundance and photophysiology, a volume of 500–2000 mL (section- and biomass-dependent) of meltwater was filtered through a sterile 0.2 µm Supor membrane disc filter (Pall Corporation), using a peristaltic pump, for subsequent 16S/18S rRNA gene amplicon sequencing. Tubing and filter housing were flushed with Milli-Q between each melting method/core section combination. Filters were stored at −20°C until transport to SIO on liquid nitrogen. At SIO, filters were stored at −80°C until extraction on June 25, 2019. Extractions were performed on a KingFisherTM Flex Purification system using a MagMax Microbiome Ultra Nucleic Acid Extraction kit (ThermoFisher Scientific). Amplicon library preparation and sequencing were conducted at the Argonne National Laboratory using a 2 × 150 bp library architecture on the Illumina MiSeq platform. Universal primers 515F and 806R (16S, V4; Walters et al., 2016) and 1380F and 1510R (18S, V9; Amaral-Zettler et al., 2009) were used for rRNA gene amplification. Illumina reads were denoised and merged using the R package dada2 (Callahan et al., 2016). The resulting amplicon sequence variants (ASVs) were used as input for the paprica v0.7.0 pipeline, a phylogenetic placement approach for determining community structure and metabolic potential (Bowman and Ducklow, 2015; Erazo et al., 2021; https://github.com/bowmanjeffs/paprica).
For Bacteria and Archaea, paprica places sample reads onto a reference tree created from the concatenated 16S and 23S rRNA genes from all completed genomes in the Genbank RefSeq database (Haft et al., 2018). Each ASV is placed to either an internal or terminal branch on the tree (Nawrocki and Eddy, 2013; Barbera et al., 2019; Czech et al., 2020). Taxonomically, ASVs assigned to terminal completed genomes are described by their complete taxonomic name, while ASVs assigned to internal estimated genomes are described by the consensus taxonomy for all downstream nodes. “Placement proportions” represent the percentage of times a unique ASV was assigned to a specific terminal branch. ASVs with lower placement proportions were placed to multiple nodes, indicating lower confidence in the final classification. Regardless, the taxonomic associations derived from terminal placements should not be taken literally; our use here indicates only that the association is the closest relative among genomes in RefSeq. The point of placement in the reference tree is then used to estimate genome size, 16S rRNA gene copy number, and GC content (Bowman and Ducklow, 2015; Karp et al., 2021). The 16S rRNA gene copy number was used to normalize the abundance data and combined with total cell counts from flow cytometry data to calculate normalized read counts for each ASV. For Eukarya, paprica places ASVs onto a reference tree created from the unique 18S rRNA genes present in the PR2 4.13.0 database (Guillou et al., 2013). The 18S rRNA gene ASVs are represented by relative abundance only and were not normalized to 18S rRNA gene copy. ASVs with an abundance of 1, samples with bad library builds (<5000 total ASV counts), and suspected mitochondria/chloroplasts were filtered out prior to all further downstream analyses.
Statistical analyses
All statistical analyses were carried out using the R programming language (version 4.1.2) in R Studio (R Core Team, 2021). To assess differences between melting methods or ice shaving on cell abundance (log10 transformed), chlorophyll a content (log10 transformed), effective quantum yield, and microbial diversity (inverse Simpson Index), we first ran a two-way crossed ANOVA with an error factor accounting for ice horizon using all samples collected (vegan package; Oksanen et al., 2015). When statistically significant differences were encountered for either of these procedural variations (with no significant interaction between them), or to further investigate fine scale variability within ice horizons (0–10 cm, 10–30 cm, and 30–50 cm), we then individually tested for significant differences at each horizon using the non-parametric Kruskal-Wallis test with Dunn pair-wise comparisons (vegan package; Oksanen et al., 2015). Only Kruskal-Wallis tests where significant differences were encountered in one or more ice horizon are presented here.
To assess impacts of melt procedures on microbial community structure, a non-metric multidimensional scaling (NMDS) ordination was run based on the Bray-Curtis dissimilarity matrix (Bray and Curtis, 1957) of 16S/18S ASV read counts/relative abundances (phyloseq package; McMurdie and Holmes, 2013). A permutational analysis of variance (PERMANOVA) was then used to individually test for significant differences between melting methods, shaved and non-shaved ice, and ice horizons within each dissimilarity matrix (999 permutations). Confidence in results was confirmed using a beta dispersion test (vegan package; Oksanen et al., 2015). All p-values were adjusted to correct for multiple comparisons using the Holm-Bonferroni method with initial significance cut-off of α = 0.05 (Holm, 1979).
We used the package DeSeq2 to assess if specific 16S or 18S rRNA gene phylotypes were differentially abundant between ice-melting methods (Love et al., 2014). We limited this analysis to procedures that significantly impacted community structure in the PERMANOVA results. DESeq2 calculates dispersion estimates using a parametric, dispersion-mean relation through a robust gamma-family generalized linear model. Size factors were estimated using the median ratio method to control for differences in library size, and the Wald test was used to determine log2fold changes in ASV abundances. P-values for significant ASVs were adjusted for false discovery rate using the Benjamini-Hochberg adjustment where values below 0.1 were considered significant (the default setting for DESeq2; Love et al., 2014). Methods for utilizing DeSeq2 on ASVs were adapted from Webb et al. (2019) and Erazo and Bowman (2021).
Results
In this section, comparisons between melting method within individual ice horizons are referred to as method-5, method-20, and method-40, where method indicates direct melt (DM), seawater melt (SM) or isohaline melt (IM) and numbers indicate the midpoint for the 0–10 cm, 10–30 cm, and 30–50 cm ice horizons, respectively. For example, samples from the 0–10 ice horizon that were buffered using the seawater-salinity melt method are designated “SM-5” (Figure 2, Table 1).
Procedural overview with averages and standard deviations of key parameters
| Ice Horizon (cm) | Melting Method (label)a | Ice Shaving | AF Abundancea (cells mL−1) | TC Abundancea (cells mL−1) | Chlorophyll aa (mg L−1) | Effective Quantum Yield (%) |
|---|---|---|---|---|---|---|
| 0–10 | Direct melt (DM-5) | Shaved | 6.7 ± 4.0 × 104 | 2.2 ± 0.88 × 105 | 131 ± 17 | 0.09 ± 0.09 |
| Non-shaved | 6.7 ± 4.3 × 104 | 2.6 ± 0.26 × 105 | 151 ± 10 | 0.07 ± 0.01 | ||
| Isohaline melt (IM-5) | Shaved | 7.6 ± 1.4 × 104 | 3.1 ± 0.70 × 105 | 151 ± 42 | 0.11 ± 0.05 | |
| Non-shaved | 4.4 ± 0.80 × 104 | 2.3 ± 0.69 × 105 | 180 ± 332 | 0.04 ± 0.06 | ||
| Seawater melt (SM-5) | Shaved | 3.9 ± 1.3 × 104 | 2.6 ± 0.10 × 105 | 182 ± 22 | 0.15 ± 0.04 | |
| Non-shaved | 7.7 ± 3.6 × 104 | 6.5 ± 1.6 × 105 | 213 ± 40 | 0.04 ± 0.04 | ||
| 10–30 | Direct melt (DM-20) | Shaved | 2.4 ± 0.69 × 103 | 2.4 ± 0.78 × 104 | 12 ± 4 | 0.12 ± 0.07 |
| Non-shaved | 3.0 ± 0.37 × 103 | 2.7 ± 0.46 × 104 | 9 ± 2 | 0.09 ± 0.07 | ||
| Isohaline melt (IM-20) | Shaved | 2.5 ± 0.35 × 103 | 3.2 ± 0.62 × 104 | 11 ± 2 | 0.06 ± 0.03 | |
| Non-shaved | 2.6 ± 0.13 × 103 | 3.4 ± 1.0 × 104 | 18 ± 14 | 0.08 ± 0.08 | ||
| Seawater melt (SM-20) | Shaved | 3.9 ± 1.5 × 103 | 3.7 ± 0.85 × 104 | 25 ± 9 | 0.16 ± 0.04 | |
| Non-shaved | 2.7 ± 0.27 × 103 | 4.9 ± 0.99 × 104 | 19 ± 4 | 0.10 ± 0.08 | ||
| 30–50 | Direct melt (DM-40) | Shaved | 1.2 ± 0.40 × 103 | 1.8 ± 1.0 × 104 | 7 ± 1 | 0.12 ± 0.09 |
| Non-shaved | 2.0 ± 0.60 × 103 | 2.7 ± 0.49 × 104 | 14 ± 3 | 0.09 ± 0.05 | ||
| Isohaline melt (IM-40) | Shaved | 1.4 ± 0.68 × 103 | 1.4 ± 0.49 × 104 | 10 ± 4 | 0.12 ± 0.09 | |
| Non-shaved | 1.1 ± 0.57 × 103 | 2.2 ± 1.1 × 104 | 11 ± 2 | 0.10 ± 0.06 | ||
| Seawater melt (SM-40) | Shaved | 6.1 ± 2.4 × 103 | 6.6 ± 2.6 × 104 | 17 ± 4 | 0.11 ± 0.07 | |
| Non-shaved | 3.0 ± 0.52 × 103 | 3.6 ± 0.75 × 104 | 28 ± 16 | 0.08 ± 0.06 |
| Ice Horizon (cm) | Melting Method (label)a | Ice Shaving | AF Abundancea (cells mL−1) | TC Abundancea (cells mL−1) | Chlorophyll aa (mg L−1) | Effective Quantum Yield (%) |
|---|---|---|---|---|---|---|
| 0–10 | Direct melt (DM-5) | Shaved | 6.7 ± 4.0 × 104 | 2.2 ± 0.88 × 105 | 131 ± 17 | 0.09 ± 0.09 |
| Non-shaved | 6.7 ± 4.3 × 104 | 2.6 ± 0.26 × 105 | 151 ± 10 | 0.07 ± 0.01 | ||
| Isohaline melt (IM-5) | Shaved | 7.6 ± 1.4 × 104 | 3.1 ± 0.70 × 105 | 151 ± 42 | 0.11 ± 0.05 | |
| Non-shaved | 4.4 ± 0.80 × 104 | 2.3 ± 0.69 × 105 | 180 ± 332 | 0.04 ± 0.06 | ||
| Seawater melt (SM-5) | Shaved | 3.9 ± 1.3 × 104 | 2.6 ± 0.10 × 105 | 182 ± 22 | 0.15 ± 0.04 | |
| Non-shaved | 7.7 ± 3.6 × 104 | 6.5 ± 1.6 × 105 | 213 ± 40 | 0.04 ± 0.04 | ||
| 10–30 | Direct melt (DM-20) | Shaved | 2.4 ± 0.69 × 103 | 2.4 ± 0.78 × 104 | 12 ± 4 | 0.12 ± 0.07 |
| Non-shaved | 3.0 ± 0.37 × 103 | 2.7 ± 0.46 × 104 | 9 ± 2 | 0.09 ± 0.07 | ||
| Isohaline melt (IM-20) | Shaved | 2.5 ± 0.35 × 103 | 3.2 ± 0.62 × 104 | 11 ± 2 | 0.06 ± 0.03 | |
| Non-shaved | 2.6 ± 0.13 × 103 | 3.4 ± 1.0 × 104 | 18 ± 14 | 0.08 ± 0.08 | ||
| Seawater melt (SM-20) | Shaved | 3.9 ± 1.5 × 103 | 3.7 ± 0.85 × 104 | 25 ± 9 | 0.16 ± 0.04 | |
| Non-shaved | 2.7 ± 0.27 × 103 | 4.9 ± 0.99 × 104 | 19 ± 4 | 0.10 ± 0.08 | ||
| 30–50 | Direct melt (DM-40) | Shaved | 1.2 ± 0.40 × 103 | 1.8 ± 1.0 × 104 | 7 ± 1 | 0.12 ± 0.09 |
| Non-shaved | 2.0 ± 0.60 × 103 | 2.7 ± 0.49 × 104 | 14 ± 3 | 0.09 ± 0.05 | ||
| Isohaline melt (IM-40) | Shaved | 1.4 ± 0.68 × 103 | 1.4 ± 0.49 × 104 | 10 ± 4 | 0.12 ± 0.09 | |
| Non-shaved | 1.1 ± 0.57 × 103 | 2.2 ± 1.1 × 104 | 11 ± 2 | 0.10 ± 0.06 | ||
| Seawater melt (SM-40) | Shaved | 6.1 ± 2.4 × 103 | 6.6 ± 2.6 × 104 | 17 ± 4 | 0.11 ± 0.07 | |
| Non-shaved | 3.0 ± 0.52 × 103 | 3.6 ± 0.75 × 104 | 28 ± 16 | 0.08 ± 0.06 |
aValues are the average of meltwater measurements (n = 3) multiplied by the dilution factor for the buffered methods, which for seawater melt (SM) samples was 2.06 and for isohaline melt (IM) samples was 1.2 (IM-5), 1.3 (IM-20), and 1.6 (IM-40). AF indicates autofluorescent cells; TC, total cells.
Cell abundance and photophysiology
Average total cell (TC) abundance and autofluorescent cell (AF) abundance did not vary significantly between shaved and non-shaved ice horizons (AF, p = 0.797; TC, p = 0.167; Table 2). However, average TC and AF abundance did vary significantly between melting methods (ANOVA test: TC, p < 0.0001; AF, p = 0.0345; Table 2), and Kruskal-Wallis analyses for each ice horizon revealed that melting method was a statistically significant factor only in the two uppermost sections of the ice column (Figure 3A and B, Table 3). At 10–30 cm, the SM method resulted in significantly higher TC abundances than the DM method (Figure 3A, Table 3), with mean TC abundance in the SM-20 samples 75% higher than DM-20 samples (4.3 × 104 and 2.6 × 104 cells mL−1, respectively; Table 1).
ANOVA summary statistics
| Melting Method | Ice Shaving | Interaction | ||||
|---|---|---|---|---|---|---|
| Parametera | F (2, xx)b | P adjc | F (1, xx)b | P adjc | F (2, xx)b | P adjc |
| AF abundance (xx = 48) | 4.912 | 0.0345 | 0.067 | 0.7965 | 1.675 | 0.3960 |
| TC abundance (xx = 48) | 13.73 | <0.0001 | 3.122 | 0.1672 | 0.550 | 0.5808 |
| Chlorophyll a (xx = 46) | 16.6 | <0.0001 | 5.179 | 0.0552 | 0.460 | 0.6345 |
| Effective quantum yield (xx = 52) | 0.112 | 1.000 | 3.878 | 0.1629 | 0.165 | 1.000 |
| Bacterial diversity (xx = 44) | 0.555 | 1.000 | 0.000 | 1.000 | 0.295 | 1.000 |
| Eukaryotic diversity (xx = 51) | 3.489 | 0.114 | 2.345 | 0.264 | 0.975 | 0.384 |
| Melting Method | Ice Shaving | Interaction | ||||
|---|---|---|---|---|---|---|
| Parametera | F (2, xx)b | P adjc | F (1, xx)b | P adjc | F (2, xx)b | P adjc |
| AF abundance (xx = 48) | 4.912 | 0.0345 | 0.067 | 0.7965 | 1.675 | 0.3960 |
| TC abundance (xx = 48) | 13.73 | <0.0001 | 3.122 | 0.1672 | 0.550 | 0.5808 |
| Chlorophyll a (xx = 46) | 16.6 | <0.0001 | 5.179 | 0.0552 | 0.460 | 0.6345 |
| Effective quantum yield (xx = 52) | 0.112 | 1.000 | 3.878 | 0.1629 | 0.165 | 1.000 |
| Bacterial diversity (xx = 44) | 0.555 | 1.000 | 0.000 | 1.000 | 0.295 | 1.000 |
| Eukaryotic diversity (xx = 51) | 3.489 | 0.114 | 2.345 | 0.264 | 0.975 | 0.384 |
aAF indicates autofluorescent cells; TC, total cells; parenthetic value, degrees of freedom.
bTwo-way crossed analysis of variance (ANOVA), testing melting method (direct, isohaline, and seawater) and ice shaving (shaved, non-shaved) with an error factor to account for ice horizon.
cAll p-values adjusted to account for multiple comparisons using the Holm-Bonferroni correction, α = 0.05 (Holm, 1979).
Microbial parameters by ice horizon, shaving, and melting method. (A) Total cell (TC) abundance, (B) autofluorescent (AF) cell abundance, (C) chlorophyll a concentration, (D) effective quantum yield, (E) bacterial alpha diversity, and (F) eukaryotic alpha diversity are plotted by ice horizon and melting procedure, with individual measurements represented by triangles (shaved) or circles (non-shaved) and colored by melting method (direct, isohaline, or seawater melt). The midline of each box represents the median, below which is the first quartile and above, the third quartile. Whiskers represent the minimum and maximum of the spread, and outliers are marked by standalone points. Cell abundances and chlorophyll a concentrations were log-transformed, and alpha diversity was calculated using the inverse Simpson Index. Significance values (ns = p > 0.05, * = p < 0.05, ** = p < 0.01) represent a one-way non-parametric Kruskal-Wallis tests run for melting method on each individual ice horizon; p-values were adjusted to account for multiple comparisons using the Holm-Bonferroni correction, α = 0.05 (Holm, 1979).
Microbial parameters by ice horizon, shaving, and melting method. (A) Total cell (TC) abundance, (B) autofluorescent (AF) cell abundance, (C) chlorophyll a concentration, (D) effective quantum yield, (E) bacterial alpha diversity, and (F) eukaryotic alpha diversity are plotted by ice horizon and melting procedure, with individual measurements represented by triangles (shaved) or circles (non-shaved) and colored by melting method (direct, isohaline, or seawater melt). The midline of each box represents the median, below which is the first quartile and above, the third quartile. Whiskers represent the minimum and maximum of the spread, and outliers are marked by standalone points. Cell abundances and chlorophyll a concentrations were log-transformed, and alpha diversity was calculated using the inverse Simpson Index. Significance values (ns = p > 0.05, * = p < 0.05, ** = p < 0.01) represent a one-way non-parametric Kruskal-Wallis tests run for melting method on each individual ice horizon; p-values were adjusted to account for multiple comparisons using the Holm-Bonferroni correction, α = 0.05 (Holm, 1979).
Kruskal-Wallis summary statistics for pairwise comparisons of melting method
| Ice Horizon | ||||
|---|---|---|---|---|
| Parameter | Test Statistics | 0–10 cm (xx = 26, yy = 19, zz = 22) | 10–30 cm (xx = 24, yy = 17, zz = 19) | 30–50 cm (xx = 18, yy = 18, zz = 18) |
| AF abundance | KWa, F (2, xx) | 0.573 | 0.986 | 11.1 |
| KWa, P | 0.751 | 0.611 | 0.00389 | |
| DNb, P adj <0.05 | NSc | NS | 0.0221, SM–DM; 0.00513, SM–IM | |
| TC abundance | KWa, F (2, xx) | 3.52 | 9.75 | 10.7 |
| KWa, P | 0.172 | 0.00765 | 0.00472 | |
| DNb, P adj <0.05 | NS | 0.00683, SM–DM | 0.0347, SM–DM; 0.00513, SM–IM | |
| Chlorophyll a | KWa, F (2, xx) | 7.03 | 9.06 | 10.2 |
| KWa, P | 0.0298 | 0.0108 | 0.00624 | |
| DNb, P adj <0.05 | 0.0242, SM–DM | 0.00882, SM–DM | 0.0148, SM–DM; 0.0148, SM–IM | |
| Bacterial diversity | KWa, F (2, yy) | 1.04 | 6.94 | 4.46 |
| KWa, P | 0.594 | 0.0311 | 0.107 | |
| DNb, P adj <0.05 | NS | 0.0265, SM–IM | NS | |
| Eukaryotic diversity | KWa, F (2, zz) | 0.708 | 5.04 | 11.4 |
| KWa, P | 0.702 | 0.080 | 0.00332 | |
| DNb, P adj <0.05 | NS | NS | 0.00738, SM–DM; 0.00985, SM–IM | |
| Ice Horizon | ||||
|---|---|---|---|---|
| Parameter | Test Statistics | 0–10 cm (xx = 26, yy = 19, zz = 22) | 10–30 cm (xx = 24, yy = 17, zz = 19) | 30–50 cm (xx = 18, yy = 18, zz = 18) |
| AF abundance | KWa, F (2, xx) | 0.573 | 0.986 | 11.1 |
| KWa, P | 0.751 | 0.611 | 0.00389 | |
| DNb, P adj <0.05 | NSc | NS | 0.0221, SM–DM; 0.00513, SM–IM | |
| TC abundance | KWa, F (2, xx) | 3.52 | 9.75 | 10.7 |
| KWa, P | 0.172 | 0.00765 | 0.00472 | |
| DNb, P adj <0.05 | NS | 0.00683, SM–DM | 0.0347, SM–DM; 0.00513, SM–IM | |
| Chlorophyll a | KWa, F (2, xx) | 7.03 | 9.06 | 10.2 |
| KWa, P | 0.0298 | 0.0108 | 0.00624 | |
| DNb, P adj <0.05 | 0.0242, SM–DM | 0.00882, SM–DM | 0.0148, SM–DM; 0.0148, SM–IM | |
| Bacterial diversity | KWa, F (2, yy) | 1.04 | 6.94 | 4.46 |
| KWa, P | 0.594 | 0.0311 | 0.107 | |
| DNb, P adj <0.05 | NS | 0.0265, SM–IM | NS | |
| Eukaryotic diversity | KWa, F (2, zz) | 0.708 | 5.04 | 11.4 |
| KWa, P | 0.702 | 0.080 | 0.00332 | |
| DNb, P adj <0.05 | NS | NS | 0.00738, SM–DM; 0.00985, SM–IM | |
aOne-way, nonparametric Kruskal-Wallis (KW) tests run with melting method factors individually for each ice horizon. DM indicates direct melt; IM, isohaline melt; SM, seawater melt.
bP-values from Dunn pairwise comparisons (DN) adjusted to account for multiple comparisons using the Holm-Bonferroni correction, α = 0.05 (Holm, 1979). Only significant values reported.
cNS indicates no significant comparisons between melting methods.
At 30–50 cm, the SM method resulted in significantly higher TC and AF abundances than either the DM or IM methods (Figure 3A and B, Table 3). Mean AF abundance in SM-40 samples was 188% higher than DM-40 and 260% higher than IM-40 samples (Table 1). Mean TC abundance in SM-40 samples was 181% higher than IM-40 samples and 128% higher than DM-40 samples (Table 1).
Based on the two-way ANOVA analysis, mean chlorophyll a concentration also varied significantly by melting method, but not by ice shaving (p < 0.0001 and p = 0.0552, respectively; Table 2). SM samples retained a significantly higher mean chlorophyll a concentration than DM samples across all ice horizons (41% in SM-5, 100% in SM-20, and 109% in SM-40; Figure 3C, Tables 1 and 3). Mean chlorophyll a concentration was also significantly higher (109%) in mean SM-40 samples than IM-40 samples (Figure 3C, Tables 1 and 3). Effective quantum yield showed no significant variation between the shave or melt treatments (Figure 3D, Table 2).
Microbial community structure
Very few archaeal phylotypes (1.6% of 16S ASVs) were detected following sequencing quality control, so the prokaryotic analysis of this study focuses solely on the bacterial community. Initially, 42% of the bacterial community was identified as cyanobacteria/chloroplasts and 22% as mitochondrial DNA. After removing these ASVs, bacterial community structure was dominated by Proteobacteria (28% Gamma, 22% Alpha, 0.6% Beta) and Flavobacteriia (33%) taxa across all samples (including all melt procedures and ice horizons, shaved and non-shaved; Figure 4A). The parametric two-way crossed ANOVA analysis across all samples yielded no significant variation in bacterial alpha diversity among 16S ASVs within either melting methods or ice shaving (Table 2). However, when differences were assessed individually within ice core sections, the non-parametric Kruskal-Wallis and Dunn Pairwise Comparison showed significantly higher alpha diversity in SM-20 than in IM-20 (Figure 3E, Table 3). These differing results suggest that while there was significant variation between melting method in the 10–30 cm ice horizon, this trend is lost to noise when analyzing all samples in the dataset. While melting methods did not separate visually in the first two dimensions of the NMDS ordination (Figure S1A), the PERMANOVA on the Bray-Curtis dissimilarity matrix of 16S ASV read counts revealed that dispersion among groups was significantly different across both melting method (p = 0.012) and ice horizon (p = 0.004; Table 4). This result indicates strong methodological variation in sequence-based bacterial community structure from melting approach. Confidence in PERMANOVA results was confirmed by the beta dispersion test (F = 0.662, Df = 2, p = 0.516).
Overview of taxonomic diversity across all samples. Percentage of unique amplicon sequence variants (ASVs) assigned as a completed, or most closely related to, genome within each taxonomic class across all samples (all treatments and ice horizons). (A) Bacterial class percentages calculated using read counts of each 16S ASV normalized to cell count. (B) Eukaryotic class percentages calculated using relative abundance of each 18S ASV.
Overview of taxonomic diversity across all samples. Percentage of unique amplicon sequence variants (ASVs) assigned as a completed, or most closely related to, genome within each taxonomic class across all samples (all treatments and ice horizons). (A) Bacterial class percentages calculated using read counts of each 16S ASV normalized to cell count. (B) Eukaryotic class percentages calculated using relative abundance of each 18S ASV.
PERMANOVA summary statistics
| Melting Method | Ice Shaving | Ice Horizon | Interaction 1:2 | |||||
|---|---|---|---|---|---|---|---|---|
| Parameter | F (2, xx)a | P adjb | F (1, xx)a | P adjb | F (2, xx)a | P adjb | F (2, xx)a | P adjb |
| 16S ASV dissimilarity matrix (44) | 4.266 | 0.012 | 3.016 | 0.066 | 46.94 | 0.004 | 1.167 | 0.304 |
| 18S ASV dissimilarity matrix (51) | 2.438 | 0.051 | 2.057 | 0.364 | 24.83 | 0.004 | 0.940 | 0.298 |
| Melting Method | Ice Shaving | Ice Horizon | Interaction 1:2 | |||||
|---|---|---|---|---|---|---|---|---|
| Parameter | F (2, xx)a | P adjb | F (1, xx)a | P adjb | F (2, xx)a | P adjb | F (2, xx)a | P adjb |
| 16S ASV dissimilarity matrix (44) | 4.266 | 0.012 | 3.016 | 0.066 | 46.94 | 0.004 | 1.167 | 0.304 |
| 18S ASV dissimilarity matrix (51) | 2.438 | 0.051 | 2.057 | 0.364 | 24.83 | 0.004 | 0.940 | 0.298 |
aPermutational multivariate analysis of variance performed on Bray-Curtis dissimilarity matrices of 16S amplicon sequence variant (ASV) read counts normalized to total cell count and of eukaryotic 18S ASV relative abundances.
bP-values adjusted using the Holm-Bonferroni correction to account for multiple comparisons, α = 0.05 (Holm, 1979).
Eukaryotic community composition across all samples was dominated by diatoms (29% Bacillariophyta), amoeboid Cercozoa (18% Filosa-Imbricatea, 10% Filosa-Thecofilosea), and Preaxostyla (17%; Figure 4B). Similar to the bacterial community, there was no significant difference in eukaryotic alpha diversity among 18S relative abundances within the melting methods or ice shaving in the two-way crossed ANOVA (Table 2). Yet, the non-parametric Kruskal-Wallis and Dunn Pairwise Comparison showed significantly lower alpha diversity in SM-40 than either DM-40 or IM-40, again indicating that while some significant variation does occur from melting method within ice horizons (Figure 3F, Table 3), it is not strong enough for significance when all samples are considered. Clusters by melting method were not identifiable visually in the NMDS ordination (Figure S1B), and the PERMANOVA on the Bray-Curtis dissimilarity matrix of transformed 18S ASV relative abundances identified no significant dispersion among groups for melting method (p = 0.051), only ice horizon (p = 0.004; Table 4). This result indicates little methodological variation in sequence-based eukaryotic community structure from melting approach. Confidence in PERMANOVA results was confirmed by the beta dispersion test (F = 0.5035, Df = 2, p = 0.62).
Ice shaving was not a significant contributor to dispersion within the Bray-Curtis dissimilarity matrix of either bacterial ASV read counts or eukaryotic ASV relative abundances (Table 4) and were not considered in further analyses.
Differential abundances within bacterial taxa
Our DeSeq2 analysis yielded many 16S phylotypes that were differentially abundant between melting methods. Comparisons between each method (SM–DM, IM–DM, and IM–SM) were run independently for each ice horizon, where positive log2-fold changes indicate ASVs that were significantly (p < 0.1) more abundant in the first of the two compared methods and negative log2-fold changes indicate the opposite (Figure S2). A total of 55 ASVs were identified as significantly differentially abundant when comparing the SM and DM methods. Only six of these ASVs were more abundant in DM samples (two belonging to class Spirochaetia, two to Flavobacteriia, one to Alphaproteobacteria, and one to Gammaproteobacteria; Figure S2A and C), while 49 were more abundant in SM samples (26% Gammaproteobacteria, 21% Alphaproteobacteria, 12% Flavobacteriia, and the remaining 41% split among 10 other classes; Figure S2A and C). When comparing IM and SM methods, 71 ASVs were identified as significantly differentially abundant, with 15 being more abundant in SM samples (27% Flavobacteriia, 20% Alphaproteobacteria, 13% Acidimicrobiia, and 40% among other classes; Figure S2D and F) and 56 more abundant in IM samples (25% Gammaproteobacteria, 20% Flavobacteriia, 14% Alphaproteobacteria, 5% unidentified bacteria, and 36% among other classes; Figure S2D–F). When comparing IM and DM methods, 90 ASVs were identified as significantly differentially abundant, with 4 being more abundant in DM samples (one each of Alphaproteobacteria, Flavobacteriia, Spirochaetia, and Mollicutes; Figure S2G and I), and 86 more abundant in IM samples (27% Gammaproteobacteria, 20% Alphaproteobacteria, 14% Flavobacteriia, 5% unidentified bacteria, and 33% among other classes; Figure S2G–I).
When considering all ice horizon and melting method comparisons (removing duplicate ASVs, meaning that the ASV was significant in more than one of the above comparisons), a total of 110 unique bacterial ASVs were identified as significantly differentially abundant (Table S2). Figure 5 presents the 50 most abundant of these ASVs, or those with the highest read counts across all samples. These data indicate that the strongest variations in read count occurred in the 0–10 cm ice horizon (Figure 5). Many of the significantly differentially abundant ASVs with high read counts across samples appeared more abundant in the buffered SM and/or IM samples than in the DM samples (Figure 5). Most of these ASVs had high placement proportions to reference genomes within Proteobacteria such as Candidatus Pelagibacter or Bacteroidetes such as Owenweeksia hongkongensis (DSM 17368), Fluviicola taffensis (DSM16823) or other unidentified Flavobacteriaceae (Figure 5). Some ASVs appeared significantly more abundant in DM samples, placing mostly to a reference genome belonging to Glaciecola amylolytica (Figure 5). Significant differential abundances between IM-5 and SM-5 (Figure S2F) were among less abundant ASVs or among ASVs placed to the same reference genome (Table S2).
Significantly differentially abundant bacterial amplicon sequence variants (ASVs) between melting method within each ice horizon. Cell count normalized 16S ASV read counts were scaled using a z-score between 0 and 1. Bray-Curtis dissimilarity distances were used to cluster taxa. Columns (samples) are colored by ice horizon and melting method. The 50 most abundant of significantly (p < 0.1) differentially abundant ASVs from each comparison are shown with duplicate ASVs (significant in more than one comparison) shown only once. The paprica taxonomic ID (closest relative among RefSeq or PR2 genomes) and placement proportion were used to identify each unique ASV (Bowman and Ducklow, 2015). Heatmaps were generated using R package pheatmap (Kolde, 2019).
Significantly differentially abundant bacterial amplicon sequence variants (ASVs) between melting method within each ice horizon. Cell count normalized 16S ASV read counts were scaled using a z-score between 0 and 1. Bray-Curtis dissimilarity distances were used to cluster taxa. Columns (samples) are colored by ice horizon and melting method. The 50 most abundant of significantly (p < 0.1) differentially abundant ASVs from each comparison are shown with duplicate ASVs (significant in more than one comparison) shown only once. The paprica taxonomic ID (closest relative among RefSeq or PR2 genomes) and placement proportion were used to identify each unique ASV (Bowman and Ducklow, 2015). Heatmaps were generated using R package pheatmap (Kolde, 2019).
Differential abundances within eukaryotic taxa
Several 18S phylotypes yielded differential abundances between melting methods, with a large proportion of those more abundant in the buffered melts (SM or IM) belonging to classes within the division Cercozoa (i.e., Filosa-Imbricatea, Filosa-Thecofilosea, Filosa-Sarcomonadea, and Filosa-Granofilosea; Figure S3). A total of 24 eukaryotic ASVs were identified as significantly differentially abundant when comparing the SM and DM methods, with 14 being more abundant in DM samples (36% Bacillariophyta, 36% Dinophyceae, and 28% other) and 10 more abundant in SM samples (50% Cercozoa, 30% Bacillariophyta, 20% other; Figure S3A–C). When comparing IM and SM methods, 39 eukaryotic ASVs were identified as significantly differentially abundant, with 16 being more abundant in IM samples (varying classes) and 23 more abundant in SM samples (34% Bacillariophyta, 22% Cercozoa, 17% other; Figure S3D–F). When comparing IM and DM methods, 80 eukaryotic ASVs were identified as significantly differentially abundant, with 38 more abundant in IM samples (13% Cercozoa, and other varying taxa) and 21 more abundant in DM samples (38% Dinophyceae, 33% Bacillariophyta, and 28% other; Figure S3G–I).
When considering all ice horizon and melting method comparisons (again removing duplicates), a total of 81 unique 18S ASVs were significantly differentially abundant (Table S3). Figure 6 presents the 50 most abundant of these significant eukaryotic ASVs across all samples. A relatively large proportion of these ASVs placed to reference genomes that belonged to the division Cercozoa (i.e., Protaspa-lineage sp., Novel clade-2 X sp., CCW10-lineage X sp., Spongomonas minima, Marimonadida XX sp., Cryothecomonas sp., Cercomonas sp., Amphiprora paludosa, Trinematidae, and Capsellina sp.). However, distinctions were more pronounced among ice horizons than melting methods, with subsets of ASVs relatively more abundant in the 0–10 cm or 30–50 cm horizon and some shared ASVs of both in the 10–30 cm horizon (Figure 6).
Significantly differentially abundant eukaryotic amplicon sequence variants (ASVs) between melting method within each ice horizon. Relative 18S ASV abundances were scaled using a z-score between 0 and 1. Bray-Curtis dissimilarity distances were used to cluster taxa. Columns (samples) are colored by ice horizon and melting method. The 50 most abundant significantly (p < 0.1) differentially abundant ASVs from each comparison are shown with duplicate ASVs (significant in more than one comparison) shown only once. The paprica taxonomic ID (closest relative among RefSeq or PR2 genomes) and placement proportion were used to identify each unique ASV (Bowman and Ducklow, 2015). Heatmaps were generated using R package pheatmap (Kolde, 2019).
Significantly differentially abundant eukaryotic amplicon sequence variants (ASVs) between melting method within each ice horizon. Relative 18S ASV abundances were scaled using a z-score between 0 and 1. Bray-Curtis dissimilarity distances were used to cluster taxa. Columns (samples) are colored by ice horizon and melting method. The 50 most abundant significantly (p < 0.1) differentially abundant ASVs from each comparison are shown with duplicate ASVs (significant in more than one comparison) shown only once. The paprica taxonomic ID (closest relative among RefSeq or PR2 genomes) and placement proportion were used to identify each unique ASV (Bowman and Ducklow, 2015). Heatmaps were generated using R package pheatmap (Kolde, 2019).
Discussion
Our results agree with previous studies that sea ice melting methodology has an impact on resulting microbial parameters (Garrison and Buck, 1986; Mikkelsen and Witkowski, 2010; Campbell et al., 2019; Roukaerts et al., 2019). They advance understanding of this impact by providing new insights into the effects of melting on bacterial and eukaryotic community structure at the 16S and 18S rRNA gene level.
Impacts of melting method were greatest on measurements of cell abundance and chlorophyll a content, where the SM method resulted in significantly higher TC and AF cell abundances (upper ice only) and chlorophyll a concentration (all ice horizons) than either the IM or DM methods (Figure 3). Significant patterns were consistent between the AF and TC abundance measurements at 30–50 cm, indicating that the SM method buffers equally well for both the photosynthetic and non-photosynthetic communities. However, the photosynthetic community (Figure 3B) may be somewhat less sensitive to osmotic stress than the broader microbial community closer to the ice-water interface (SM-5 samples not visually higher and SM-20 samples not significantly higher as observed for TC abundances in Figure 3A). When taken together, the greater cell and chlorophyll retention in the SM method of our study agrees well with the results of Campbell et al. (2019), who found that large volume dilutions (8:1 filtered seawater to ice) resulted in highest values of photosynthetic parameters and production in bottom first year sea ice. Both our study and that by Campbell et al. (2019) were based on sea ice collected during Arctic spring and indicate that large volume dilutions with buffers at or close to seawater salinities (29–35 ppt) result in “best” results for the parameters measured. Additionally, because our study used a sterile NaCl solution as opposed to filtered natural seawater, this conclusion is reached without the potential bias that a natural seawater buffer might bring with its addition of dissolved organic carbon and nutrients to the sample.
However, a wintertime study that included measurements of bacterial abundance in land-fast first year sea ice at a site similar to ours (near Utquiaġvik, AK) yielded somewhat contrasting results. Ewert et al. (2013) found that direct melts resulted in significant (55%) cell loss only in the uppermost 3 cm of this colder ice (approximately −12°C). While their results agree with ours that sea ice communities living closer to the ice-seawater interface appear to be less prone to cell loss, they observed higher bacterial abundances in this upper ice using a near isohaline buffer method, melting the samples with brine salinities of approximately 160 ppt to a meltwater salinity of 100 ppt (Ewert et al., 2013). In our study, the IM buffer method did not result in significantly higher cell counts for any ice horizon (Figure 3A), including the coldest upper ice we collected (−5.6°C), and the SM buffer treatment outperformed both IM and DM methods for cell retention in this upper sea ice (Figure 3A and B). This performance likely stems from the season and associated differences in the in situ ice conditions.
Our samples were collected during a relatively warm period when the air temperature at the time of sampling (−11°C) was above the average April high (−13°C) according to the US Climate database. The underlying seawater temperature was not assessed here, but the ice temperatures for all cores were near or above −5°C. Constrained by a two-phase equilibrium, the dominant processes responsible for brine channel formation are heat and mass transport (Worster and Wettlaufer, 1997), and some brine loss or ice-ocean exchange may have been occurring, as previously observed during this season (Niedrauer and Martin, 1979; Buck et al., 1998; Morawetz et al., 2017). If so, the high eukaryotic diversity seen in the 10–30 cm ice horizon could be from infiltration or cell movement, preventing the dominance of a single species (Stoecker et al., 1992; van Leeuwe et al., 2018). This potential seawater infiltration, in addition to the relatively large size (10–20 cm length) and warm temperatures (Table S1) of our core sections, could potentially explain the outperformance of the larger volume SM method for our samples. The IM method may have been more osmotically stressful than helpful due to inadequate buffering by the low volumes of hypersaline brine required to match the salinity conditions within our sampled ice. Either way, considering the in situ ice conditions during sampling and the total ice volume of subsamples is important when considering the melting method to apply (Miller et al., 2015).
In addition to total cell abundances and photophysiology, our results suggest that melting method coupled to differential osmotic cell lysis has the potential to impact estimates of bacterial community structure from 16S rRNA gene sequencing. With several taxa significantly more abundant in the buffered SM and IM methods than the unbuffered DM method, where only a few taxa are significantly more abundant in comparison to the buffered methods, our analyses indicate that buffered melts retain a broader spectrum of the bacterial community. However, unlike cell count and chlorophyll a concentration, neither buffered method (SM or IM) appears particularly better than the other at ASV retention for our sample set. Differences in ASV read count between the IM and SM comparisons were less striking than the differences between the IM and DM or SM and DM comparisons (Figure 5), which accounted for many of the ASVs that were significantly differentially abundant across multiple comparisons (Table S2). While direct measurements of alpha diversity were only significantly impacted in the 10–30 cm layer (Figure 3E), dispersion among groups for all samples was significantly different between melting methods (Table 4), indicating the potential for community bias throughout the ice column. Indeed, the DeSeq2 analysis resulted in more ASVs that were differentially abundant in the 0–10 cm core section than in the upper ice sections for each melt treatment comparison (Figure S2). This finding contrasts with our cell abundance results, indicating that, while cell loss is reduced at the ice-seawater interface, either the larger community size or stronger seawater influence results in a community subset that is less well adapted to osmotic stress and therefore significantly less abundant using the DM method. This variability in bacterial community structure across melt treatments indicates that cell losses between buffered and unbuffered methods can result in biased abundances of specific subsets of the bacterial community.
Most of the 110 significant ASVs (Table S2) were identified as belonging to classes that comprise the bulk of the community found in our samples (Figure 4A) and commonly found in sea ice environments, such as Alpha- and Gammaproteobacteria and members of the Bacteroidetes (Boetius et al., 2015). This finding reduces concern that any differential abundances in phylotype were caused by a potential contamination of the sterile NaCl solution during transit to the field. Although the taxonomic classifications used in this study are not exact and merely represent the closest phylogenetic neighbor based on complete genomes in reference databases, they provide insight into the ecological functions that may be affected by preferential cell lysis due to melting method. Many of the ASVs significantly more abundant in either the IM or SM methods were identified by paprica as members of the ubiquitous Pelagibacterales order (i.e., Candidatus Pelagibacter), known for genomic streamlining and their ability to survive across low nutrient regimes (Giovannoni et al., 2005), or members of the diverse phylum Bacteroidetes associated with processing complex organic compounds and byproducts of primary production (e.g., Flavobacteriaceae such as Owenweeksia hongkongensis DSM 17368 or Cryomorphaceae such as Fluviicola taffensis DSM 16823; Hahnke et al., 2016; Figure 5). Why these specific subsets of the microbial community may have been particularly affected by buffering method is unclear, but in studies based on direct melting, the biased removal of such cells from bacterial community structure may have ramifications for various ecological measurements (carbon or nutrient turnover, etc.) and interpretations when comparing bacterial community composition across different sea ice environments. Other ASVs, identified by paprica as belonging to Glaciecola amylolytica of the Gammaproteobacteria, were significantly more abundant in the DM samples (Figure 5). Originally isolated from Antarctic sea ice (Bowman et al., 1998), Glaciecola is an aerobic, halophilic genus since found in many marine environments (Xiao et al., 2019); perhaps its sea ice members are more suited to surviving the osmotic stress associated with steep salinity shifts.
Contrary to previous studies (Garrison and Buck, 1986; Mikkelsen and Witkowski, 2010), the melting methods we tested did not appear to bias eukaryotic microbial community structure, as dispersion among groups in our 18S ASVs did not differ significantly among the methods. Nevertheless, the results from our DeSeq2 analysis showed a total of 81 unique ASVs as significantly differentially abundant between melting methods (Table S3). In all but the IM–SM melt comparison, the highest number of differentially abundant ASVs were again found in the 0–10 cm core section at the ice-seawater interface (Figure S3). The IM–DM comparison resulted in the greatest number of significantly differential ASVs, but many of these were shared with the SM–DM comparison (Table S3). While the taxonomic spread of significant ASVs was generally representative of the total community composition (Figure 4), one trend was that members of the division Cercozoa, small single-celled protozoan eukaryotes, were more likely to be significantly higher in abundance in the buffered melts than the direct melt. This enrichment could, again, be due to the ice conditions at time of sampling, as well as potentially analysis technique. The relative abundance of protozoa in Arctic sea ice, such as Cercozoa species, is highest in early Spring (van Leeuwe et al., 2018). As melt begins, the diatom community at the ice-seawater interface is the first to “slough off,” while flagellate species endure longer in bottom sea ice communities (Tamelander et al., 2009; Torstensson et al., 2015; van Leeuwe et al., 2018). Previous melt studies have used microscopic analyses to conduct taxonomic assignment and define biases among osmotic buffers and direct melt methods, yet numbers of heterotrophic protozoa are often underestimated in field microscopic analyses because they lack chlorophyll and are small (van Leeuwe et al., 2018). By taking a sequencing approach, the inclusion of more of these protists may have complicated the trends seen in previous studies.
Across all analyses, mechanically shaving the sea ice before melting it had no direct impact on our downstream measurements. However, this procedure added to the overall sample processing time at room temperature, increasing the potential for brine loss prior to melting. While ice-shaving did shorten melting time by a few hours for the buffered melts (shaved DM samples melted more slowly due to snow-like insulation), our core sections were large enough that the degree of difference was smaller than previously determined using thinner ice sections (Campbell et al., 2019).
Conclusions
Based on these results, we conclude that variations in melting methods for sea ice have a significant impact on resulting measurements of bacterial and eukaryotic community structure, while first shaving the ice does not. The degree of bias by melting method depends on the type of measurement pursued and depth horizon of the ice. For spring-time collections of coastal first year sea ice, our results suggest that large (equivalent) volume dilutions with buffers at or near local seawater salinity best minimize osmotic stress and retain the most cell material. However, use of either of the osmotic NaCl buffers we tested—isohaline or seawater salinity—would be sufficient to reduce biased taxonomic representation in microbial community composition relative to directly melted ice. While cell abundance varied across melting method most strongly in the upper ice (10–30 cm and 30–50 cm core sections), the potential for community level bias at the ice-seawater interface (0–10 cm core section) appears to be higher, particularly for subsets of the bacterial community within the Proteobacteria and Bacteroidetes taxonomic groups. Eukaryotic community level biases among the melting methods were less apparent using our sequencing approach than in previously reported microscope-based work. Ultimately, these results confirm that the “best” melt procedure should be chosen based on research goals, sea ice sampling conditions and final analysis techniques. Repeating similar experiments across varying ice types, seasons, and microbial community compositions would help to identify more conditions in which current methodological practices bias measurements in sea ice and potentially the ensuing biogeochemical and ecological interpretations of this important changing polar habitat.
Data accessibility statement
All data are available on NCBI SRA at BioProject PRJNA779070: https://www.ncbi.nlm.nih.gov/bioproject/779070.
Supplemental files
The supplemental files for this article can be found as follows:
Figure S1. Impacts of melting method and ice shaving on total microbial community structure. Individual samples presented in the first two dimensions of a non-metric multi-dimensional scaling (NMDS) ordination created from (A) bacterial (16S) ASV read counts normalized to cell count and (B) eukaryotic (18S) ASV relative abundances. Samples are distinguished by melt method (color) and whether or not they were shaved prior to melt (shape). (Figure S1.png)
Figure S2. Ratio intensity plots of differentially abundant bacterial ASVs between melt method and ice horizon. The log2-fold change in 16S amplicon sequence variant (ASV) read counts normalized to total cell count between (A) DM-40 and SM-40, (B) DM-20 and SM-20, (C) DM-5 and SM-5, (D) IM-40 and SM-40, (E) IM-20 and SM-20, (F) IM-5 and SM-5, (G) IM-40 and DM-40, (H) IM-20 and DM-20, and (I) IM-5 and DM-5. The labeling system includes melt method (SM for seawater melt, IM for isohaline melt, and DM for direct melt) and ice horizon (5, 20, and 40 for 0–10 cm, 10–30 cm, and 30–50 cm, respectively). When comparing two melt methods, a positive log2-fold change represents higher abundance in the first method listed and a negative log2-fold change indicates higher abundance in the second method. Significantly differential abundances (p < 0.1) are displayed in color by taxonomic class. (Figure S2.png)
Figure S3. Ratio intensity plots of differentially abundant eukaryotic ASVs between melt method and ice horizon. The log2-fold change in relative 18S amplicon sequence variant (ASV) read counts between (A) DM-40 and SM-40, (B) DM-20 and SM-20, (C) DM-5 and SM-5, (D) IM-40 and SM-40, (E) IM-20 and SM-20, (F) IM-5 and SM-5, (G) IM-40 and DM-40, (H) IM-20 and DM-20, and (I) IM-5 and DM-5. The labeling system includes melt method (SM for seawater melt, IM for isohaline melt, and DM for direct melt) and ice horizon (5, 20, and 40 for 0–10 cm, 10–30 cm, and 30–50 cm, respectively). When comparing two melt methods, a positive log2-fold change represents a higher abundance in the first method listed and a negative log2-fold change indicates a higher abundance in the second method. Significantly differential abundances (p < 0.1) are displayed in color by taxonomic class. (Figure S3.png)
Table S1. Final dilutions for isohaline melt treatments. (TableS1.doc)
Table S2. Phylogenetic placement of significantly differentially abundant bacterial amplicon sequence variants (ASVs) between melt method comparisons. (TableS2.xls)
The first four columns represent taxonomic identification of each significant ASV and the next six columns display the results from the DeSeq2 analysis (Love et al., 2014), including the log2-fold change and significant (p < 0.1) adjusted p-values between the seawater melt (SM), direct melt (DM), and isohaline melt (IM) comparisons. When comparing two melt methods, a positive log2-fold change represents a higher abundance in the first method listed and a negative log2-fold change indicates a higher abundance in the second method.
Table S3. Phylogenetic placement of significantly differentially abundant eukaryotic amplicon sequence variants (ASVs) between melt method comparisons. (TableS3.xls)
The first four columns represent taxonomic identification of each significant ASV and the next six columns display the results from the DeSeq2 analysis (Love et al., 2014), including the log2-fold change and significant (p < 0.1) adjusted p-values between the seawater melt (SM), direct melt (DM), and isohaline melt (IM) comparisons. When comparing two melt methods, a positive log2-fold change represents a higher abundance in the first method listed and a negative log2-fold change indicates a higher abundance in the second method.
Acknowledgments
The authors would like to acknowledge Avishek Dutta and Jesse Wilson (University of California San Diego) for their mentorship, along with Natalia Erazo, Elizabeth Connors, and Benjamin Klempay (University of California San Diego) for their input on data analysis.
Funding
This work was funded by NSF OPP 1821911 (USA), VR 2018-04685 (Sweden), and Formas 2018-00509 (Sweden). E Chamberlain was funded by the NSF Graduate Research Fellowship and the University of California San Diego Cota Robles Fellowship. JP Balmonte was funded by a Carl Tryggers Foundation Fellowship. A Torstensson was funded by the Swedish Research Council (2017-06205).
Competing interests
To the best of our knowledge, no authors have a conflict of interest that would impact this research. J Bowman is an associate editor at Elementa. He was not involved in the review process of this article.
Author contributions
E Chamberlain conducted the laboratory analyses, statistical analyses, developed the first draft of the manuscript, and designed all figures. J Bowman conducted the bioinformatics analyses and provided mentorship. J Bowman, A Fong, and P Snoeijs Leijonmalm developed the initial experimental design, while J Bowman, E Chamberlain, JP Balmonte, and A Torstensson collected samples and executed experimental setup at the field site. All authors contributed to revising the final manuscript and approve the submitted version for publication.
References
How to cite this article: Chamberlain, EJ, Balmonte, JP, Torstensson, A, Fong, AA, Snoeijs-Leijonmalm, P, Bowman, JS. 2022. Impacts of sea ice melting procedures on measurements of microbial community structure. Elementa: Science of the Anthropocene 10(1). DOI: https://doi.org/10.1525/elementa.2022.00017
Domain Editor-in-Chief: Jody W. Deming, University of Washington, Seattle, WA, USA
Associate Editor: Lisa A. Miller, Institute of Ocean Sciences, Fisheries and Oceans Canada, Sidney, BC, Canada
Knowledge Domain: Ocean Science
Part of an Elementa Special Feature: Insights into Biogeochemical Exchange Processes at Sea Ice Interfaces (BEPSII-2)
University of California Press