Changes in Sea-Ice Protist Diversity With Declining Sea Ice in the Arctic Ocean From the 1980s to 2010s

The large declines in Arctic sea-ice age and extent over the last decades could have altered the diversity of sea-ice associated unicellular eukaryotes (referred to as sea-ice protists). A time series from the Russian ice-drift stations from the 1980s to the 2010s revealed changes in community composition and diversity of sea-ice protists from the Central Arctic Ocean. However, these observations have been biased by varying levels of taxonomic resolution and sampling effort, both of which were higher in the early years at drift stations on multiyear sea ice (MYI) in the Central Arctic Ocean. We here combine the Russian ice-drift station data with more recent data to (1) identify common sea-ice protists (in particular diatoms) in drifting sea ice of the Central Arctic Ocean; (2) characterize the potential change in such communities over 35 years in terms of species number and/or community structure; and (3) relate those shifts to relevant environmental factors. In terms of relative abundance, pennate diatoms were the most abundant sea-ice protists across the Arctic, contributing 60% on average of counted cells. Two pennate colony-forming diatom species, Nitzschia frigida and Fragilariopsis cylindrus, dominated at all times, but solitary diatom species were also frequently encountered, e.g., Cylindrotheca closterium and Navicula directa. Multiyear sea ice contained 39% more diatom species than first-year ice (FYI) and showed a relatively even distribution along entire sea-ice cores. The decrease in MYI over the last decades explained the previously reported decreases in sea-ice protist diversity. Our results also indicate that up to 75% of diatom species are incorporated into FYI from the surrounding sea ice and the water column within a few months after the initial formation of the ice, while the remaining 25% are incorporated during ice drift. Thus, changing freeze-up scenarios, as currently witnessed in the Central Arctic, might result in long-term changes of the biodiversity of sea-ice protists in this region.


INTRODUCTION
The age and extent of Arctic sea ice have dramatically declined over the last decades (Stroeve and Notz, 2018) with likely negative consequences for the diversity of flora and fauna that inhabit sea ice (Melnikov, 2005;Bluhm et al., 2017). The liquid-filled network of brine channels and pockets in sea ice is inhabited by a high diversity of organisms ranging from bacteria and Archaea to unicellular and multicellular eukaryotes, termed sympagic pro-and eukaryotes (Bluhm et al., 2017). Unicellular sympagic eukaryotes, called sea-ice protists here, are a phylogenetically diverse group which comprises photo-, mixo-and heterotrophic taxa Melnikov et al., 2002;Poulin et al., 2011). Mixo-and heterotrophic taxa are commonly represented by species within the dinoflagellates and ciliates while the collective term ice algae is frequently used for phototrophic protists in sea ice which are generally dominated by diatoms. Species number estimates range from 1027 to 1276 taxa across the Arctic Bluhm et al., 2017). The diversity of sea-ice protists is influenced by geographic location (Niemi et al., 2011;Hardge et al., 2017a), season and the age of sea ice. The older the ice, the more complex its structure, leading to increased diversity of the sea-ice inhabiting flora and fauna (Melnikov, 2009;Hardge et al., 2017b). Within the sea ice, different communities are recognized both on horizontal (from local patchiness to geographic differences) and vertical (along the ice column or ridge) dimensions (Syvertsen, 1991;Horner et al., 1992;Fernández-Méndez et al., 2018). Bottom and subice communities are characterized by a dominance of marine pennate diatoms and the mat-forming centric diatom Melosira arctica (Horner et al., 1992;Różańska et al., 2009;Fernández-Méndez et al., 2014;Poulin et al., 2014;Campbell et al., 2018), while surface melt pond communities may contain different freshwater taxa (Kilias et al., 2014), but usually in low biomass due to low nutrient concentrations on sea ice (Garrison et al., 2003). Brackish water melt ponds that have melted through the ice can sustain higher biomass through nutrient exchange with the underlying water column (Lee et al., 2011;Mundy et al., 2011) and are characterized by algal communities dominated by diatoms, including large algal aggregates (von Quillfeldt, 1997;Lee et al., 2011;Assmy et al., 2013;Fernández-Méndez et al., 2014. Diatoms, as well as the prymnesiophyte Phaeocystis pouchetii, can also be found at the snow-ice interface when the ice becomes flooded (McMinn and Hegseth, 2004;Fernández-Méndez et al., 2018). Landfast ice algal communities are distinct from offshore pack-ice communities , reflecting the age (first-year) and structure (generally flat) of landfast sea ice as well as shallower water depth. Algal biomass is not uniformly distributed in sea ice, with patchiness related to snow depth, distribution of brine channels and ice melt on smaller scales (e.g., Mundy et al., 2007;Campbell et al., 2018;Lange et al., 2019) and different nutrient regimes on both small and larger scales (Gradinger, 2009;Dalman et al., 2019).
Changes in sea-ice extent and structure, and enhanced melting, affect organisms living inside the ice matrix. Light conditions under the ice are modulated seasonally by day length and locality, snow depth and other properties, ice thickness as well as particle content in the ice (Leu et al., 2015;Katlein et al., 2019). During the melt season, ponds develop on top of the ice and increase light transmission from 5-15% below white ice to 40-70% below ponds Katlein et al., 2019). A continuation of the observed decline in sea-ice extent and thickness will increase the amount of light penetrating into the Arctic Ocean (Nicolaus et al., 2012), which will further enhance melting and alter the upper ocean ecosystem (Flores et al., 2019). In particular the snow thickness on top of the ice controls light penetration and, thus, the accumulation of ice algal biomass, with highest biomass under thin snow cover (Leu et al., 2015). However, if snow and ice cover are very thin, the ice algae may receive damaging levels of irradiation during spring , leading to under-ice blooms of phytoplankton (Arrigo et al., 2012;Assmy et al., 2017). Seasonal warming and desalination of sea ice during the melt season cause sloughing off of algae. Thin snow cover on the sea ice can cause early meltout of the ice algal bloom (Leu et al., 2015). Thus, maximum biomass may be observed under intermediate thickness of snow cover (Mundy et al., 2005). In the Central Arctic Ocean, much of the primary production is often generated by sea-ice algae rather than phytoplankton (Gosselin et al., 1997;Fernández-Méndez et al., 2015), while they contribute an important but relatively small fraction in landfast ice and seasonal ice on shelf seas (Gradinger, 2009). Nevertheless, due to their early bloom, they form a significant food source for grazers (Søreide et al., 2010).
The overall response of ice algae to climatological forcing is complex due to the anticipated changes in ice characteristics but also to the increase in Arctic precipitation (Bintanja and Selten, 2014), the timing of precipitation relative to open water and ice-covered seasons (Merkouriadi et al., 2017) and whether the precipitation falls as snow or rain (Bintanja and Andry, 2017). In the Central Arctic, the sea ice has become thinner (Kwok, 2018;Stroeve and Notz, 2018) and the freshwater content and stratification of the upper water column have increased at least in the Amerasian Basin due to higher volumes of riverine runoff along the Arctic coast (Prowse et al., 2015;Polyakov et al., 2018). Increased snow load may lead to negative freeboard, giving rise to infiltration communities (McMinn and Hegseth, 2004;Fernández-Méndez et al., 2018), as frequently recorded from Antarctic sea ice (e.g., Robinson et al., 1997;Kristiansen et al., 1998). In late summer and autumn, the inherent thinner sea ice leads to melt ponds, which subsequently melt through making a connection with the underlying water which results in the development of new habitat and growth for sea-ice algae (Lee et al., 2011). The net effect of changes in the sea-ice regime on ice algal primary production for the Arctic remains uncertain, with large regional and latitudinal differences in seasonal seaice extent and thickness (Barber et al., 2015;Leu et al., 2015;Tedesco et al., 2019). Some modeling studies indicate a decrease in ice algal areal production on a pan-Arctic scale (Dupont, 2012;Duarte et al., 2017) while others forecast increasing values (Matrai and Apollonio, 2013).
Summer sea-ice extent in the Arctic Ocean has declined by over 30% since the satellite record began in 1979 (Vaughan et al., 2013;Meier et al., 2014;Perovich et al., 2018), which is over a similar period as this synthesis study. With declining multiyear ice (MYI) extent, the first-year ice (FYI) coverage has increased in relative proportion and importance, and large parts of the Arctic are currently devoid of sea ice for extended periods of time (Arrigo et al., 2011). Effects of reduction in sea ice on ice algal diversity are uncertain, but consequences for iceassociated taxa seem inevitable. This may especially be true in the Central Arctic where the transition from dominance of longlived to short-lived sea ice has been most prominent (Stroeve and Notz, 2018). A decline in the number of sympagic eukaryote species has been suggested to occur between the 1980s and 2000s (Melnikov, 2009), and a change in community structure over four decades was also suspected (Bluhm et al., 2017). Since the dataset used in Bluhm et al. (2017) included few recent ice cores from the Central Arctic Ocean, we here expand their study by including data from 2000-2015. We build a metaanalysis on observations of generally broad distribution ranges of ice algal taxa in comparable habitats (Okolodkov, 1992;Poulin et al., 2011). Here, we aim to (1) identify common sea-ice algae and other single-celled eukaryotes (collectively termed sea-ice protists) of the Central Arctic Ocean; (2) identify and characterize the change in those communities over nearly four decades in terms of species number and/or community structure, and (3) relate those shifts to environmental variables.

Ice Cores
The data considered in this meta-analysis are based on 224 icecore samples originating from 14 field studies covering the time period from 1980 to 2015 (Table 1). There were 165 ice cores including the entire ice column (=whole cores) in the dataset. Each individual ice core including at least the ice-water interface (=bottom) section was considered as one sample (= replicate) in data analyses. The dataset consists of 101 MYI, and 123 FYI samples (Figure 1). While the focus of this study was on the Central Arctic basins, samples from north of Svalbard (near Yermak Plateau) were included to balance the decrease in sampling effort and level of taxonomic identification from the 1980s to the 2015.
Cores from different years were pooled to decadal scale for statistical analyses (with year 2000 being considered as part of the 1990s). Further, cores were classified based on ice type (MYI and FYI), month, geographic sampling region (Amerasian Basin, Siberian Shelf-Slope, Transpolar Drift, and North of Svalbard), solar elevation angle [a complementary angle of solar zenith angle; see Bluhm et al. (2018)] and field-measured ice thickness.

Datasets
Three separate datasets from different research institutes were included in the study ( Table 1). These datasets differed in methodology, sample preservation and taxonomic identification protocols. Main data sources originated from the Arctic and Antarctic Research Institute (AARI), the Alfred Wegener Institute Helmholtz Center for Polar and Marine Research (AWI) and the Norwegian Polar Institute (NPI). The AARI data originated from Russian transpolar drift expeditions as well as the icebreakers CCGS Des Groseilliers and Akademik Fedorov from 1980 to 2011 (Table 1 and Figure 1). Sampling procedures have been described in detail in Melnikov et al. (2002) and Melnikov (2005). Sampling effort and level of taxonomic identification vary within the dataset, with generally more detailed sampling during the early expeditions (Table 1). Multi-year and FYI ice cores were collected with a 12-cm AARI-type ice auger and/or a 10cm fiberglass-barrel CRREL-type corer. Cores were first sectioned into 10 or 20-cm segments, then transferred into one or twoliter plastic containers and finally melted in such containers at room temperature (20-22 • C) for 4-5 h, with no addition of filtered seawater (Rintala et al., 2014). Subsamples of 100-200 mL from each ice section were preserved in formaldehyde buffered Unit indicates original unit in the dataset (per = abundance percentage, rel = relative abundance from 1 to 5, and L −1 = cells per liter). MYI indicates the percentage of samples from multiyear ice, the rest being samples from first-year ice. Standard deviation in brackets is given for ice thickness (Ice) and water depth (Depth). Total refers to the total number of ice cores (=samples), and Whole the number of whole cores covering the entire ice column. Allocation of ice cores to ice type (first-year ice, multiyear ice) by decade. Numbers in brackets indicate the total number of ice cores for each ice type. The map was created using the PlotSvalbard package (Vihtakari, 2018) with GEBCO one-minute grid and Natural Earth data for bathymetry (https://www.gebco.net/data_and_products/ gridded_bathymetry_data/gebco_one_minute_grid/) and land shapes (http://www.naturalearthdata.com/downloads/50m-physical-vectors/).
with sodium acetate (final formaldehyde concentration of ca. 1%). Samples for cell enumeration and species identification were settled in Zeiss-type settling chambers for at least 12 h before counting with a Zeiss inverted light microscope (Utermöhl, 1931). Russian taxonomic experts identified ice protist cells containing pigment from each ice core sample to the lowest possible taxonomic level. Horizontal transects across the bottom of the chamber were counted at 450× magnification for small, abundant organisms. The number of transects was dependent on the relative number of cells present in the chamber, but usually 1/10 of the chamber bottom was counted. A single transect through the center of the chamber was counted at 300× magnification for large, rare protists. The AWI data were based on samples for species enumeration taken during the RV Polarstern PS78 TransArc expedition from 15 August to 23 September 2011. Three replicate ice cores within 1 m 2 of sea ice were drilled with a 9-cm diameter corer (Kovacs Enterprises, Roseburg, United States) at each of stations PS78_218, PS78_222, and PS78_227 (Schauer, 2012). The bottom 5 cm was cut and melted in filtered sea water at 4 • C. Sub-samples for microscopy were preserved in hexamethylenetetraminebuffered formaldehyde (final concentration 0.5%) and stored in brown glass bottles. For microscopic analyses, an aliquot of 20-50 mL was transferred to a settling chamber where the cells were allowed to settle for 48 h. Based on studies of Edler (1979), at least 400 cells of the dominant species or groups were counted with an inverted microscope. Ice protist cells were identified into groups and size classes of low taxonomic resolution (Supplementary Table S1).
The NPI ice core data were obtained during two sampling campaigns. First, the Norwegian Polar Institute's Centre for Ice, Climate and Ecosystems (ICE) cruise on RV Lance in the marginal ice zone (MIZ) north of Svalbard from 28 April to 14 May 2011 [see Nomura et al. (2013) for details]. Second, the Norwegian young sea ICE (N-ICE) campaign between January and June 2015 when RV Lance was frozen into sea ice at about 83 • N of Svalbard and allowed to drift with the pack ice to the ice edge at about 80 • N (Granskog et al., 2016;Olsen et al., 2017a). During N-ICE, the age of the sampled ice (FYI and MYI) was determined from salinity profiles and ice thickness (Olsen et al., 2017b) and corroborated by the oxygen isotopic composition of the ice (Granskog et al., 2017). Only FYI cores were collected during the ICE cruise, and both ice types were collected during N-ICE where MYI was second-year ice. Ice cores were retrieved with 9 or 14-cm diameter corers (Mark II coring system, Kovacs Enterprises). During the ICE cruise, only the bottom 3 cm of the ice cores was analyzed for sea-ice protist taxonomy, while during N-ICE most cores (97% , Table 1) contained the entire ice column often with a 10-cm bottom section followed by 10-95 cm long segments up to the ice-atmosphere interface (=top). The sections were melted overnight at room temperature in opaque plastic containers with lids without addition of filtered seawater (Rintala et al., 2014). Samples for ice protist taxonomy were collected in 100 mL brown glass bottles and fixed with glutaraldehyde and hexamethylenetetramine-buffered formaldehyde at final concentrations of 0.1% and 1%, respectively. Samples were stored cool (5 • C) and dark until analysis at the Institute of Oceanology Polish Academy of Sciences (IO PAN) by Magdalena Różańska-Pluta and Agnieszka Tatarek (N-ICE), and at the Norwegian Polar Institute by Philipp Assmy (ICE), who was trained at IO PAN in Arctic sea-ice protist identification prior to analysis of the ICE samples. Depending on the cell density of the sample, a volume of 10, 25, or 50 mL was settled in sedimentation chambers (Hydro Bios, Kiel, Germany) for 48 h. Cells were identified and enumerated using a Nikon inverted light and epifluorescence microscope (Nikon TE300, Ti-U and Ti-S, Tokyo, Japan) following Throndsen (1995) using magnifications 100-600× depending on the size of the organisms examined. A minimum of 50 cells of the most abundant species were counted, resulting in 95% confidence intervals being ±28% from the mean estimate (i.e., precision) assuming a normal distribution (Edler, 1979;Edler and Elbrächter, 2010). In order to aid species identification under light microscopy, selected samples from the ICE cruise were prepared for scanning electron microscopy (SEM). One set of raw samples was washed six times with distilled water to remove salt while another set of samples was additionally treated with potassium-permanganate and hydrochloric acid to remove organic matter. The latter treatment was used to clean the diatom frustules for better identification. The samples were then dried on round cover slips (10 mm in diameter), mounted on 12.5 mm diameter aluminum stubs (Plano GmbH, Wetzlar, Germany), sputter-coated with gold-palladium and observed with a Quanta FEG 200 SEM (FEI, Hillsboro, United States).

Taxonomic Nomenclature
Sea-ice protists were identified to the lowest possible taxonomic rank, which ranged widely from phylum, class, order, family, genus, and species to variety or forma levels. Further, the taxonomic nomenclature has changed considerably within the four decades of data coverage in this study causing difficulties in comparing ice protist communities over time. Taxon names reported in the original datasets were corrected, updated and unified using a three-step protocol. First, all reported taxon names were passed through the World Register of Marine Species database (WoRMS Editorial Board, 2018) to confirm the validity of a name using the taxize package (Chamberlain and Szöcs, 2013) for R (R Core Team, 2018). Second, the returned taxon names were validated by protist taxonomy experts within the author team, and when needed, checked against the AlgaeBase database 1 , which represents more up-to-date classification of protist taxa than WoRMS, but is not programmatically available due to copyright restrictions. Finally, the taxon names were manually edited to make the taxonomic ranks across datasets as comparable as possible. Since the taxonomic ranks are not consistent within the compiled dataset, we use the term "taxa" when referring to all sea-ice protists and the term "species" when referring to species level or lower. Thus, our data analyses also include variety and forma as separate species entries.

Data Analysis
Most studies reported species abundances as absolute abundances (cells L −1 , expressed as L −1 from here on) except for the early 1980-1981 AARI data that provided percentage abundances (PA), and PAICEX-2007 data expressed as relative abundances (RA) on a scale from 1 to 5 (Table 1). Ice cores with less than five taxa (total of six cores) were removed from analyses because these cores were clearly outliers in the dataset representing conditions that were not comparable to other ice cores.
The original cell counts for ice core sections were converted to integrated cells per square meter of sea-ice estimates (Ab, unit expressed as m −2 from here on) for each taxon in an entire ice core as follows: Where y ij is the cell count in m −3 (1000 × cells L −1 ) for ice-core section i and taxon j from a species abundance matrix Y = y ij , of size (n × p) with sections of an ice core as rows (i = 1......n) and taxa as columns (j = 1.....p); h i the height of the ice-core section i in meters; and n the number of ice-core sections within a core.
Percentage abundances (PA) or relative abundances (RA) were first summed up by ice core and then converted to average abundance percentages (AP) by dividing with the overall sum of values within an ice core (see below): Where y j+ is row sum (i.e., the sum over all ice-core sections within a core) for taxon j; and Y ++ is overall sum of the entire percentage or relative abundance matrix for an ice core. The AWI and ICE-cruise datasets had 1-2 orders of magnitude higher total cell count estimates (with medians of 8.09 and 5.75 × 10 9 cells m −2 , respectively) compared to AARI (median 4.22 × 10 7 cells m −2 ) and N-ICE (median 4.49 × 10 8 cells m −2 ). These differences may be explained by seasonal and spatial variability in sampling, and possibly by the differences in sample preservation, counting and abundance calculation methods. Further, the cell count estimates were strongly heteroscedastic and non-normally distributed when grouped using dataset identifier, ice type, year or decade. Despite the differing abundance values and metrics, the relative contributions of taxa within samples were likely comparable and therefore used in this study as explained below.

Abundance Metrics
Two abundance metrics were used: (1) Average abundance percentage (AP): the mean percentage contribution of a species to the total abundance of a sample, calculated as an arithmetic mean of percentages (Martin et al., 1946;Bluhm et al., 2018): Where y ij is the integrated, percentage or relative abundance value for taxon j in ice core i, and Y i+ the sum of all taxa (row sum) in ice core, and (2) Frequency of occurrence (FO): the proportion of samples containing one or more cells of a given taxon. The effect of a dataset which systematically did not identify a taxon was removed from the overall AP and FO estimates. The AP estimates were converted to proportions before data analyses. We use set terminology (qualitatively) connected to the metrics. Abundance percentage (AP): abundance, abundant. Frequency of occurrence (FO): frequency, occurrence, common, rare, encountered. Number of taxa (or species): diversity, diverse, species-rich, uniform.

Statistical Methods
Higher than genus level taxonomic ranks were removed from the dataset for diversity and community analyses. A genus was excluded from species counts in diversity analyses if there were other species of the genus present in an ice core. Varieties were treated as separate taxa.
Species richness and vertical distribution of sea-ice protists were examined using generalized linear mixed models (GLMMs). The discrepancies in species identification, sampling effort and location made it difficult to compare species counts over decades or ice types. Consequently, the diversity analyses were restricted to diatoms only since these taxa are easier to identify, less affected by different sample fixation approaches and, thus, presumably less biased than flagellates, dinoflagellates and ciliates.
The explanatory variables selected for the comparison (region, dataset and season to correct for unbalances in the compiled dataset as well as ice type, decade and ice thickness as actual predictor variables) were all correlated. Since ice type (FYI versus MYI) was the most important factor describing diatom diversity, and ice type often formed interactions together with other explanatory factors, the analyses were run separately for FYI and MYI using each ice core as random intercepts to remove the bundled correlations. This procedure simplified the model fitting and removed model convergence and over-dispersion problems encountered using other variables to correct for biases in the dataset. The results echoed the understanding of the dataset established during the data exploration phase and did give similar Central estimates than more complicated models. The exception to the simplified GLMM fitting were the species richness models for vertical distribution, which were corrected by the proportion of protists in an ice-core section by using the proportions as random intercepts. Poisson family log-link function was used to linearize the GLMMs for species richness (count data), while binomial distribution was used for proportion data (Bolker et al., 2009;O'Hara and Kotze, 2010). The models were fitted using the lme4 (Bates et al., 2014) package for R statistical programming environment (R Core Team, 2018) and the Laplace Approximation routine (Bolker et al., 2009). Model estimates and confidence intervals were back-transformed to counts using the effects package (Fox, 2003). Changes in diatom diversity with ice thickness were curve fitted by local polynomial regression (LOESS fit). Multiple comparisons among variable levels were conducted using Tukey tests and Holm-adjusted p-values using the multcomp package (Hothorn et al., 2008). Due to the complex biased dataset, the GLMM results should be interpreted with caution especially if the significance level is close to the alpha limit (0.05).
The community composition was examined using principal component analyses (PCA) with square-root transformed abundance proportion matrix [ AP j / 100% ; equal to Hellinger transformation in Legendre and Gallagher (2001)] using the rda function from the vegan package (Oksanen et al., 2017). Higher than genus level taxonomic ranks were removed from the dataset prior to analysis, leading to genus, species and variety/forma being considered "taxa" (= columns) in the PCAs. Sampling region, original dataset and examined ice-core length explained much of the inertia in the unconstrained community PCA. Since these factors were not the explanatory variables of interest, they were removed by conditioning the PCA orientation (also called partial PCA). Explanatory variables were fitted to the PCA ordinations using the envfit function from the vegan package, and R 2 values together with graphical presentation were used to examine the explanatory power of each variable.

Most Common Taxa
The total number of sea-ice protist species encountered in the combined dataset was 201 or 221 if varieties/forma were counted separately (Supplementary Table S1). These taxa originated from 120 genera. Pennate diatoms was the group with most species (with species and varieties/forma included) followed by centric diatoms, flagellates, dinoflagellates and ciliates ( Table 2). The low species richness for flagellates and ciliates (13% of total number) likely reflects the difficulty of identifying these groups to species level using light microscopy and their fragility when it comes to melting procedure and sample preservation.
Pennate diatoms were present in almost all samples and were by far the most abundant group (measured as average abundance percentage, AP) in the entire dataset followed by flagellates, dinoflagellates, centric diatoms, chlorophytes, xanthophytes, and chrysophytes (Table 3 and Supplementary Table S2 for overview for all taxonomic groups). Flagellates and ciliates were frequent, although not very abundant, in all datasets except for the AARI dataset, from which flagellates were less abundant and ciliates missing. This likely reflected the method of direct melting of ice-core sections before fixation (I.A. Melnikov, pers. obs.). Pennate diatoms had higher abundance in MYI compared to FYI, and a similar pattern was present for chlorophytes (Figure 2). Only unique species, varieties and forma are included in cell counts.
Frontiers in Marine Science | www.frontiersin.org Dinoflagellates were generally more abundant in FYI than in MYI, although the diversity was rather similar. Dominant species for each of the major taxonomic groups were: pennate diatoms: Nitzschia frigida and Fragilariopsis cylindrus; centric diatoms: Attheya septentrionalis; dinoflagellates: Polarella glacialis; flagellates: Groenlandiella brevispina; and ciliates: Mesodinium rubrum ( Table 3). The relative contribution of Fragilariopsis cylindrus, Pseudo-nitzschia delicatissima complex and Polarella glacialis increased towards a FYI regime, whereas melt pond specialists, chlorophytes Chlamydomonas nivalis (also known as snow algae; see Procházková et al. (2019) for its taxonomic status) and Trochiscia cryophila, were abundant in MYI with no record from FYI in this dataset. The genus Nitzschia Hassall was common during the MYI years because Nitzschia was not identified to species level in the early AARI dataset and AWI samples. The centric diatom Melosira arctica had much higher frequency of occurrence in MYI, but because of its colonial form (i.e., colonies rather than individual cells were counted), the abundance counts were rather low. The same reason for low abundance also applied to Attheya spp., which are ephiphytic diatoms with patchy distribution. Attheya septentrionalis, which is ephiphytic on M. arctica, can be among the most common members of a sea-ice community, despite low occurrence of their most common supporting algae (von Quillfeldt et al., 2003).

Diversity
Diatom diversity was highest in MYI samples from North of Svalbard, followed by the Amerasian Basin samples ( Figure 3A). There were no substantial differences in diversity among regions in FYI. The N-ICE samples had higher diatom diversity than the AARI samples in MYI, while the differences in FYI were unclear due to variable occurrences in samples, but AARI data appeared to have lower species richness than the other datasets ( Figure 3B). Ice thickness (x) was generally correlated with ice type (probability of ice being MYI p (x) = e −15.17+11.53x 1+e −15.17+11.53x , p = 0.002 for the slope, where x is ice thickness in meters), but did not explain variability in diatom species counts alone ( Figure 3C). Ice type, on the other hand, was the best explanatory variable of diatom species number: diatom diversity was 39% higher in MYI than in FYI, based on respective Central estimates from GLMM of 16.3 and 11.7 diatom species (Figure 3D). The GLMM analysis (Supplementary Tables S3, S4) did not indicate long-term changes in diatom diversity for FYI, but diversity decreased in MYI from the 1980s to the 2000s (Figure 3D). From 2000 to 2015, the diversity in MYI increased again based on samples collected North of Svalbard. The decline in diatom diversity over decades in MYI was also detectable in the Transpolar Drift region, although low sample sizes and multiple confounding factors made the relationship nonsignificant (Supplementary Figure S1).
The diatom diversity was influenced by seasonality in FYI (Figure 4A), and this effect was best explained (and in contrast to the overall dataset not biased by sampling effort or method) by a subset of the N-ICE data that was collected from newly formed sea ice north of Svalbard from May to June 2015 ( Figure 4B). The dataset demonstrated an increase in diatom diversity from the initial three species on average in beginning of May to approximately 20 species on average in the beginning of June within the same year and ice floe.

Community Composition
The sea-ice protist community, while partially overlapping, was clearly separated by dataset and sampling region as indicated by the principal component analysis (Table 4 and Figure 5). The percentage variability explained (32.9%) by the first two axes in PCA reflects that there were more species than stations in the analysis (thus, >30% inertia explained is considered high for this type of data). Fragilariopsis cylindrus, Navicula pelagica, Navicula transitans, Pseudo-nitzschia delicatissima complex, and Polarella glacialis were more abundant in N-ICE and ICE samples from North of Svalbard than in AARI and AWI samples from the Central Arctic. This caused partial separation in datasets ( Figure 5B) and of regions ( Figure 5C). Fragilariopsis oceanica and Nitzschia polaris were more abundant in the AARI dataset than in other datasets, while some pennate diatom genera such as Nitzschia and Navicula Bory were generally more abundant in the Central Arctic than North of Svalbard due to previously mentioned differences in taxonomic resolution. Decade, ice type , through datasets (B), ice thickness and type, with lines as LOESS fit (C), and through the four decades (D) included in the meta-analysis ("All" implies decades combined for first-year ice and multiyear ice). Light and dark blue dots indicate diatom diversity in first-year ice (FYI) and multiyear ice (MYI), respectively. Red dots indicate GLMM model estimates, with their 95% confidence intervals. Letters on top of x-axis represent multiple Tukey comparisons using the GLMMs, Holm correction and p-value limit of ≤0.05. and month were largely overlapping in the principal component analysis (Figures 5D-F). Since these factors were correlated with regional and dataset effects, their actual effect on the community structure was unclear (Table 4). Essentially, the samples from the 2010s overlapped with all other decades (Figure 5D), FYI encompassed that of MYI (Figure 5E), and the months of April/May encompassed all other months (Figure 5F). Conditioning the PCA reduced the explanatory power of the analysis, but also indicated that time (year and decade, correlated with ice thickness and type) and ice thickness (correlated with ice type) may have influenced the community composition (Table 4 and Supplementary Figure S2).

Vertical Distribution of Sea-Ice Protists in Cores
The ice-water interface (=bottom) contained the highest average proportion of sea-ice protist cells in FYI while the average

Taxonomic Inventories and Sampling Biases
The species richness of microalgae and other protists in sea ice is generally high, but the inventory of these singlecelled eukaryotes is inconsistent among studies, which makes assessment of temporal changes challenging. The first pan-Arctic inventories based on morphological identifications reported 1027 single-celled eukaryotes inhabiting Arctic sea ice Daniëls et al., 2013). The subsequent synthesis by the Circumpolar Biodiversity Monitoring Program (CBMP) Sea Ice Biota Expert Network documented that increased effort still increases the inventory, which currently includes more than 200 additional taxa for a total of 1276 sympagic microalgae and other protists (Bluhm et al., 2017). This is considerably more than the 221 species, varieties and forma counted in this study, but many of the taxa in the inventory above were not eukaryotic algae (e.g., phototrophic bacteria), or not identified to species level, or with variable taxonomic resolution. Standard microcopy counts performed in this study did not take into account the diversity of smaller mixo-and heterotrophic microbial eukaryotes from sea ice. Molecular techniques, e.g., 18S metabarcoding, indicate that the actual sea-ice protist community is even more diverse than the morphological taxonomic inventory suggests, including many protists that are difficult to identify with microscopy (Comeau et al., 2013;Kilias et al., 2014;Hardge et al., 2017a,b). Thus, the increased number of protist species occurring in sea ice in recent compilations is largely driven by changing scientific methods of species identification and not by changing ice regimes.
Differences between the datasets largely affect a true evaluation of long-term development of unicellular eukaryote diversity in Arctic sea ice. This is partly related to the different sampling and analyses approaches in the various studies, but also because different regions of the Arctic Ocean were sampled over time. This involves direct ice melt used for the NPI and AARI datasets and ice melt with the addition of filtered sea water as salinity buffer in the AWI dataset. Such differences in melting approaches cannot only affect activity estimates of sea-ice biota (e.g., Campbell et al., 2019) but also estimates of abundance and composition (e.g., Garrison and Buck, 1986;Rintala et al., 2014). While diatoms and other taxa with hard casings (e.g., silicoflagellates) are not impacted by osmotic stress in their morphology, delicate forms like naked flagellates might change their form (making them unrecognizable) or even dissolve or explode. Thus, with regard to our 35-year comparison, diatom occurrences provide the most robust data, while interpretation of changes in flagellate diversity could be biased due to the applied methodologies. In summary, the original dataset (and associated laboratory analysis) had the largest effect on the community composition and number of sea-ice protist taxa in our analysis, which calls for a more standardized method between several laboratories with regard to taxonomic analyses of microalgae and other eukaryotes in sea ice (e.g., Bluhm et al., 2017).

Diversity Patterns and Potential Environmental Drivers
Changes in diversity during the sampling period were linked to reductions in sea-ice thickness and concentration in the Arctic Ocean during the four decades. Multiyear sea ice has the highest diversity (species number) of sea-ice protists in the Arctic (Hardge et al., 2017a,b;van Leeuwe et al., 2018;this study). Thus, the dramatic decline of MYI over the last decades (Stroeve and Notz, 2018) has most likely led to decline in seaice flora and fauna diversity (Melnikov et al., 2002;Gradinger et al., 2010;Hardge et al., 2017a;Olsen et al., 2017b). Because of   declining MYI cover, 80% of the cores in the 2010s were from FYI, which typically had lower protist diversity. This observation is concurrent with the generally lower diatom diversity in FYIdominated Antarctic sea ice compared to Arctic sea ice (Lizotte, 2003). Based on our GLMM analysis of diatom species, we conclude that the diversity of sea-ice protists has decreased over 35 years through the relationship to MYI. This was also supported by the decline in the number of diatom species in MYI from the Transpolar Drift from 1980s to the 2010s and indicated by Melnikov (2018), who determined that the diversity of centric diatoms in the North Pole Region had decreased from 12 species in 2007 to three species in 2015, while pennate diatoms fluctuated during the same period, but showing no clear trend.
Regional differences were also apparent in the diversity patterns, specifically for MYI but not for FYI. Pennate diatoms were the most diverse group as is commonly observed in Arctic sea ice Leu et al., 2015;van Leeuwe et al., 2018). With the effect of dataset removed, the sampling region greatly influenced the sea-ice protist diversity in MYI with most diatom species in the area North of Svalbard followed by the Amerasian Basin as compared to the lower diatom species numbers in the Transpolar Drift and Siberian Shelf-Slope. This pattern likely relates to the overall current regime in the Arctic seas (Bluhm et al., 2017) and suggests a contribution by advection of Atlantic-origin species to the area around Svalbard and Pacificorigin species into the Amerasian Basin, which results in higher biodiversity in these regions. This advection effect could also explain the apparent resurgence in diatom species numbers in MYI from 2000 to 2015 after the pronounced decline in diatom diversity from the 1980s to the 2000s. The 2010s data all stem from the Atlantic-influenced region North of Svalbard, which is known to harbor characteristic protist communities (Metfies et al., 2016). Abelmann (1992) encountered the highest diatom concentrations in MYI in the Transpolar Drift between 83 and 86 • N. A mechanistic explanation for this was the incorporation of protists over the Siberian Shelf and further accumulation during freezing and melting processes as the ice floes drift across the Arctic Ocean (Abelmann, 1992;Assmy et al., 2017;Hardge et al., 2017a). Backtracking used to determine origin and approximate age of sea ice as it drifts across the Arctic Ocean towards Fram Strait (Hop and Pavlova, 2008) showed that the recent Arctic warming interrupts the transport of ice rafted matter within the Transpolar Drift (Krumpen et al., 2019), which could explain the observed decline of diatom species numbers in MYI from the Transpolar Drift (Supplementary Figure S1).

Vertical Distribution and Seasonality of Algal Communities
Both algal community structure and biomass vary vertically within the ice sheet. Generally, bottom maxima of ice algae are often observed across the Arctic in terms of abundance, biomass, and activity (Duarte et al., 2015;Leu et al., 2015;van Leeuwe et al., 2018), with some exceptions (von Quillfeldt et al., 2003). The bottom 10 cm contains most of the ice algae, as indicated by our median values of 60%. In our estimates for vertical ice algae distribution, we calculated the percentage contribution of cells for each ice-core section and used these relative abundances to examine the vertical distribution, unlike in other studies that used absolute cell abundances (Gradinger, 1999(Gradinger, , 2009). While our method might have given too much weight to some cores that contained few ice algae, the method is not biased by a few very high abundances during the peak ice algal bloom.
Seasonality affects ice protist diversity as well as their blooms and production (Barber et al., 2015). During the ice algal bloom, the bottom ice communities are predominantly represented by colonial pennate diatoms, e.g., Nitzschia frigida and Fragilariopsis cylindrus, while some solitary cells are also frequently encountered, e.g., Cylindrotheca closterium and Navicula directa. This effect of seasonality was strongest on diatom diversity in FYI while its effect on the number of diatom species was rather low in MYI. This difference can be attributed to the fact that MYI already starts with a seeding stock of iceassociated species incorporated during previous growth seasons  (Olsen et al., 2017b), which was also reflected in the more evenly distributed diatom abundance along the ice column in MYI ( Figure 6B). Pennate diatoms with a benthic lifestyle are particularly well adapted to the sea ice-environment and MYI provides a more stable and persistent habitat than the more ephemeral FYI, reflected in the high number of benthic diatoms unique to MYI (Supplementary Table S6). Benthic diatoms (diversity described by e.g., Karsten et al., 2012Karsten et al., , 2019 can be incorporated into sea ice during ice formation and deep convection on the shallow shelves (Abelmann, 1992) as has been shown for shallow regions such as the Chukchi Sea (von Quillfeldt et al., 2003) and the Laptev Sea (Tuschling et al., 2000).
The much higher diversity in benthic than pelagic diatom species could partly explain the generally higher diatom abundance in MYI compared to FYI. In contrast, during the colonization of FYI from the water column as the sea ice forms (Gradinger and Ikävalko, 1998;Różańska et al., 2008) or by exchange between the water column, melt ponds and the sea ice during brine rejection (Hardge et al., 2017b), many pelagic protist taxa are incorporated into the newly formed sea ice. However, the majority of species do not thrive in sea ice, with the exception of the cryopelagic species mentioned below, which is supported by the high proportion of typical pelagic taxa unique to FYI while they are notably absent among the taxa unique to MYI   Table S6). This fits with the observation that the protist composition in newly formed sea ice resembles that of the underlying water column it was formed from, but as the ice becomes older, successional patterns tend toward dominance of typical ice-associated pennate diatoms . Indeed, colonization of FYI by typical ice-associated taxa within the first month of its formation can account for approximately 75% of the diversity found in MYI as reflected in the steep initial increase in diatom diversity during May/June in the N-ICE dataset from an Atlantic-influenced environment ( Figure 4B). This suggests that successional patterns in sea ice tend towards a typical ice-associated community within a few weeks of its formation. As many of the typical ice-associated diatom taxa are usually not very abundant in the water column, colonization from adjacent MYI floes (Olsen et al., 2017b) or resuspension of benthic or sedimented sea-ice diatoms during deep winter mixing over the shallow shelves (Abelmann, 1992;Tuschling et al., 2000) are likely important seeding sources for FYI.

(Supplementary
Seasonal accumulation of biomass can also be found in the intermediate and surface sea-ice sections (Duarte et al., 2015;van Leeuwe et al., 2018), including recent observations of snowinfiltration communities at the snow-ice interface (Fernández-Méndez et al., 2018). The ice cores containing high abundances above the bottom 10-cm section were collected during the spring, summer and autumn (April to September). Using epifluorescence microscopy, Gradinger (1999) recorded vertical differences in species composition with diatoms dominating in the bottom layers and the more mobile and smaller flagellates in the ice interior. A data compilation carried out by van Leeuwe et al. (2018) further indicated vertical differences in sea-ice protists along the ice core. These findings are contrary to our studies where diatom diversity was more homogenous throughout the core, but we agree that sampling of only the bottom layer of ice cores may underestimate both biomass and diversity of protists (probably by 15-20%; Gradinger, 2009;van Leeuwe et al., 2018).
A Changing Sea-Ice Cover -The Future Protist Diversity?
The regime shift from a MYI-dominated towards a FYIdominated Arctic Ocean (Stroeve and Notz, 2018) will not only result in a decline in sea-ice protist diversity, but will also change the relative composition of the sea-ice protist community (Hardge et al., 2017a,b), with apparent regional differences. Based on the frequency of occurrence and abundance data (Supplementary Table S2) one can infer preferences for MYI versus FYI for some of the more prominent species, which we have grouped according to known habitat predilection (Figure 7). The melt pond specialists, green microalgae Chlamydomonas nivalis and Trochiscia cryophila, show a clear preference for MYI, which is consistent with the fact that melt ponds on MYI are generally fresher than on FYI (Kim et al., 2018). Melt ponds on FYI tend to melt all the way through as the season progresses and incorporate microalgae from the water column and the bottom ice assemblage (Lee et al., 2011;Hardge et al., 2017b).
Among the sea-ice specialist taxa, Nitzschia frigida is generally abundant in sea ice (both FYI and MYI), which has also been confirmed by 18S rDNA (Hardge et al., 2017a), while Attheya septentrionalis and particularly Melosira arctica seem to prefer MYI. The notion that N. frigida will continue to thrive under the new FYI regime in the Arctic can also be inferred from the prominence of its sibling species N. stellata, forming similar FIGURE 7 | Schematic habitat preference of sea-ice protists and associated taxa/groups. Pelagic and benthic diatom taxa are more associated with first-year sea ice, and most often incorporated during freezing. When sea ice melts the protists become part of the vertical flux. arborescent colonies in Antarctic FYI (Scott et al., 1994). The cryopelagic species Fragilariopsis cylindrus thrives both in the water column and in sea ice and seems to have a preference for FYI, which is consistent with its abundance in Antarctic sea ice (Scott et al., 1994). Furthermore, F. cylindrus has been used as a paleoproxy to reconstruct seasonal ice cover in the Southern Ocean over the geological past (Zielinski and Gersonde, 1997), but has also been assessed as an indicator for cold water primarily rather than sea ice (von Quillfeldt, 2004). Although Pseudonitzschia H. Peragallo is usually characterized as a pelagic diatom genus reflected in the generally low abundance percentage in FYI and MYI, the frequent occurrence of the P. delicatissima complex in FYI in the more recent ICE and N-ICE datasets suggests that species of this genus could adapt to a more cryopelagic lifestyle under the thinner and more ephemeral sea-ice regime.
Preferences for FYI and low or no occurrence in MYI can be attributed to species advected to the Arctic Ocean. This includes the prymnesiophytes Phaeocystis pouchetii and Emiliania huxleyi as well as the phototrophic ciliate Mesodinium rubrum and typical pelagic species with representatives in the centric diatoms, as shown here by Chaetoceros gelidus (formerly C. socialis), Conticribra weissflogii (formerly Thalassiosira weissflogii), Thalassiosira antarctica var. borealis, Porosira glacialis, and dinoflagellates. Emiliania huxleyi is a species of Atlantic and Pacific origin that is advected into the Arctic during extensive blooms in the Norwegian and Barents seas (Hegseth and Sundfjord, 2008;Oziel et al., 2020) and the Bering Sea (Iida et al., 2012). Recent genomic analysis indicates that E. huxleyi should be included in the ubiquitous coccolithophore genus Gephyrocapsa Kamptner (Bendif et al., 2019). Cells of the other mentioned species usually sediment out of the surface water column at the end of the phytoplankton spring bloom (Morata et al., 2011;Wassmann and Reigstad, 2011) and typically survive as resting spores in surface sediments (Eilertsen et al., 1995;Hegseth et al., 2019). Deep convective mixing during sea-ice formation over shallow coastal seas can result in resuspension of resting spores and/or vegetative cells and incorporation into newly formed sea ice, but these species do not seem to thrive in sea ice . During the pelagic bloom typical phytoplankton species can be incorporated into sea ice during surface flooding and the formation of infiltration communities at the snow-ice interface (McMinn and Hegseth, 2004;von Quillfeldt et al., 2009;Fernández-Méndez et al., 2018). These have been commonly reported from the FYI-dominated Antarctic, but seem currently to be restricted to the Atlantic sector of the Arctic Ocean (Fernández-Méndez et al., 2018).

CONCLUSION AND OUTLOOK
Sea-ice protists are present in different compartments of the sea-ice ecosystem, including internal, ice-water interface, snowinfiltration layer and melt pond communities (Figure 7). The division between ice-associated and pelagic algae is not distinct since pelagic algae can be incorporated into sea ice. Protists of Atlantic or Pacific origin may also be present, such as Mesodinium rubrum and Emiliania huxleyi (Hegseth and Sundfjord, 2008;Iida et al., 2012;Olsen et al., 2019). This meta-analysis confirms key findings from our recent data analysis (Bluhm et al., 2017) in several regards. First, sea ice is inhabited by a microalgal community which is taxonomically diverse and often dominated by diatoms in both abundance and species richness. Second, many microalgae have a widespread distribution across the Central Arctic. Finally, monitoring strategies can only yield unequivocal results on temporal trends if conducted with consistent methodology. Standardization of biodiversity sampling by ice cores has been suggested multiple times, for instance in the N-ICE project (Norwegian Polar Institute) and Nansen Legacy Sampling Protocol (2020).
A key finding is the higher biodiversity found in MYI in this pan-Arctic dataset. Since MYI is projected to continue to decrease and possibly disappear altogether (Stroeve et al., 2012), and icefree summers may occur in the near future (Overland and Wang, 2013), a continued decrease in sea-ice protist diversity can be anticipated. In the future, there may be less diversity with respect to types of ice communities and those that still exist will have more or less the same species as in the water column. We can also expect changes and maybe enrichment in biodiversity due to increased inflow of Pacific and Atlantic species, enhancing their role in inflow areas (Carmack and Wassmann, 2006), including the contribution of brackish-water species. When some of us started working in the study area more than 30 years ago, the sub-ice community had the highest biomass while now this community is usually hard to find, at least in the Atlantic sector of the Arctic. Currently, interstitial communities typically appear at the bottom of the core and thick algal mats are less frequent below sea ice, but can still be observed within the Arctic Basin (Boetius et al., 2013). Melt pond and some sea-ice specialists are also likely to decline, while cryopelagic and pelagic species of both Atlantic/Pacific and Arctic origin are likely to increase in relative importance under the new Arctic sea-ice regime.

AUTHOR CONTRIBUTIONS
BB, PA, RG, IP, and IM contributed the datasets. PA, MP, CQ, IP, LO, and LZ contributed the taxonomic expertise with regard to revisions and tables included. MV and HH have done data processing and figures. HH led the project. HH, MV, and PA wrote a first draft of the manuscript, with subsequent input and revisions from co-authors.

FUNDING
Funding for data analysis and synthesis was provided by the former Centre for Ice, Climate and Ecosystems at the Norwegian Polar Institute, Norwegian Environment Agency, as part of the CBMP-Marine project, and the Ministry of Foreign Affairs, Norway (project ID Arctic). This work was supported by grants of the RF #0149-2018-0009, RFBR #18-05-00099, and RFBR/RGS #17-05-41197, and the Research Council of Norway (projects #244646 and #268286) and the PACES (Polar Regions and Coasts in a Changing Earth System) programme of the Helmholtz Association. Digitization of the Russian data records was financed by CBMP and the BMF-Project TRANSDRIFT ("System Laptewsee" -Das transpolare System des Nordpolarmeeres). MP was supported by the Canadian Museum of Nature, Ottawa, Canada. BB, RG, and HH were supported by the Arctic SIZE project co-funded by UiT The Arctic University of Norway and the Tromsø Research Foundation (project number 01vm/h15). The publication charges for this article have been funded by a grant from the publication fund of UiT The Arctic University of Norway.

ACKNOWLEDGMENTS
This study builds upon the State of the Arctic Marine Biodiversity Report (SAMBR, 2017) developed by the Circumpolar Biodiversity Monitoring Programme (CBMP), under the Conservation of Arctic Flora and Fauna (CAFF) Arctic Council working Group. Samples were provided from expeditions to the Arctic Ocean by the Arctic and Antarctic Research Instituted (AARI), Alfred Wegener Institute (AWI), Centre for Ice, Climate and Ecosystems (ICE), and Norwegian young sea ICE expedition (N-ICE2015) of the Norwegian Polar Institute (NPI).