Strong Seasonality in Arctic Estuarine Microbial Food Webs

Microbial communities in the coastal Arctic Ocean experience extreme variability in organic matter and inorganic nutrients driven by seasonal shifts in sea ice extent and freshwater inputs. Lagoons border more than half of the Beaufort Sea coast and provide important habitats for migratory fish and seabirds; yet, little is known about the planktonic food webs supporting these higher trophic levels. To investigate seasonal changes in bacterial and protistan planktonic communities, amplicon sequences of 16S and 18S rRNA genes were generated from samples collected during periods of ice-cover (April), ice break-up (June), and open water (August) from shallow lagoons along the eastern Alaska Beaufort Sea coast from 2011 through 2013. Protist communities shifted from heterotrophic to photosynthetic taxa (mainly diatoms) during the winter–spring transition, and then back to a heterotroph-dominated summer community that included dinoflagellates and mixotrophic picophytoplankton such as Micromonas and Bathycoccus. Planktonic parasites belonging to Syndiniales were abundant under ice in winter at a time when allochthonous carbon inputs were low. Bacterial communities shifted from coastal marine taxa (Oceanospirillaceae, Alteromonadales) to estuarine taxa (Polaromonas, Bacteroidetes) during the winter-spring transition, and then to oligotrophic marine taxa (SAR86, SAR92) in summer. Chemolithoautotrophic taxa were abundant under ice, including iron-oxidizing Zetaproteobacteria. These results suggest that wintertime Arctic bacterial communities capitalize on the unique biogeochemical gradients that develop below ice near shore, potentially using chemoautotrophic metabolisms at a time when carbon inputs to the system are low. Co-occurrence networks constructed for each season showed that under-ice networks were dominated by relationships between parasitic protists and other microbial taxa, while spring networks were by far the largest and dominated by bacteria-bacteria co-occurrences. Summer networks were the smallest and least connected, suggesting a more detritus-based food web less reliant on interactions among microbial taxa. Eukaryotic and bacterial community compositions were significantly related to trends in concentrations of stable isotopes of particulate organic carbon and nitrogen, among other physiochemical variables such as dissolved oxygen, salinity, and temperature. This suggests the importance of sea ice cover and terrestrial carbon subsidies in contributing to seasonal trends in microbial communities in the coastal Beaufort Sea.

Microbial communities in the coastal Arctic Ocean experience extreme variability in organic matter and inorganic nutrients driven by seasonal shifts in sea ice extent and freshwater inputs. Lagoons border more than half of the Beaufort Sea coast and provide important habitats for migratory fish and seabirds; yet, little is known about the planktonic food webs supporting these higher trophic levels. To investigate seasonal changes in bacterial and protistan planktonic communities, amplicon sequences of 16S and 18S rRNA genes were generated from samples collected during periods of ice-cover (April), ice break-up (June), and open water (August) from shallow lagoons along the eastern Alaska Beaufort Sea coast from 2011 through 2013. Protist communities shifted from heterotrophic to photosynthetic taxa (mainly diatoms) during the winter-spring transition, and then back to a heterotroph-dominated summer community that included dinoflagellates and mixotrophic picophytoplankton such as Micromonas and Bathycoccus. Planktonic parasites belonging to Syndiniales were abundant under ice in winter at a time when allochthonous carbon inputs were low. Bacterial communities shifted from coastal marine taxa (Oceanospirillaceae, Alteromonadales) to estuarine taxa (Polaromonas, Bacteroidetes) during the winter-spring transition, and then to oligotrophic marine taxa (SAR86, SAR92) in summer. Chemolithoautotrophic taxa were abundant under ice, including iron-oxidizing Zetaproteobacteria. These results suggest that wintertime Arctic bacterial communities capitalize on the unique biogeochemical gradients that develop below ice near shore, potentially using chemoautotrophic metabolisms at a time when carbon inputs to the system are low. Co-occurrence networks constructed for each season showed that under-ice networks were dominated by relationships between parasitic protists and other microbial taxa, while spring networks were by far the largest and dominated by bacteriabacteria co-occurrences. Summer networks were the smallest and least connected, suggesting a more detritus-based food web less reliant on interactions among microbial taxa. Eukaryotic and bacterial community compositions were significantly related to trends in concentrations of stable isotopes of particulate organic carbon and nitrogen, among other physiochemical variables such as dissolved oxygen, salinity, and temperature. This suggests the importance of sea ice cover and terrestrial carbon subsidies in contributing to seasonal trends in microbial communities in the coastal Beaufort Sea.

INTRODUCTION
Aquatic microorganisms drive global cycling of carbon, nitrogen, and many other elements by carrying out key ecosystem functions including primary production, organic matter remineralization, and transformations of inorganic compounds (Falkowski et al., 2008;Ferrera et al., 2015;Worden et al., 2015). The efficiency with which microbes perform these functions is undoubtedly influenced by their physical and chemical environment (Gilbert et al., 2012), but also by interactions with each other within microbial communities (Azam and Malfatti, 2007;Fuhrman et al., 2015;Guidi et al., 2016). The composition and function of microbial communities varies strongly with seasonal changes in coastal ecosystems including day length, solar radiation, temperature, and salinity (Gilbert et al., 2012;Cram et al., 2015;Bunse and Pinhassi, 2017), and in polar regions these seasonal changes are particularly extreme, with additional complexities including ice cover and wide variations in river runoff (Holmes et al., 2012). Climate change is warming the Arctic approximately two times faster than lower latitudes (Serreze and Barry, 2011), and is amplifying seasonal variations in temperature (Serreze and Barry, 2011), ice extent (Stroeve et al., 2012), and river flow (McClelland et al., 2006;Morison et al., 2012). Moreover, increased river runoff in spring is accelerating coastal ice melt (Whitefield et al., 2015), particularly along the extensive Arctic continental shelf, where the interplay between these variables influences the timing and magnitude of biological production (Arrigo and van Dijken, 2015;Marchese et al., 2017), the species composition of primary producers (Li et al., 2009;Ardyna et al., 2014), and, in turn, higher and lower trophic levels (Wassmann et al., 2011;Vernet et al., 2017). Establishing the baseline relationship between microbial communities in Arctic coastal waters and their physical and chemical environment is key to understanding and predicting how they will respond to continued climate-induced changes to the Arctic system.
Most investigations of seasonality in microbial community composition and function in the Arctic Ocean have focused on offshore regions in the Chukchi and Canadian Beaufort Seas, the Norwegian Coast, and the plumes of very large Arctic rivers (Alonso-Sáez et al., 2008;Garneau et al., 2008;Ghiglione et al., 2012;Marquardt et al., 2016;Onda et al., 2017). Less is known about shallow estuarine environments on Arctic coastlines, despite their importance to coastal fisheries (von Biela et al., 2013) and as breeding habitat for over 157 species of migrating birds . Nearly one-half the Alaskan Beaufort Sea coast and one-third of the Chukchi Sea coast is skirted by an irregular and discontinuous chain of barrier islands that enclose shallow (< 6 m deep) lagoons (Dunton et al., 2006;Schreiner et al., 2013). Seasonal changes in these lagoons are different than in the open Arctic Ocean. For example, the magnitude of seasonal temperature fluctuations is larger in the lagoons, ranging from as low as −2.1°C in the winter to over 10°C in the summer, while that in much of the rest of the Arctic Ocean does not exceed 0-4°C (Timmermans and Ladd, 2018). Salinity fluctuations are also larger in the lagoons, in some cases ranging from hypersaline in winter due to sea ice brine rejection to nearly fresh conditions in spring due to river inputs (Harris et al., 2017). The organisms inhabiting these lagoon systems must be capable of surviving rapid changes in physical and chemical conditions.
Several studies have demonstrated that organic carbon from terrestrial runoff subsidizes lagoon food webs in the Arctic (Dunton et al., 2006(Dunton et al., , 2012Bell et al., 2016;Mohan et al., 2016;Harris et al., 2018). These subsidies likely enter food webs via heterotrophic bacterial and protistan communities; however, the extent to which terrestrial subsidies influence the composition of microbial communities in these lagoons remains unknown. One study in a lagoon near Barrow, Alaska, used experimental incubations to show a change in Arctic marine bacterial community composition and an increase in production in response to tundra-derived organic matter amendments (Sipler et al., 2017). Understanding how coastal microbial populations incorporate terrestrial organic matter and use terrestrially derived nutrients is paramount to refining our understanding of pathways for the integration of terrestrial carbon into coastal Arctic marine systems. A first step in achieving this is to characterize how microbial populations in terrestrially influenced Arctic waters change seasonally and in response to inputs of riverine material.
In this study, we describe seasonal variation in prokaryotic and protistan community composition in coastal lagoons of the Alaskan Arctic Ocean, and identify potential controls on microbial population dynamics, including organic matter source and prokaryotic-eukaryotic associations. This work was carried out in the context of a larger interdisciplinary study aimed at understanding how terrestrial inputs control physical (Harris et al., 2017), biogeochemical (Connelly et al., 2015;Mohan et al., 2016), and ecological (Nolan et al., 2011;Dunton et al., 2012;Harris et al., 2018) properties of lagoon ecosystems along the Alaskan Beaufort Sea coast.

Sample Collection
Water samples (2 L) for microbial community analyses were collected from several sites within lagoons and outside barrier islands along the Alaskan Beaufort Sea coast in August 2011, and April, June, and August 2012, Jago (JA), Angun (AN), and Nuvagapak (NU), and one site outside the barrier islands near Barter Island (BP) were sampled in all three seasons (Figure 1, BP was not sampled in August 2011). Two more lagoons, Tapkaurak (TA) and Demarcation Bay (DE), and three additional sites outside the barrier islands, near the Hulahula River (HU), Bernard Spit (BE), and Demarcation Point (DP), were also sampled in August (Figure 1). Severe weather limited sample collection to KA, JA, AN, BP, and BE in August 2013. Samples were collected from one to two stations per site in April and June, and two to three stations per site in August of each year. BP had only one station in all seasons. Most sites were less than 4 m deep, with the exception of BE and DP, which were ~9-10 m deep. Samples were collected approximately 10 cm below the bottom off the ice cover in April (ice thickness 1.3-1.7 m) and from the top 2 m of the water column in June using a peristaltic pump, and by submerging hand-held sample bottles to ~0.5 m below the water surface in August. River endmembers were collected from the Canning, Jago, and HulaHula rivers in August 2011 and from Canning and Jago rivers in August 2012.
Samples were also collected for a suite of environmental measurements including particulate organic carbon (POC) and nitrogen (PON) concentrations and stable isotope ratios (POC δ 13 C and PN δ 15 N), chlorophyll a (Chl a) concentration, dissolved organic carbon and nitrogen concentrations, dissolved inorganic nitrogen (DIN = NO 3 + NH 4 ) concentrations, and oxygen stable isotope ratios of water (H 2 O-δ 18 O). Sample processing methods and measurements of particulate parameters (POC, PON, Chl a) are discussed in Connelly et al. (2015). Methods for dissolved parameters follow procedures described in McClelland et al. (2014). A YSI Sonde was used for measuring temperature, salinity, and dissolved oxygen from depths sampled (in addition to other depths throughout the water column). See Harris et al. (2017) for details of physical measurements and oxygen stable isotope ratios.

Microbial Sample Processing, DNA Extraction, and Polymerase Chain Reaction Amplification
After collection, samples were kept under shade during transit back to the Arctic National Wildlife Refuge field station in Kaktovik, Alaska. Within hours of collection, 2 L of water was filtered onto a 0.22-μm Sterivex filter (Millipore) using a peristaltic pump and preserved with 1 ml of DNA extraction buffer (100 mM Tris, 100 mM NaEDTA, 100 mM phosphate buffer, 1.5 M NaCl, 1% CTAB) and kept frozen until extraction. Prior to filtration, duplicate 14-ml samples were collected from the sample bottles, fixed with glutaraldehyde (2% final concentration), and frozen for estimation of bacterial abundance using flow cytometry.
Prior to extraction, Sterivex filter cartridges were cracked open with pliers and filters were removed using an ethanolflamed scalpel. The DNA extraction buffer from the cartridge was decanted into a sterile 2-ml microcentrifuge tube and the filter was subsequently cut into multiple pieces on a sterile cutting board and placed in the same tube. Samples were then subjected to three freeze-thaw cycles, followed by enzymatic lysis with Lysozyme (0.2 mg/ml final concentration) and Proteinase K (2 mg/ml final concentration) at 37°C for 30 min and continued digestion and lysis with the addition of SDS (1% final concentration) at 65°C for up to 2 h. Samples were then extracted two times with an equal volume of Phenol:Chloroform:Isoamyl alcohol (25:24:1) and nucleic acids were precipitated using 100% isopropanol (0.6 × volumes of the resulting supernatant) for 2 h up to overnight. Samples were then pelleted at 18,000 RCF for 30 min, rinsed, and re-pelleted two times with 70% ethanol, and dried down in a roto-evaporator. Once dry, samples were resuspended in 250 ml of nuclease-free water.
For community composition analysis, we amplified the V4 region (515F, GTGCCAGCMGCCGCGGTAA and 806R, GGACTACHVGGGTWTCTAAT) of the 16S rRNA gene for prokaryotic composition, and the V9 region (1391F, GTACACA CCGCCCGTC and EukBr, TGATCCTTCTGCAGGTTCA CCTAC) of the 18S rRNA gene for eukaryotic composition for sequencing on the Illumina MiSeq platform using Earth Microbiome Project protocols (http://www.earthmicrobiome. org/protocols-and-standards/16s/, but with only 30 PCR cycles). However, a known mismatch in the 16S primers with Thaumarcheaota, a dominant phylum of the marine Archaeal community, precluded us from drawing conclusions about Archaeal community composition. Each sample was amplified Frontiers in Microbiology | www.frontiersin.org three times, pooled, quantified using Picogreen, and then, for each amplicon, pooled at equimolar concentrations (100 μmol each). The 16S sample pool and 18S sample pool were each cleaned using a MoBio Ultraclean PCR Clean-Up Kit and quantified using Picogreen. Amplicon pools were sequenced at Argonne National Lab (the 16S sample library composed of August 2011, and April and June 2012 samples) or the Oregon State University Center for Genome Research and Biocomputing (all 18S sample libraries and an additional 16S library from August 2012 and all 2013 samples) 2 × 150 bp paired-end reads. Gene amplicon sequences have been deposited in NCBI Sequence Read Archive (SRA) bioproject accession number PRJNA530074, under run accessions SRR8832739-SRR8833063 (16S rRNA gene) and SRR8837972-SRR8838296 (18S rRNA gene) 1 .

Bacterial Abundance Measurements
Cell counts were performed using a BD Biosciences FACSCalibur Flow Cytometer at UMCES Horn Point Laboratory (2011 and 2012 samples) and Oregon State University (2013 samples). Single samples were counted for 2011 sites, while duplicate samples were counted and averaged for all sites after 2011. In the field, 14 ml of seawater was preserved with glutaraldehyde (2% final concentration) and frozen. In the lab, samples were thawed and 1.5-ml aliquots were stained overnight in the case of 2011 and 2012 samples with 20 μl of 1:200 SYBR Green I. The next day, samples were spiked with 15 μl (25 μl in 2013) of a sonicated beadstock created from PeakFlow Flow Cytometry Reference Beads (Life Technologies, Inc.) for internal reference. Samples from 2013 were stained and counted on the same day. Data were collected using the program CellQuest Pro (BD Biosciences) in logarithmic mode based on side scatter (SSC) and green fluorescence (FL1) with a target rate of 100-1,000 events s −1 for a total of 20,000 events for 2011-2012 samples and for a set period of time for 2013 samples (average 78,000 events). See Meyer et al. (2014) for additional methodological details, including how cell concentration was calculated from counted events.

Sequence Analysis
Reads that were successfully paired using fastq-join (Aronesty, 2011) were quality filtered with an expected error rate of 0.5, dereplicated (derep_fulllength), and abundance sorted (sortbysize) using UPARSE v. 8 (fastq_filter; Edgar, 2013). Singleton sequences were removed in the latter step to prevent them from seeding clusters when clustering sequences into operational taxonomic units (OTUs). Reads were then clustered into OTUs (cluster_otus in UPARSE pipeline) at 97% similarity. A de novo chimera check is inherent in the cluster_otus algorithm and chimeric sequences were removed during OTU clustering. Reference-based chimera filtering was performed using UPARSE (uchime_ref) with the Gold Database 2 as reference. Reads (including singletons) were subsequently mapped back to OTUs using UPARSE (usearch_global) and an OTU table created. Taxonomy of the representative sequences was assigned in QIIME v. 1.9 (assign_taxonomy.py; Caporaso et al., 2010) using the RDP classifier trained to the Greengenes database (v. 13.8, http://greengenes.secondgenome. com/) for 16S amplicons or the Silva database (v. 119; Quast et al., 2013;Yilmaz et al., 2014) for 18S amplicons. Any remaining singletons and OTUs occurring in only one sample were removed in QIIME (filter_otus_from_otu_table.py). Sequences identified as Archaeal, chloroplast, and mitochondrial were also removed from 16S reads. For the 18S rRNA gene library, we removed clades known to have multicellularity, as well as unclassified reads, in order to focus on protists. After these quality control steps, the average number of reads per sample was 22,326 for 16S amplicons (range 3,651-73,169 sequences per sample) and 43,093 sequences for 18S amplicons (range 6,720-103,750 sequences per sample).

Statistical Analyses
Given recent insights that rarefying microbiome datasets may not be the best method for comparing samples (McMurdie and Holmes, 2014), we chose not to randomly subsample OTU tables for the bulk of our analyses, with the exception of alpha diversity estimates. For alpha diversity measurements, the 18S rRNA gene OTU table was rarefied to 6,700 sequences per sample, and the 16S rRNA gene OTU table to 3,650 sequences per sample. Alpha diversity was calculated as Chao1 Diversity Index to measure species richness (Chao, 1984), Simpson's Evenness Measure (Smith and Wilson, 1996) to measure evenness, and Phylogenetic Diversity, which incorporates phylogenetic differences among species in the calculation of diversity (Faith, 1992;Caporaso et al., 2011). For beta diversity analyses, comparisons with environmental data, and indicator species analysis, OTU tables were normalized using proportional abundance of each OTU within each sample. To verify that using proportional abundance did not substantially change our conclusions compared to using OTU tables that were subsampled, we ran a subset of the analyses described in this paper using rarefied OTU tables and found no significant difference in results or conclusions.
Microbial community structure was assessed using nonmetric multidimensional scaling (NMDS) calculated using the metaMDS function in the Vegan package for R (Oksanen et al., 2019). Variability in bacterial and eukaryotic community composition among samples was calculated using Bray-Curtis dissimilarity. Permutational multivariate analysis of variance (PERMANOVA; Anderson, 2017) and Analysis of Similarity (ANOSIM; Clarke, 1993) calculated using the adonis and anosim functions in the Vegan package for R (Oksanen et al., 2019) were used to test for differences among sample groupings determined a priori (e.g., by season, inside versus outside of barrier islands). PERMANOVA provides a pseudo-F-ratio, a value of p for the group-wise tests for differences (as you would get from a standard ANOVA), and the percent of variation in the community dataset explained by the grouping. ANOSIM provides an R value ranging from 0 to 1 with higher values indicating stronger differences between or among groups, and a significance value for the ANOSIM R value based on 999 permutations.
The degree to which physico-chemical data explained the variation in bacterial and eukaryotic communities was assessed Frontiers in Microbiology | www.frontiersin.org using three methods. First, a Procrustes analysis was used to compare ordinations of community and physico-chemical data (Peres-Neto and Jackson, 2001) yielding correlation and significance values. Second, envfit in the Vegan package of R was used to decipher which variables were contributing to the structure of community nonmetric multidimensional scaling ordinations by fitting vectors of significant physico-chemical variables onto community NMDS ordinations. Finally, redundancy analysis (RDA) was used to quantify the percent of variation in bacterial or eukaryotic community composition explained by the physico-chemical environmental characteristics. Bacterial and Eukaryotic OTU tables were Hellinger-transformed prior to use in the RDA. Before running the RDA, physico-chemical variables for the model were selected to reduce multicollinearity using correlation matrices. The absence of substantial multicollinearity in this subset of variables was verified using the vif.cca function available in the Vegan package for R. RDA was run using the Vegan package for R.
Indicator species analysis (Dufrêne and Legendre, 1997) was used to identify bacterial and eukaryotic taxa that significantly contributed to seasonal differences in the coastal Beaufort Sea microbial community. In order to distinguish between river indicator species and lagoon indicator species, four sample groups were used for this analysis: River, April, June, and August. The indval program in the labdsv package for R was used to run the Dufrêne-Legendre Indicator Species Analysis, and OTUs having an indicator value (IV) > 0.7 and p < 0.005 were considered significant indicators. Monthly indicators were then further broken down into two groups, high-abundance indicators, having an average relative abundance of greater than 0.5% of the total average population for that month, and low-abundance indicators, which were significant but had an average relative abundance of less than 0.5%. The taxonomic composition of only high-abundance indicators was further scrutinized. The relationship between the distribution of highabundance indicators and the Beaufort Sea environment was examined using Spearman correlations, with values of p adjusted using the Benjamini-Hochberg correction (Benjamini and Hochberg, 1995). Correlations were calculated using the Hmisc package for R, while the calculated values of p were adjusted using the base R stats package.

Co-occurrence Network Analysis
Microbial association networks were generated for each month, across all years, using CoNet (Faust et al., 2012). In order for an OTU to be included in the network it had to be present in 25-33% of the samples (April minocc = 4, June minocc = 5, and August minocc = 14). In April and June, a percentage slightly higher than 25% was used because the number of correlations was very large and computation time was too great, preventing completion of network calculations at a minimum occurrence of 25%. Pairwise scores were computed for both Bray-Curtis similarity and Spearman correlation. Associations with a Spearman correlation above 0.7 or below −0.7 and a Bray-Curtis similarity of above 0.6 or below 0.4 were retained. For each measure and edge, 1,000 permutations (with renormalization for correlation measures) and bootstrap scores were generated, following the ReBoot routine. Values of p were calculated as described in Weiss et al. (2016) and measure-specific values of p were merged using Brown's method. Associations were corrected using the Benjamini-Hochberg's false discovery rate (Benjamini and Hochberg, 1995) and edges with merged values of p below 0.05 were retained. Edges had to be significant using both similarity measures to be kept. Network statistics were calculated in Cytoscape 3.6.1 (Smoot et al., 2011). Chord diagrams, created using the R package circlize, were used to display significant associations among the 15 most abundant taxa groups across all three networks (Gu et al., 2014).

Environmental Conditions
April waters were ice-covered and cold (average of −2°C; Harris et al., 2017), with high salinity and inorganic nutrients, and low Chl a, dissolved oxygen, pH, organic matter, and bacterial abundance (Figure 2A). June waters, sampled during ice breakup, were also cold but had the highest organic matter concentrations, the highest SUVA 254 (a measure of DOC aromaticity), and the lowest salinity because of freshwater input from rivers ( Table 1). August waters were warmer (average of 8.9°C), with lower concentrations of inorganic nutrients and organic matter, and higher values of H 2 O-δ 18 O, POC δ 13 C, and PN δ 15 N. Ranges of these variables fluctuated interannually, but seasonal patterns of change in coastal Beaufort Sea waters were the same from year to year.

Alpha Diversity
We identified 17,340 bacterial OTUs and 9,583 protistan (unicellular eukaryotes, including fungi) OTUs. For bacterial OTUs, river samples had the highest species richness and phylogenetic diversity (FDR-corrected p < 0.005). Evenness in river samples was not significantly different from lagoon coastal waters in April and August. Among lagoon and coastal samples, richness was highest in June (FDR-corrected p < 0.01, Supplementary Figure S1), phylogenetic diversity was lowest in August (FDR-corrected p < 0.01), and evenness was highest in August (FDR-corrected p < 0.01) and lowest in June (FDR-corrected p < 0.005). There was no interannual variability in richness and evenness in April or June, but evenness was significantly greater in August 2012 and 2013 than in 2011 (FDR-corrected p < 0.05, Supplementary Figure S2). Bacterial richness and evenness were the same between sites within and outside of barrier islands except in August when evenness outside the barrier islands was lower (FDR-corrected p = 0.0035, Supplementary Figure S3).
As with bacteria, eukaryotic species richness was greatest in rivers, but unlike bacteria the coastal eukaryotic communities had the lowest richness and phylogenetic diversity in June and the highest in April (Supplementary Figure S1). No significant differences in richness were observed among the eukaryotic communities when grouped by month or location.
There was also no interannual variability in eukaryotic richness, phylogenetic diversity, or evenness except in August 2011 when richness and phylogenetic diversity values were significantly lower than later years (FDR-adjusted p < 0.05, Supplementary Figure S2). Eukaryotic richness and evenness were the same between sites within and outside barrier islands in June, but richness was greater outside the islands in April (FDR-adjusted p < 0.05, Supplementary Figure S3) and evenness was greater inside the islands in August (FDR-adjusted p < =0.021, Supplementary Figure S3).

Beta Diversity
With respect to beta diversity, both bacterial and eukaryotic communities differed seasonally, and all coastal communities differed from river communities (Figures 2B,C). We found that 43% of the variance in bacterial communities and 27% Frontiers in Microbiology | www.frontiersin.org of the variance in eukaryotic communities was accounted for by seasonal coastal and river group differences (p = 0.001, using PERMANOVA). Pairwise seasonal differences, as quantified using ANOSIM tests, were greater among eukaryotic communities (0.62-0.97 for eukaryotes and 0.4-0.96 for bacteria, p = 0.001) with the exception of April-June comparisons, which had a higher ANOSIM R value for bacteria than for eukaryotes (0.93 vs. 0.88, p = 0.001). In total, 1,207 bacterial OTUs and 712 eukaryotic OTUs were shared among all coastal microbial communities. Among coastal communities, June had the greatest number of unique bacterial OTUs (49.8%) but the fewest unique eukaryotic OTUs (38.9%). August had the most unique eukaryotic OTUs (55%). June coastal communities shared the greatest percentage of OTUs with river communities among all coastal water-river water comparisons (38.8% for bacteria and 50.9% for eukaryotes). August communities also shared a substantial percentage of their OTUs with river communities (34.1% for bacteria and 31.7% for eukaryotes), but April samples shared 20% or less of their OTUs with river communities.

Indicator Species
Indicator species analysis was used to determine which OTUs significantly contributed to differences among seasons. We focused on OTUs with indicator values greater than 0.7, values of p less than 0.001, and average relative abundance greater than 0.5% of the average community composition for the month for which the OTU was an indicator. In April, 23 bacterial indicator taxa made up 36% of communities in these ice-covered waters ( Figure 4A). Many of these indicator taxa belonged to the Gammaproteobacteria order Oceanospirillales and the Bacteroidetes order Flavobacteriales. April indicator taxa also included members of the methanotrophic order Methylococcales, and several chemolithoautotrophic taxa including iron-oxidizing Zetaproteobacteria and sulfur-oxidizing SAR324. In June, 13 indicator taxa composed 50% of communities in these highly productive, lower salinity waters ( Figure 4B). Most of these taxa belonged to the Bacteroidetes phylum, and to the Betaproteobacteria genus Polaromonas and methylotrophic order Methylophilales. June indicators also included Alphaproteobacteria related to Loktanella sp. In August, 22 indicator taxa made up 31% of the bacterial communities in these late summer, nutrientpoor waters ( Figure 4C). Many of these indicators belonged to the Alphaproteobacteria family Rhodobacteraceae, including Phaeobacter and Octadecabacter spp. August indicators also included Gammaproteobacteria from oligotrophic marine clades (OM60, OM182, SAR86, and SAR92). As with the bacteria, a small number of eukaryotic indicator taxa made up a large fraction of the average April, June, and August communities (Figure 4). In April, 13 eukaryotic indicator taxa from four phyla made up 31% of the April eukaryotic community ( Figure 4D) and included several OTUs closely related to the parasitic order Syndiniales, and several marine stramenopiles (MAST) belonging to groups 1, 7, and 8. In June, a more diverse set of 17 indicator taxa made up 58% of the eukaryotic community ( Figure 4E). June indicators were dominated by several diatoms closely related to Chaetoceros, Skeletonema, and Melosira sp., but also included a diverse community of taxa from the class Dinophyceae (dinoflagellates), phyla Ciliophora (ciliates), Chlorophyta (green algae), and Cercozoa. In August, an even broader array of indicators was observed, with 23 OTUs representing 32% of the eukaryotic community ( Figure 4F). August indicators were dominated by Chlorophyta, Dinophyceae, and Ochrophyta OTUs, but included other taxa such as haptophytes, cryptophytes, and ciliates. While diatoms made up the bulk of the Ochrophyta indicators in June, this was not the case in August, when Ochrophyta indicators instead belonged to the Dictyochophyceae, Chrysophyceae, and Pelagophyceae.

Environmental Drivers of Coastal Beaufort Sea Microbial Communities
Several methods were used to investigate relationships between physico-chemical parameters and microbial community composition. First, Procrustes analysis showed that both bacterial and eukaryotic community composition were significantly correlated with variation in Beaufort Sea lagoon environmental conditions (Corr BAC = 0.7278, Corr EUK = 0.5923, sig = 0.001). Second, physico-chemical vectors that correlated significantly with bacterial and eukaryotic NMDS ordinations were overlain onto NMDS plots to determine the environmental gradient that correlated with the variations in community structure (Supplementary Figure S10). For the bacterial community, the first NMDS axis was negatively correlated with salinity and POC δ 13 C (r 2 < −0.7) and positively correlated with SUVA 254 , POC, DOC, and Chl a suggesting that separation of June communities from August and April along this NMDS axis represents gradients in terrestrial input and productivity. The second NMDS axis was strongly correlated with temperature (r 2 < −0.9) and negatively correlated with nitrate and ammonium (r 2 > 0.9, Supplementary Figure S10A), suggesting that separation of August from April communities is driven by seasonal changes in temperature and nutrients. Finally, redundancy analysis was used to quantify the amount of variation explained by these physico-chemical variables (Supplementary Figure S11). Approximately 70% of the variation in bacterial community composition could be explained by the Beaufort Sea environment.
For the eukaryotic community, the orientation of samples on the NMDS was rotated slightly relative to the bacterial NMDS, but the major trends were essentially the same. April and August communities were separated from June communities along gradients in terrestrial input and productivity (Supplementary Figure S10B). The C:N ratio of particulate organic matter (POM) and concentrations of DOC, nitrate, and ammonium were all strongly correlated this first NMDS axis (r 2 > 0.7), suggesting that the August and April communities existed in waters with more degraded organic matter and lower nutrients than June communities. Like bacterial communities, April and August eukaryotic communities separated along a temperature and nutrient gradient, but productivity (e.g., Chl a) and organic matter source components (SUVA 254 , Chl a, Salinity, and POC δ 13 C) were also important correlates with the second NMDS axis for the eukaryotic ordination (Supplementary Figure S10B). Finally, redundancy analysis showed that approximately 55% of the variation in eukaryotic community composition could be explained by the Beaufort Sea environment (Supplementary Figure S11).
Correlations between high-abundance indicator taxa and measured physico-chemical variables reflected those of the whole communities ( Figure 5). April indicator taxa in ice-covered waters correlated with less productive and cold, nutrient-rich conditions (negative correlations POC, PN, Chl a, bacterial

Microbial Co-occurrence Network Properties
Network Structure Co-occurrence networks for the three seasons were strikingly different in size and topology ( Table 2). The percentage of bacterial and eukaryotic OTUs retained in these networks after frequency filtering and subsequent co-occurrence analyses varied from 7 to 36% (Supplementary Table S1), causing the networks to vary substantially in size. Among the taxa retained in these networks, significant combinations represented 1.4-14.9% of the possible combinations (Supplementary Table S1), and average path length, or the number of nodes needed to link two nodes was short, ranging from 3.1 to 4.3 ( Table 2). Network diameter, or the longest distance in a network, ranged from 10 edges (June) to 13 edges (August). Network density, a normalized measure for the average connectivity within a network, was the highest in June (0.057) and lowest in August (0.014). The average node degree, or the average number of connections for each node, was very high in June (111) compared to April (33) and August (9). Overall network complexity, as estimated by connectance (the fraction of all possible links that are realized in a network; Williams et al., 2002), was highest in June (0.0282) and lowest in August (0.0071) ( Table 2). The number of connected components, or a set of nodes in the network graph for which there is always an interconnecting path (Corel et al., 2016), was the lowest in April (5) and highest in August (36, Table 2). If a network has only one connected component, all nodes can be linked to any other node in the network either directly or indirectly. The presence of more than one connected component indicates that some groups of nodes are segregated from the main network, not significantly correlated with any nodes in that graph.  Indicator taxa were most abundant in their corresponding seasonal networks, and river indicators were most abundant in the June network (Supplementary Figure S12). In all cases, these indicator taxa were outnumbered by non-indicator taxa, but in most cases indicator taxa had higher average node degrees than non-indicator taxa (Supplementary Figure S13), indicating that they were more highly connected within the networks. This was the case for seasonal indicators in the April and August networks, and for river indicator taxa in the June network underscoring the importance of river microbes in spring surface waters in these coastal lagoons.

Network Composition
Taxa that were abundant during each season were also abundant in corresponding seasonal networks (Supplementary Figure S14). The April and August networks shared many of the same abundant taxa (Supplementary Figure S16) including bacterial taxa Rhodobacterales, Alteromonadales, Oceanospirillales, and Flavobacteriales, and eukaryotic taxa Rhizaria, Diatomea, Ciliophora, Syndiniales, and Dinophyceae. The June network included some of the same abundant eukaryotic taxa with the addition of Chrysophyceae, but featured a different set of bacterial taxa including Legionellales and Actinobacteria (Supplementary Figure S16).
All three networks were dominated by positive edges (indicating co-presence), with far fewer negative edges (indicating mutual exclusion; Table 1). Also, organism-organism associations were far more abundant than those between organisms and environmental variables (Supplementary Figure S17). The June network was dominated by co-occurrences between prokaryotes (Figure 6), while April and August networks had a more even distribution of bacteria-bacteria (Bac-Bac), bacteria-eukaryote (Bac-Euk), and eukaryote-eukaryote (Euk-Euk) associations (Supplementary Figure S18).
We examined the positive and negative correlations among the 10 most abundant taxonomic groups within each network, which represented a total of 15 different groups (Figure 6). In April, Euk-Euk edges were dominated by Syndiniales taxa, including marine alveolate (MALV) Groups I and II and Amoebophrya. These taxa correlated most frequently with themselves, with dinoflagellates including Gymnodiniphycidae and Peridiniphycidae, and with ciliates, including Oligotrichia and Choreotrichia (Figure 6). Syndiniales also co-occurred with more bacterial taxa than other abundant eukaryotic taxa, and were positively correlated with Flavobacteriales and Deltaproteobacteria.
There were far fewer Euk-Euk and Bac-Euk edges than Bac-Bac edges in June, and the most abundant eukaryote in the June network, Chrysophyceae, was not the most abundant eukaryotic taxa in June samples (i.e., Diatomea). Chrysophyceae correlated mostly with themselves but also had significant correlations with Betaproteobacteria, Deltaproteobacteria, and Legionellales (Figure 6). By comparison, Diatomea had far fewer correlations within the network. The June network was dominated by Bac-Bac edges involving Betaproteobacteria, Legionellales, and Rhizobiales, all of which are commonly associated with freshwater and brackish environments (Figure 6). Deltaproteobacteria also had a high number of edges in June networks.
The August network featured many of the same taxonomic groups as the April network, but the connections among the nodes were different. Most Euk-Euk associations were positive correlations among Syndiniales and other dinoflagellates. By comparison, correlations involving ciliates were less frequent. Many Bac-Euk associations were negative, particularly those involving protist groups Syndiniales and Rhizaria and bacterial groups Alteromonadales, Flavobacteriales, Betaproteobacteria, and Rhodobacterales (Figure 6). In contrast, most Bac-Bac associations in the August network were positive.
The most connected taxa in the three networks were not always the most abundant taxa, suggesting that, in many cases, more abundant taxa may not require mutualistic interactions to thrive and can become abundant without the "help" of other microbial taxa, while the reverse may be true for the less abundant but highly connected taxa (Supplementary Figure S15). In the April network, the highest node degrees were associated with relatively rare bacterial OTUs related to the gammaproteobacterium HTCC2188, Thiotrichales, and Gemmatimonadetes, and eukaryotic taxa Developayella and MAST-3 OTUs, a single Goniomonas OTU and a Pirsonia OTU (Supplementary Figure S19). Similarly, in the June network, taxa with the highest node degrees were relatively rare OTUs from the Enterobacteriales, WS3, and SR1 taxa, all with average node degrees >300 (Supplementary Figure S19). However, OTUs representing the abundant taxa Rhizobiales and Betaproteobacteria also had high average node degrees (>200). In the August network, many of the abundant taxa had high average node degrees, including Oceanospirillales and Alteromonadales (Supplementary Figure S19). Also, in August, several taxa with high average node degrees featured a large fraction of negative correlations including the bacteria Saprospirales, and eukaryotes MAST-9 and Palpitomonas (Supplementary Figure S19).

DISCUSSION
Coastal waters along the North Slope of Alaska are important feeding and breeding grounds for many species of migratory birds Taylor et al., 2010) and fish, such as Arctic char and Arctic cod (Craig, 1984) that are critical to native subsistence fisheries (Dunton et al., 2012). Maintenance of healthy lagoon and coastal ecosystems is crucial to sustaining these higher trophic levels. The base of food webs in these ecosystems is occupied by several interacting and species-rich microbial communities that perform many important ecosystem services, including organic matter degradation, nutrient regeneration, and carbon fixation via photosynthesis. In coastal systems, these communities provide a critical pathway for the incorporation of terrestrial organic matter and nutrients into estuarine and marine food webs (e.g., Carlsson et al., 1993;McCallister et al., 2004), especially in the Arctic, where terrestrial inputs are high (Whitefield et al., 2015). In the Beaufort coastal lagoons, microbial communities must maintain ecosystem functions despite huge seasonal changes in environmental conditions. This study demonstrates that microbial communities in these lagoons respond to seasonal changes through annually repeating seasonal shifts in species composition of both prokaryotic and microbial eukaryotic communities.

Terrestrial Subsidies
Previous studies have shown that terrestrial inputs of organic matter help fuel food webs along the Alaska Beaufort Sea coast (Dunton et al., 2006(Dunton et al., , 2012Harris et al., 2018); yet, no studies to this point have characterized relationships between the microbial communities living within these coastal waters and the organic matter inputs to them. Characterized by relatively low concentrations of POC, PON, and pigments, especially Chl a (Connelly et al., 2015), April waters in these coastal lagoons had a low contribution of phototrophic microbial taxa. The suspended organic material present was highly processed, with Arcs that remain within a bar denote significant correlations among OTUs within that taxon.
high ratios of phaeopigments to Chl a and elevated saturated fatty acid proportions (Connelly et al., 2015), suggestive of a heterotroph-dominated system. Indeed we show that April communities were dominated by high proportions of small heterotrophs (e.g., MAST) and parasitic Syndiniales clades in eukaryotic communities ( Figure 2B). Furthermore, the refractory nature and low concentrations of organic matter favored relatively high proportions of chemoautotrophs in bacterial communities (Figure 2A, Supplementary Figures S4, S5). OTUs belonging to family Oceanospirillaceae, members of which have been implicated in hydrocarbon degradation (Satomi and Fujii, 2014), were also in high relative abundance in April. Arctic peat contains hydrocarbons (Yunker et al., 1993) and aromatics that likely contribute to the DOM in these coastal waters, as observed further west in the coastal Chukchi Sea (Sipler et al., 2017). The ability to degrade what is commonly considered more refractory components of organic matter (Yunker et al., 2002) may give members of Oceanospirillaceae a competitive advantage at the end of winter, after all of the fresh phytoplankton-derived organic matter has been degraded. By June, the peak of the spring freshet had passed and ice break-up was well underway. POM analyses pointed to a much more productive system characterized by carbon inputs from both terrestrial sources and autochthonous phytoplankton (Connelly et al., 2015). Eukaryotic microbial communities were dominated by diatoms, and bacterial communities by a mixture of freshwater bacteria and a distinct estuarine community that presumably grew to dominate these communities in each year of the study (Figure 3B). Bacteroidetes, including Cyclobacteriaceae and Flavobacterium spp., and Betaproteobacteria, particularly Polaromonas sp. were abundant in the June surface waters of the lagoons. Cyclobacteriaceae and Flavobacterium sequences were found to be enriched in low-salinity waters of the Columbia River estuary in Oregon and generally showed tolerance to a wide range of salinities (Smith et al., 2017). Polaromonas is a euryhaline bacterial taxa that can survive across a wide range of salinities (Veillette et al., 2011). Interestingly, this taxa was observed to be enriched in sea ice brackish brines (salinity 2.4-9.6) in the central Arctic Ocean but not in the surface seawater below the sea ice (salinity 33.3-34.9), which was thought to indicate that they were unable to survive the salinity shock during brine rejection (Fernández-Gómez et al., 2019). However, closer to the coast we observed that Polaromonas appears to survive this transit from sea ice into surface seawater; perhaps lower salinity surface waters resulting from river input coincident with sea ice melt provide a refuge for these taxa.
By the middle of the open-water period in August, POM was characterized by elevated proportions of terrestrial and dinoflagellate fatty acids relative to those of diatoms (Connelly et al., 2015), which was validated by a shift from a diatomdominated community in June to a dinoflagellate-dominated community in August (Figure 3B, discussed in detail below). Coincident with these changes in OM source and decreased inorganic nutrient concentrations, the bacterial community came to resemble a typical coastal ocean community, becoming enriched in bacterial clades commonly considered to be oligotrophic, including SAR86, SAR92, and OM182 (Figure 4). Many of these clades were also present in April but in lower proportions. SAR92 and OM182 belong to the oligotrophic marine Gammaproteobacteria (OMG), while SAR86 is more distantly related and possesses an even more streamlined genome (Spring et al., 2015). OM182 and SAR86 were observed to become more abundant in late summer and fall in brackish waters of the Baltic Sea (Hugerth et al., 2015), aligned with an oligotrophic lifestyle. SAR92 is common in coastal waters at high (Ghiglione et al., 2012) and lower latitudes (Teeling et al., 2012), often in association with or just following phytoplankton blooms. Since Chl a was lower in August than in June, it appears that SAR92 can also persist in coastal waters of the Beaufort Sea well after peak primary production. Altogether, these parallels between our study and Connelly et al. (2015), coupled with the high percentage (55-70%) of community variation explained by physico-chemical measurements, suggest strong linkages between organic matter source and microbial community composition, and are consistent with similar seasonal changes in POM pigments and phytoplankton communities observed further east near the Mackenzie River plume (Morata et al., 2008).

Photosynthetic Protists
In much of the Arctic Ocean, diatoms are the most abundant primary producers in spring (Figure 5), while smaller picoeukaryotes dominate the photoautotroph community the remainder of the year Marquardt et al., 2016). This also occurs in the Beaufort coastal lagoons in June despite significant river influence and lower salinity. Chaetoceros and Thalassiosira, two abundant taxa in June, are dominant diatoms in under-ice blooms on the Chukchi Shelf (Arrigo et al., 2012) and in pelagic spring blooms across the Arctic (Poulin et al., 2011), including in the Beaufort Sea (Balzano et al., 2012b). Melosira and Navicula, also abundant in June, are common sea-ice associated diatom taxa (Booth and Horner, 1997;Poulin et al., 2011) that are thought to seed pelagic phytoplankton communities in spring (Michel et al., 1993;Hardge et al., 2017). In the Beaufort coastal lagoons, Melosira was only abundant in June, suggestive of a sea-ice source, but Navicula was present in all seasons demonstrating that, while sea-ice may be a source for Navicula, members of this genus persist in the water column and contribute to the pelagic phytoplankton community (Hardge et al., 2017). Other notable primary producers in spring included the chlorophytes Carteria, Chlamydomonas, and chrysophyte Ochromonas which are commonly considered to be freshwater and snow genera, but their presence has been reported in Arctic coastal waters influenced by the Mackenzie River (Balzano et al., 2012a) and in sea ice and melt ponds elsewhere in the Arctic (Kilias et al., 2014). Given that members of these genera appear to survive across a wide range of salinities, these euryhaline phototrophs may become increasingly important in coastal Arctic waters, and across the Arctic as a whole, with the forecasted freshening of the Arctic Ocean (McPhee et al., 2009;Morison et al., 2012).
As spring progressed into summer, the composition of the primary producers shifted from large cells (diatoms) to smaller picophytoplankton, predominately prasinophytes Micromonas and Bathycoccus (Supplementary Figure S9). Micromonas was dominant in 2011 (average 11% vs. 0.9% for Bathycoccus) when summertime waters were relatively cold (7.9°C) and salty (27.5 PSU). Bathycoccus was dominant in 2012 and 2013 (4.7 and 6.5% vs. 0.6 and 0.5% for Micromonas) when waters were warmer and fresher (9-11.2°C, 19.6-22.1 PSU). Micromonas is commonly thought to be the most abundant Arctic prasinophyte , but Bathycoccus was more abundant on the river-influenced Mackenzie Shelf (Monier et al., 2015) and during polar sunset and polar night in the Amundsen Gulf Region (Joli et al., 2017), possibly due to differences in low-light survival strategies. In both cases, light, grazing, and nutrients were hypothesized to drive this taxonomic shift. Both Micromonas and Bathycoccus have been shown to be capable of osmotrophy (Hernández-Ruiz et al., 2018), but Bathycoccus appears to have a stronger preference for amino acids as a carbon source relative to bicarbonate, suggesting that Bathycoccus is particularly adapted to organic matter utilization. Low-light adaptation and the ability to consume organic matter may explain the success of Bathycoccus in 2012 and 2013, when river input and organic matter concentrations were higher.
Under the ice in April, photoautotrophs were much less abundant, and consisted mainly of the prasinophytes discussed above, stramenopiles related to Bolidophyceae (2.7%) and dictyochophyte Pedinellales (3.4%), which was also abundant in August (4.1%). Bolidophyceae and Pedinellales have been observed elsewhere in the Beaufort Sea (Balzano et al., 2012a;Terrado et al., 2013), under the ice in the Central Arctic Ocean (Pedinellales only; Hardge et al., 2017), and in Canadian High Arctic sea ice (Piwosz et al., 2013). Seventeen percent of the bolidophyte cells investigated from sea ice were found to have at least one bacterium in their food vacuole (Piwosz et al., 2013) and thus their presence under ice may be sustained through heterotrophy rather than photosynthesis.

Chemolithoautotrophs and Methylotrophs
Ice-covered waters were dominated by bacterial taxa known to thrive in low-organic matter conditions, such as chemolithoautotrophs including Zetaproteobacteria (4.6%), Deltaproteobacteria clade SAR324 (1.3%), and methylotrophs including Methylococcales (2%) (Figure 3). Zetaproteobacteria (4.6%) are mat-forming Fe(II) oxidizers that are closely related to the chemotrophic iron-oxidizing genus Mariprofundus (Singer et al., 2011). Their presence in April waters was consistent with the orange tint observed on several April sample filters, but was surprising given that this is, to our knowledge, the first evidence of this microbial taxa in coastal Arctic waters. Zetaproteobacteria have been observed in iron-rich hydrothermal vents of the Loihi Seamount (McAllister et al., 2011), and coastal waters in Maine, USA (McBeth et al., 2011). Our study extends their distribution to include the coastal Arctic Ocean. These lagoons receive large pulses of iron during the spring snow melt (Rember and Trefry, 2004), and iron concentration in Arctic freshwaters increases through spring and summer (Pokrovsky et al., 2016) and may be enhanced by permafrost thaw (Barker et al., 2014).
SAR324 (1.3%) have genes for sulfur and alkane oxidation and have the capacity to degrade short-chain fatty acids, among other metabolic strategies (Sheik et al., 2014). Interestingly, Connelly et al. (2015) observed the highest proportional abundance of short-chain fatty acids in April waters compared to June or August. SAR324 were also found to be proportionally more abundant in surface waters under-ice than in open waters off Point Barrow (Sipler et al., 2017) and were shown to be important in nitrogen cycling in the winter (Connelly et al., 2014). Chemoautotrophic production under-ice may help sustain biological communities during the long winter, as it is thought to in other continually ice-covered systems (Boyd et al., 2014;Vick-Majors et al., 2016).
Methylococcales are exclusively methylotrophs and type I methanotrophs that have been observed to thrive in association with iron-oxidizing microbial mats in freshwater systems (Quaiser et al., 2014). Our data suggest that iron-oxidizing and methaneoxidizing bacteria may also live in close association in ironrich, coastal marine waters. Methane is present in shallow sediments throughout the Beaufort Sea shelf (Coffin et al., 2013), and dissolved methane is highly concentrated in Beaufort Sea water (Lorenson et al., 2016), particularly in shallow waters. This methane is mainly generated by microbial degradation of organic matter (Lorenson et al., 2016), but may also arise from permafrost-associated methane gas hydrates (Shakhova et al., 2017), which are present throughout the Beaufort shelf region (Riedel et al., 2014). We did not measure methane concentrations in our water samples, but wintertime under-ice methane concentrations to the west of our sample region were 3-28 times greater than in summer (Kvenvolden et al., 1993). While both of the iron-and methane-oxidizing taxa in our samples are aerobic, they prefer to live at oxic-anoxic interfaces to allow for the presence of both oxygen and reduced electron donors, which may have been available given the presence of low oxygen levels in some lagoons in April (Connelly et al., 2015). Overall, the presence of these bacterial functional groups suggests that iron, methane, nitrogen, and sulfur cycling become relatively important under the ice in these coastal lagoons as more labile organic matter is progressively depleted through the long Arctic winter.

Parasites
Heterotrophic protists play an important part in marine food webs as grazers of phytoplankton and bacterioplankton, and as food for zooplankton. In the Central Arctic Ocean, their biomass can rival or exceed that of phototrophic protists (Sherr et al., 1997). Heterotrophic and parasitic protists were relatively abundant in the Beaufort coastal lagoons in all seasons, but were particularly dominant in April when sea ice and snow attenuated light penetration into surface waters, limiting the abundance of photosynthetic protists (Figure 3). Thus, heterotrophy and parasitism likely dominated the protistan lifestyle in April waters. This is corroborated by the presence of abundant sequences related to ciliates, heterotrophic dinoflagellates, parasitic Syndiniales, and MAST taxa in the under-ice community (Figure 3). Heterotrophic protists were relatively less abundant in spring, but became dominant again later in summer, following the spring bloom and depletion of macronutrients (Table 1).
Frontiers in Microbiology | www.frontiersin.org Syndiniales, including the MALVs, are a globally distributed parasitic group within the Alveolates (Guillou et al., 2008;de Vargas et al., 2015), that constitute a substantial component of the global marine interactome (Lima-Mendez et al., 2015), and are generally considered to have a broad host range from other protists to metazoans (Guillou et al., 2008). We observed clear seasonality in Syndiniales, especially Groups I and II, with the greatest relative abundance under ice in April (20%), and lower relative abundance in June (4%) and August (6.3%). Syndiniales followed a similar abundance pattern in Franklin Bay further east in the Beaufort Sea (Terrado et al., 2009), in a high-Arctic Fjord (Marquardt et al., 2016), and in Antarctic waters (Cleary and Durbin, 2016). Oxygen may influence the distribution of these two groups in the water column, with Group I preferring low-oxygen waters or sediments and Group II preferring oxygenated waters (Guillou et al., 2008), although both groups were abundant in suboxic and anoxic fjord waters in British Columbia (Torres-Beltrán et al., 2018). Low dissolved oxygen was measured in some of the lagoons under the ice (Connelly et al., 2015), yet Group II were the most abundant Syndiniales in this season, which further suggests that oxygen is not the only driver of their distribution and that other factors such as host availability and host stress under winter conditions (Cleary and Durbin, 2016) control the abundance and diversity of Syndiniales.

Grazers
Heterotrophic flagellates like marine stramenopiles are ubiquitous in the global ocean (Lovejoy et al., 2006;de Vargas et al., 2015) and, as bacterivores, represent important links in marine microbial food webs, transferring carbon from bacteria to higher trophic levels like zooplankton (Monier et al., 2013;Worden et al., 2015). More abundant in the lagoons during less productive months, MAST clades MAST-1 and MAST-6 had the highest relative abundances in our dataset with MAST-1A and MAST-1C most abundant in April (6% of 18S rRNA genes). These clades were also found to be abundant in near-ice or under-ice stations along a transect from the Labrador Sea west to the Beaufort Sea . In August, MAST-6 was the most abundant MAST clade (2.2%). This clade has rarely been reported in the Arctic, but that is likely due to the fact that it is missed by the PCR primer set commonly used to assess protist diversity in Arctic waters (e.g., Thaler and Lovejoy, 2015). Using CARD-FISH, MAST-6 was found in first year sea ice in the Canadian Arctic Archipelago, with 20% of the cells containing at least one bacterium in their food vacuoles, suggesting that MAST-6, like MAST-1, are bacterivorous (Piwosz et al., 2013). In the Baltic Sea, MAST-6 cells were found to have both phytoplankton and bacteria in their food vacuoles suggesting that they are both algivorous and bacterivorous (Piwosz and Pernthaler, 2010). This clade of marine stramenopiles has been observed to prefer sediments across several coastal stations around Europe (Logares et al., 2012); however, our observations show that they are also important in pelagic systems in the Arctic.
Dinoflagellates had the highest relative abundance in August (24%), followed by April (14%) and June (10%); however, like with MAST, the dominant taxa varied by season. While both April and August were dominated by the Gyrodinium sp., Gymnodinium sp. was also abundant in August. While difficult to identify microscopically (Lovejoy, 2014;Kubiszyn and Wiktor, 2016), these two genera of naked heterotrophic dinoflagellates are abundant in 18S rRNA gene surveys of Arctic waters (Comeau et al., 2011;Marquardt et al., 2016). The dinoflagellate population in June was dominated by Pelagodinium sp. (7.6%), a member of the Suessiaceae. Pelagodinium is thought to be a symbiont of Foraminifera (Siano et al., 2010), but forams were a very small fraction of the protist communities, especially in June (<0.001%). The highest abundance of Foraminifera was observed in April (0.1%), but still was small compared to the relative abundance of Pelagodinium sp. in June. Given these observations, it is possible that we detected this symbiont during the free-living stage of its life cycle or that it is also a symbiont of other taxa abundant in June.
While typically less abundant than dinoflagellates in the Arctic, ciliates represent another important group of grazers in marine systems (e.g., Sherr et al., 1997). In line with microscopically obtained abundance estimates, dinoflagellate 18S rRNA gene sequences were always at least twice as abundant as ciliates regardless of season in a high-Arctic fjord (Marquardt et al., 2016). We observed more seasonality in the ratios of these two groups of protists, with ciliate sequences twice as abundant as dinoflagellate sequences in April and June, but less abundant in August (17%, compared to 24% dinoflagellates). This suggests that ciliates may play a more important role in coastal lagoon food webs than in other Arctic systems.
Oligotrich ciliates were the most abundant group of ciliates across all seasons, including Strombidium and Laboea (the latter only in August; 1.4%). OTUs classified as Strombidium were two times more abundant in June and August than April, while all months had a large percentage of reads that could not be classified beyond Oligotrichia, similar to other studies of ciliates in polar waters (Onda et al., 2017). Strombidium was found to be abundant in surface waters elsewhere in the Arctic, especially in the spring and summer, perhaps in part due to a mixotrophic lifestyle (Stoecker et al., 2017). Strombidium and Laboea have been observed to temporarily retain and gain energy from the chloroplasts from ingested diatoms or other phototrophic prey (Dolan and PÉrez, 2000). This could provide an energetic advantage over ciliates that rely solely on phagotrophy. April had higher relative abundances of ciliates belonging to the Mesodiniidae (3.1%) and Oligohymenophorea (specifically Scuticociliatia; 2.5%). Oligohymenophorea are strictly bacterivorous (Vaqué et al., 2008), but some Mesodiniidae species are mixotrophic, bordering autotrophic, with a preferred diet of cryptomonads as a source for harvested chloroplasts (McManus and Santoferrara, 2013). Cryptomonads were most abundant in April and August, in line with the distribution of Mesodiniidae OTUs. Mixotrophy was found to be the primary metabolism of ciliates in the oligotrophic waters of Fram Strait (Seuthe et al., 2011). The dominance of several potentially mixotrophic groups of ciliates in this study suggests that mixotrophy is also important in the coastal lagoons of the Beaufort Sea and could contribute to the overall productivity of these waters.

Community Connectivity and Microbial Food Web
We used co-occurrence network analysis to investigate prokaryotic and eukaryotic community connectedness in each season using data from multiple years. We were unable to find similar seasonal networks for comparison in marine systems because most marine networks have been grouped by depth (Lima-Mendez et al., 2015;Milici et al., 2016) or were generated for entire time series datasets without seasonal breakdown (e.g., Chow et al., 2013). Still, our network sizes and clustering coefficients ( Table 2) were within the range of these marine co-occurrence networks.
In one of the only aquatic microbial time series studies that performed seasonal network analysis, it was observed that network complexity of lake microbial communities was the greatest in spring, compared to summer or autumn (Kara et al., 2013). Similarly, we found that the June network was the largest, most connected, and most complex, having the highest clustering coefficient, highest connectance, and the lowest average path length, while the August network was the least complex and the complexity of under-ice network in April was intermediate. Kara et al. (2013) also noted that spring and autumn samples had the lowest and highest diversity, respectively, resulting in a negative relationship between network complexity and diversity. We did not observe this same relationship in the lagoons, possibly due to the contribution of freshwater bacteria and protists to the marine system and the formation of a diverse brackish microbial community in June. Given the same seasonal trends in network complexity between the lake study and our study, but differences in diversity-network complexity relationships, it is possible that seasonal trends in aquatic microbial network characteristics may be independent of the number of taxa present and driven more by ecosystem productivity, with food web complexity being highest during periods of high production.
River indicator taxa had the highest average node degree in June networks (Supplementary Figure S16) and common freshwater and brackish bacterial taxa accounted for the bulk of the significant correlations in the June networks (e.g., Legionellales, Betaproteobacteria, Deltaproteobacteria, and Rhizobiales, Figure 6). Whether these taxa were actively interacting or just happen to be passively co-existing is impossible to determine from this analysis, but it is important to note that river-impacted, nearshore systems may follow different diversity and connectivity patterns than lake or open ocean systems because mixing of freshwater and marine communities may elevate microbial diversity and form more complex microbial networks. The prevalence of freshwater and brackish microbial nodes in the spring network underscores the importance of these taxa in the coastal Arctic ecosystem during the spring freshest. Interestingly, two bacterial candidate phyla with small genomes, including WS3 and SR1 (Kantor et al., 2013;Farag et al., 2017), represented network "hubs" in June, with average node degrees >300. Microbes with streamlined genomes have also been observed to be network hubs in other studies of freshwater and marine systems (Peura et al., 2015;Milici et al., 2016), relying on interactions with other taxa in order to obtain metabolites that they cannot synthesize themselves.
Protists represented a small fraction of the nodes in the June networks, possibly due to the higher relative abundance of photoautotrophic protists during ice break-up. Most of the significant protistan relationships in June were between protists and bacteria rather than with other protists and were dominated by Chrysophytes (which can by mixotrophic; Beisser et al., 2017), diatoms, and heterotrophic groups of Rhizaria (Nakamura and Suzuki, 2015). These protists could be obtaining a necessary metabolite produced by co-located bacteria or, if heterotrophic or mixotrophic, could be grazing on bacteria, which can be enhanced during phytoplankton blooms (Hyun and Kim, 2003). Still it is important to note that the total number of significant protist-bacteria edges (11,256) in June exceeded those in April or August, but was far smaller than the number of Bac-Bac edges (94,983) and thus the relative contribution was less than in other months.
The April network was the next largest and contained the largest number (and percentage) of eukaryote nodes and Euk-Euk edges, especially among nodes belonging to the Syndiniales, other dinoflagellates, and ciliates ( Figure 6A, Supplementary  Figure S19). These relationships are in line with the known hosts for this group of parasites (Guillou et al., 2008;Torres-Beltrán et al., 2018) and support our hypothesis that parasitism was an important component of the under-ice food web. Extreme winter conditions in polar systems may increase parasitism due to environmental stress (e.g., low light, low in situ production). Syndiniales were also the most abundant group in co-occurrence networks generated as part of the Tara Oceans project from all regions sampled except for the Southern Ocean (the Arctic Ocean was not sampled as part of Tara Oceans; Lima-Mendez et al., 2015). As was observed at lower latitudes, these parasitic OTUs were most commonly correlated with other Syndiniales OTUs and with Dinophyceae OTUs. Syndiniales also correlated with several radiolarians, consistent with direct observations of similar associations through single-cell sequencing of radiolarians from a Norwegian fjord (Bråte et al., 2012), and with correlations between these protist groups under the ice north of Svalbard (Meshram et al., 2017). But unlike at lower latitudes (Lima-Mendez et al., 2015), correlations between Syndiniales and Ciliophora OTUs were common under the ice in the April. It is not possible to determine if this means that parasitism of ciliates is more prevalent in Arctic waters than at lower latitudes, but these observations suggest that correlations between putatively parasitic OTUs and presumed hosts under sea ice warrant continued investigation.
The August microbial community was less connected and more fragmented than other seasons, with a higher number of connected components, highest average path length (number of nodes needed to link individual nodes), and the lowest connectance (fraction of all possible links that are realized). A similar pattern was found for a Wisconsin lake in which the autumn network had the smallest network size and highest path length (Kara et al., 2013). The breakdown of microbial networks between June and August may be driven by physical processes such as increased mixing and reduced water column stratification, which was weaker in August than in June (Harris et al., 2017). Another possible explanation is that the August food web is not as reliant on fast energy transfer from large, fast-growing phytoplankton but rather on a slower energy transfer characteristic of a more detrital food web, the latter of which is often characterized by longer path lengths and generally weaker links (Rooney and McCann, 2012).
The August network featured the greatest percentage of negative relationships between nodes, suggesting more antagonistic relationships in the microbial food web in the summer than the spring (Figure 6). This is supported by our FIGURE 7 | Diagram depicting major seasonal microbial community changes in the Beaufort Lagoons ecosystem, with arrows depicting hypothesized pathways for carbon flow from one group of organisms to another. The size of the arrow indicates relative magnitude of this hypothesized carbon flow, and the size of the text indicates relative size of the carbon pool, based on the relative abundance of microbial taxa each month. Many symbols after  or courtesy of the Integration and Application Network (ian.umces.edu/symbols/). observation of high relative abundance of heterotrophic protist sequences in August. Furthermore, the concentration of phaeopigments was the highest in August (Connelly et al., 2015), which also supports high grazing on and microbial remineralization of phytoplankton-derived POM. If microbial populations were more focused on degrading or grazing on phytodetritus in August, this could result in fewer significant correlations among taxa and also longer chains within the food web as compared to other seasons. Bacterial networks associated with POM tended to be smaller than free-living networks in the Atlantic Ocean (Milici et al., 2016). Indeed common particle-associated bacteria (Flavobacteriales, Alteromonadales, and Rhodobacterales) made up a large component of our August networks (e.g., Buchan et al., 2014). Syndiniales, Dinophyceae, and ciliates were also abundant in the August network with a much more even distribution that in April, highlighting that both grazers and parasites were integral components of the August food web. The most connected protists belonged to the heterotrophic flagellate Palpitomonas (Yabuki et al., 2010) and the bacterivorous MAST-9, with a high percentage of significant negative correlations suggestive of an important predatory role.

CONCLUSIONS
We observed strong seasonal changes in the composition and connectivity of bacterial and protistan communities in nearshore waters of the eastern Alaskan Beaufort Sea (Figures 2, 3, 7). Environmental conditions beneath sea ice favored parasitism and chemoautotrophy, including the surprising finding of Zetaproteobacteria. The presence of an increased relative abundance of chemoautotrophs suggests that iron, methane, nitrogen, and sulfur cycling are important under the ice during a time when the food web is often considered to be less productive. In the spring, we observed the formation of a complex and highly connected, brackish microbial community highlighting the importance and influence of terrestrial inputs into coastal marine ecosystems. Given the freshening of the Arctic Ocean, these microbes may become increasingly important in Arctic marine food webs in years to come. Nutrient depletion over the course of the summer favored a shift from a diatomdominated food web to one characterized by an increased relative abundance of heterotrophic and mixotrophic protists, especially dinoflagellates, as well as picophototrophs Micromonas and Bathycoccus and other small phytoflagellates. Bacterial communities became increasingly enriched in common marine oligotrophic clades typically considered to have lower carbon demands and an increased ability to consume more recalcitrant organic matter. This shift to a more detrital food web in the late summer yielded a smaller and less connected network with longer paths between organisms than in April or June.
The Arctic is currently experiencing a number of physical changes that can have far-reaching effects on Arctic Ocean food webs. Surface air temperatures are warming at twice the rate of the rest of the globe, sea ice age and thickness continues to decline, and summer sea surface temperatures continue to show a warming trend year after year (Osborne et al., 2018). These changes are no doubt amplified in shallow, coastal systems such as our study system. These warmer temperatures also result in changes in precipitation and runoff patterns. We now have a baseline understanding of microbial communities in this region from which to predict community responses to a changing Arctic Ocean; one characterized by unique brackish communities with diatom blooms in the spring followed by long periods of nutrientpoor conditions in shallow waters inhabited by small grazers, picophototrophs, and oligotrophic bacterial clades. Continued long-term observations in this region are necessary to validate these predictions and assess their effects on higher trophic levels.

AUTHOR CONTRIBUTIONS
CK collected and processed samples, performed sequence and data analysis, and led writing. BC contributed to writing and interpretation. BC, JM, and KD conceptualized the project. JM and KD lead field efforts and collected samples. All authors contributed to manuscript revision, read, and approved the submitted version.