Water samples for DNA and flow cytometry analysis were collected within the NGA LTER site at various depths, including surface, ten meters, and at the depth of the deep chlorophyll maximum (DCM) (Cohen, 2022). Three transects: the Seward line, Middleton Island line, and Kodiak Island line were sampled to characterize the spatial variation of water properties throughout the NGA ( Figure 1 ). Sampling of the Seward and Middleton Island lines was conducted in summers 2018–2021 while the Kodiak Island line was sampled summers 2018–2019 and 2021 ( Supplementary Table 1 ). In this study, summer refers to sampling conducted between June 25 and July 19 of each year, consistent with the seasonal cruise schedule of the NGA LTER (Strom et al., 2006; Aguilar-Islas et al., 2016). The 2018 cruise was conducted onboard R/V Woldstad and the other cruises were onboard R/V Sikuliaq. Regions were a priori selected based on distance from shore and their proximity to currents and bathymetric features ( Figure 1 ). Nearshore stations were within 30 nautical miles or 55.6 km from shore in a region influenced by the ACC. This region is characterized by low-salinity, iron- and copper-rich waters driven by freshwater discharge and terrestrial inputs from systems like the Copper River plume (Stabeno et al., 2016; Ortega et al., 2025; Reister et al., 2024). Shelf stations greater than 30 nautical miles from shore and over the continental shelf as defined here as areas shallower than 300 meters depth. Offshore stations were defined as those over the continental slope (> 300 m bottom depth) and are influenced by the Alaska Current - Alaskan Stream and mesoscale eddies, which facilitate offshore-shelf exchange and contribute to nutrient patchiness (Okkonen et al., 2001; Strom et al., 2010; Coyle et al., 2019). These spatial divisions are consistent with hydrographic and nutrient patterns described in previous studies of the NGA, including variations in salinity structure, current systems, and nutrient distributions (Stabeno et al., 2016; Aguilar-Islas et al., 2016; Strom et al., 2006; Reister et al., 2024; Ortega et al., 2025). The map of sample locations was produced using Ocean Data View (Schlitzer, 2022).

Whole seawater was collected for flow cytometry in cryovials and preserved with a final concentration of 0.5% glutaraldehyde. Samples were fixed in the dark for ten minutes, flash-frozen in liquid nitrogen, and stored at -80°C until analysis. EasyCheck Beads were used to calibrate the flow cytometer, following a 20x dilution and running three replicates (Cytek, Fremont, USA). For DNA analysis, samples were stained with SYBR Green (100x dilution), incubated in the dark for ten minutes, and analyzed alongside unstained controls using a Guava 5ST easyCyte flow cytometer equipped with a 488-nm laser and detectors for side scatter, forward scatter, red, green, and yellow fluorescence (Luminex, Austin, USA). Following known optical properties, heterotrophic bacteria, and Synechococcus were manually gated using Guava InCyte (Luminex, Austin, USA).

Approximately 1–5 liters of seawater was filtered through a 0.2 µm pore Sterivex filter (Fisher Scientific, Hampton, USA) and frozen at -80°C until extraction. DNA was extracted from Sterivex filters using DNeasy PowerWater Sterivex kit (Qiagen, Valencia, CA, USA, Cat. No. 14600-50-NF) following the manufacturer’s protocol. Microbial community composition was analyzed by amplifying the 16S rRNA and 18S rRNA gene sequences using TaggiMatrix primer multiplexing (Glenn et al., 2019). Prokaryotic communities were amplified using the 16S (V4 region) 515F and 806R primers (Parada et al., 2016). Eukaryotic communities were amplified using the 18S (V9 region) 1389F and 1510R primers (Amaral-Zettler et al., 2009). PCR was performed using the HiFi HotStart PCR Kit (Kapa Biosystems, Wilmington, MA, USA) on a Veriti 96-Well Thermal Cycler (Applied Biosystems, Waltham, MA, USA). The thermocycling profile consisted of an initial denaturation at 98°C for 1 minute, followed by 25 cycles of 98°C for 20 seconds, 55°C for 20 seconds, and 72°C for 30 seconds, with a final extension at 72°C for 5 minutes. Amplicons were cleaned using a HighPrep PCR Clean-up System (MagBio, Gaithersburg, MD, USA), and a second PCR was performed to add Illumina adapters and dual indices. Samples were sequenced using the Illumina MiSeq Platform (Illumina, San Diego, USA). A total of 232 prokaryotic and 198 eukaryotic samples were sequenced (PRJNA887083). The difference in sample numbers between prokaryotic and eukaryotic datasets reflects sample loss due to sequencing failure or low read counts that did not meet quality thresholds after rarefaction.

Raw sequences from Illumina were demultiplexed using Mr. Demuxy and processed using QIIME2 (version 2020.8) (Cock et al., 2009; Bolyen et al., 2019). Sequences were quality filtered (Phred score of 20) and denoised using DADA2 (Callahan et al., 2016). Operational taxonomic units (OTUs) were de novo clustered to 97% nucleotide identity and assigned taxonomy with the feature-classifier (Bokulich et al., 2018) using the Silva version 138 and PR² reference databases (Guillou et al., 2012; Quast et al., 2013). OTUs were used instead of ASVs because preliminary analyses showed that minor nucleotide differences (particularly among eukaryotic sequences) led to over-splitting of biologically similar taxa, which were then excluded from downstream analysis. This seemed to disproportionately reduce eukaryotic representation in the networks and biased results toward prokaryotes. Clustering at 97% identity allowed for more balanced representation of microbial domains while maintaining ecological relevance for network construction. Non-target 16S and 18S sequences, such as those identified as metazoan, chloroplast, and mitochondrial sequences were removed before further analysis. Samples were rarefied to a depth of 5,000 reads for the 16S rRNA gene and 2,500 reads for the 18S rRNA gene. These rarefaction depths are consistent with other marine studies using amplicon community analysis for eukaryotes, which can yield fewer sequences per sample due to differences in primer specificity and template abundance (Logares et al., 2014; Massana et al., 2015; Needham and Fuhrman, 2016). Rarefaction curves confirmed that these thresholds captured the majority of observed diversity ( Supplementary Figure 3 ).

For spatial comparisons of microbial community composition and co-occurrence network construction, samples were grouped by region (nearshore, shelf, and offshore) ( Supplementary Table 1 ). Although individual samples were collected and sequenced at distinct depths, including surface, ten meters, and a deep chlorophyll maximum (DCM), samples from each depth were grouped by region during network construction to evaluate broader spatial patterns. This approach maintained depth-level resolution during sequencing while enabling regional-scale ecological analysis.

Microbial taxon-taxon co-occurrence networks were constructed using relative read abundance data derived from rarefied OTU tables. Relative abundance was used to normalize differences in sequencing depth and emphasize proportional relationships among taxa across samples. While relative abundance helped normalize sequencing depth, we recognize that rRNA gene copy number varies widely, especially among eukaryotes. We did not apply copy number correction because our goal was to compare co-occurrence patterns not absolute abundance and robust correction methods for mixed environmental datasets are still limited. This approach is consistent with other cross-domain network studies (Lima-Mendez et al., 2015; Faust et al., 2015; Tipton et al., 2018). Networks were constructed for samples with significant Spearman rank correlations (ρ > 0.8), with p-values adjusted for multiple testing using the Bonferroni correction method (adjusted p< 0.05) in R using the Hmisc package (Harrell and Dupont, 2019). Networks were visualized in Cytoscape (Shannon et al., 2003) and subnetworks and a merged network were constructed in Cytoscape using MetScape, a Cytoscape plugin (Basu et al., 2017). Nodes represent all OTUs with significant correlations, which are represented by edges. Network and node-based measurements were calculated using NetworkAnalyzer, a built-in Cytoscape plugin (Cytoscape version 3.8.2; Assenov et al., 2008). Degree centrality (connectivity) was calculated as the number of edges linked to each node (Diestel, 2005). Potential keystone microbes were identified as the top ten OTUs with the highest degree centrality values in each network. Degree centrality reflects the number of significant co-occurrence relationships a taxon has and is commonly used as a proxy for ecological influence in microbial networks (Berry and Widder, 2014; Banerjee et al., 2018; Herren and McMahon, 2017). While additional centrality metrics such as closeness were calculated, degree centrality was prioritized due to its interpretability, simplicity, and established use in identifying highly connected taxa that likely structure microbial communities. Closeness centrality of each node was calculated as the reciprocal of the average shortest path length (Newman, 2003), where the average shortest path length is the mean length of the shortest path between the node and all other nodes. Co-occurrence fragmentation (f) was calculated for each network as:

where f is the fragmentation index, CL is the number of disconnected subgraphs (components), and N is the total number of nodes in the network. This metric describes network stability, with values closer to 0 indicating a more cohesive network, and values closer to 1 indicating higher fragmentation and reduced connectivity (Berry and Widder, 2014; Yang et al., 2024; Pilosof et al., 2014). A full set of network statistic definitions can be found in Supplementary Table 2 .

Non-metric multidimensional scaling ordination (NMDS) and significance of the multiple regression using environmental fit (envfit) with Bonferroni correction were conducted to test significant differences in microbial communities using the R package vegan (version 2.6-4, R version 4.2.1) (Dixon, 2003). Multiple response permutation procedure (MRPP) was conducted to test significant differences between microbial community composition using vegan (version 2.4-2) (Dixon, 2003). Differences in network measurements between sampling regions were compared using the non-parametric Wilcoxon rank sum test with the ggpubr package (Kassambara, 2020).

Summer spatial variability in microbial communities and their interactions across the Northern Gulf of Alaska was investigated using seawater samples collected from three regions: nearshore, shelf, and offshore ( Figure 1 ). DNA was extracted and amplified for the 16S and 18S rRNA genes for taxonomic identification of microbial communities. After rarefying samples to 5,000 reads for 16S and 2,500 reads for 18S, 175 prokaryotic and 145 eukaryotic samples remained for downstream analysis. These depths were selected based on rarefaction curves ( Supplementary Figure 3 ) and are consistent with marine microbial diversity studies (Logares et al., 2014; Massana et al., 2015; Needham and Fuhrman, 2016). Summary statistics on sequencing output, OTU counts, and taxonomic composition are provided in Supplementary Table 3 . Barplots summarizing the relative abundance of major prokaryotic and eukaryotic taxa across regions are presented in Supplementary Figure S1 . Both prokaryotic (16S) and eukaryotic (18S) communities differed significantly by region ( Figure 2 , Bonferroni adjusted p< 0.05). The robustness of regional differences in community composition were further tested with multiple response permutation procedure (MRPP) and communities were found to differ significantly by region in both prokaryotes (p = 0.001) and eukaryotes (p = 0.001). Additionally, shifts in prokaryotic and eukaryotic communities were significantly correlated with environmental variables: temperature, salinity, nitrate, and silicic acid concentrations ( Figure 2 , Bonferroni adjusted, p< 0.05). MRPP analysis also showed significant differences by depth and year in both datasets (p ≤ 0.008). The strongest depth-related shifts occurred between surface and DCM samples for both prokaryotes and eukaryotes. Interannual differences were also detected, with particularly distinct clustering among eukaryotic communities across years. These patterns align with previous reports of temporal variability in the NGA (Cohen, 2022). Taxonomic shifts across depth and year were primarily driven by changes in the relative abundance of Alphaproteobacteria, Gammaproteobacteria, Syndiniales, and diatoms, while groups like Bdellovibrionota and Group I Syndiniales remained relatively stable across environmental gradients.

Shifts in the prokaryotic community were characterized by shifts in Alphaproteobacteria and Gammaproteobacteria. Alphaproteobacteria were the most abundant offshore (19.58%) and least abundant nearshore (14.45%) ( Supplementary Figure S1A ). Pelagibacter ubique (SAR11), a prominent Alphaproteobacteria, varied by region, with clades I and II highest nearshore, clade IV highest offshore, and lowest levels on the shelf ( Supplementary Figure S1C ). Gammaproteobacteria were most abundant nearshore (21.91%) and less abundant offshore and shelf. SAR86 within Gammaproteobacteria had the highest abundance offshore. Bdellovibrionota, also known as Bdellovibrio and like organisms (BALOs), showed relatively consistent abundance across nearshore and offshore regions (2.53%).

Shifts in summertime eukaryotic communities were characterized by changes in dinoflagellates, diatoms, and prymnesiophytes. Dinoflagellates (phylum Dinoflagellata) were most abundant on the shelf (42.59%), with slightly lower relative abundance offshore (41.73%) and nearshore (37.07%) ( Supplementary Figure S1B ). Within this group, the genus Dinophysis, a mixotrophic dinoflagellate, showed highest relative abundance offshore and lowest nearshore. The class Syndiniales, composed of parasitic dinoflagellates, was the most abundant across all regions. Group I Syndiniales were more common nearshore, while Groups II, III, and unclassified Syndiniales were more abundant offshore ( Supplementary Figure S1D ). Gyrista, a recently defined phylum that includes diatoms and other stramenopiles, was most abundant offshore (24.57%), with lower relative abundance nearshore (13.59%) and on the shelf (12.30%). Within this group, pennate diatoms (class Bacillariophyceae) were more abundant nearshore (5.79%) and were primarily composed of the genus Pseudo-nitzschia ( Supplementary Figure S1E ).

Differences in microbial interactions by region were investigated by constructing microbial co-occurrence networks for each region using OTUs as nodes connected by statistically significant Spearman rank correlations (edges) ( Figure 3 ; Table 1 ). The nearshore network contained the greatest number of nodes (n = 352), followed closely by the shelf network (n = 350), while the offshore network had the fewest (n = 226). In contrast, the shelf network had the highest number of edges (n = 1682), with the most positive correlations (n = 1552). The nearshore network had slightly fewer total edges (n = 1288), but still a high number of positive correlations (n = 1139). The offshore network had the fewest edges overall (n = 1127) and the highest proportion of negative correlations (n = 242). Eukaryotes dominated all three networks, representing the largest proportion of nodes in each region ( Figure 3 ).

Degree centrality was not significantly different between nearshore and offshore networks ( Figure 4A ). However, the shelf network had a significantly higher degree centrality than the other two networks. Closeness centrality was significantly different between all three networks ( Figure 4B ). The offshore network had significantly higher closeness centrality compared to nearshore and shelf. The average shortest path length of nodes was significantly higher in nearshore compared to both shelf and offshore ( Figure 4C ). Additionally, offshore and shelf had slightly lower fragmentation values (0.60, 0.59) compared to nearshore (0.63), although no statistical test was performed since fragmentation is a single summary metric per network ( Table 1 ).

Key microbes are those likely to play critical roles within a given network, such that their removal could disrupt the network by severing interactions among microbes. To identify these key microbes within each interaction network, we focused on the top ten microbes with the highest degree centrality, or the most connected nodes in each network ( Table 2 ). This cutoff was chosen as a consistent way to highlight the most highly connected and potentially influential taxa across networks, following previous microbial network studies (Berry and Widder, 2014; Banerjee et al., 2018). These microbes are likely to play critical roles within the networks, such that their removal could disrupt the network by severing the interactions between many microbes. Nearshore and shelf networks were the most similar, having five microbes in common. While nearshore and offshore only had three microbes in common. Only shelf and offshore networks had bacteria in the top ten highest degree centrality. Two of these bacteria were in common, both being within the Proteobacteria phyla including the Gammaproteobacteria SAR86 and the Alphaproteobacteria SAR11 clade IV. Two OTUs were in the top ten in all three networks. Both are within the eukaryote domain: with closest hits to the cryptophyte P. prolonga and the haptophyte Phaeocystis. These shared taxa are known to be ecologically important in marine systems and may contribute to different processes depending on local environmental conditions. For example, SAR11 and SAR86 are ubiquitous oligotrophic heterotrophic bacteria that may fulfill similar roles in utilizing dissolved organic matter (DOM), while cryptophytes and Phaeocystis are high-latitude primary producers associated with complex trophic interactions and bloom dynamics (Giovannoni, 2017; Schoemann et al., 2005).

Microbes with the top ten highest degree centrality in the nearshore network were all eukaryotes ( Table 2 ; Figure 3 ). Three microbes were only found in the top ten most connected taxa in the nearshore network: the haptophyte Chrysochromulina sp., diatom Arcocellulus cornucervis, and the green algae Micromonas commode. The shelf network, in contrast, had four bacteria in the top ten highest degree centrality OTUs. Two of the four bacteria were Alphaproteobacteria, SAR11 clade IV and SAR11 clade Ia, while the third was a Gammaproteobacteria, SAR86, and the fourth, NS5 marine group differed as a Flavobacteria within the Bacteridota phylum. The offshore network differed the most, having the least common top ten most connected microbes to both nearshore and shelf networks. Bacteria were highest in the offshore network making up half of the ten most connected. All five bacteria were within the Proteobacteria phyla: three Alphaproteobacteria (AEGEAN-169, SAR116, SAR11 clade IV) and two Gammaproteobacteria (SAR86). An unidentified diatom within the Bacillariophyceae class and a ciliate within the Oligotrichida order were also within the top ten most connected microbes in the offshore network.

Microbial community composition in both prokaryotes and eukaryotes differed by region in the NGA ( Figure 2 ). Previous research in the tropical Pacific Ocean and Sargasso Sea has observed similar distinct nearshore, shelf, and offshore prokaryotic communities (Wang et al., 2019; Tucker et al., 2021) driven by cross-shelf gradients in salinity and temperature. In the NGA, shifts in prokaryotic and eukaryotic communities were significantly correlated with shifts in salinity and temperature ( Figure 2 ). Prokaryotic communities exhibited interannual variability, consistent with previous observations of community shifts during the 2019 marine heatwave, which favored small celled phytoplankton (Cohen, 2022; Strom, 2023). Other studies in the NGA have reported regional differences in eukaryotic plankton communities, including protists such as dinoflagellates and microzooplankton like ciliates, often linked to cross-shelf salinity gradients and freshwater input from rivers (Coyle and Pinchuk, 2005; Beamer et al., 2016; Strom et al., 2024). We observed a similar trend in our amplicon data, with both prokaryotic and eukaryotic communities clustering by region ( Figure 2 ). Nitrate and phosphate concentrations were highest in offshore and shelf regions and lowest nearshore, indicating stronger macronutrient limitation in surface nearshore waters ( Supplementary Table 1 ). This pattern is consistent with prior observations (Strom et al., 2006; Aguilar-Islas et al., 2016; Stabeno et al., 2016; Coyle et al., 2019) and supports the idea that regional niche partitioning of microbial communities is influenced by water mass characteristics.

The ACC has a low-salinity core and transports micronutrients from terrestrial runoff, including inputs from the Copper River plume, to nearshore waters (Ortega et al., 2025). Although rich in micronutrients like iron and copper, the ACC is often limited by macronutrients such as nitrate during summer (Strom et al., 2006; Aguilar-Islas et al., 2016; Coyle et al., 2019). Past studies have observed chain-forming diatoms dominating these nearshore waters, particularly during summer stratification (Strom et al., 2006). Our amplicon data supported this pattern, with high relative read abundance of Bacillariophyceae diatoms, including Pseudo-nitzschia spp. and Arcocellulus spp., in nearshore samples ( Supplementary Figure S1E ). Pseudo-nitzschia is a chain-forming pennate diatom frequently associated with elevated nutrient concentrations and high productivity (Zhang et al., 2021; Moreno et al., 2022), traits characteristic of the ACC (Strom et al., 2006; Stabeno et al., 2016). Some Pseudo-nitzschia species are also known to produce domoic acid, though toxin levels in the NGA are generally low. In contrast, offshore waters showed high relative abundance of Actinocyclus, a large centric diatom. While often associated with nutrient-rich coastal regions, Actinocyclus can also appear in offshore waters when mesoscale eddies or cross-shelf exchange events introduce nutrient pulses to otherwise iron-limited surface waters (Stabeno et al., 2004; Crawford et al., 2007; Ladd et al., 2007). This regional distinction in diatom assemblages likely reflects niche partitioning shaped by nutrient availability along cross-shelf gradients and water masses.

To examine microbial interactions shared across environmental gradients, we constructed a merged network of significant correlations and nodes found in all three regions. Co-occurrence networks infer potential interactions based on patterns of co-presence or mutual exclusion, but these patterns can result from various ecological processes, including direct interactions, shared environmental preferences, or indirect associations (Fuhrman, 2009; Weiss et al., 2016). Positive correlations may suggest mutualism, shared environmental preferences, or metabolic exchange, but do not necessarily reflect direct interactions (Hibbing et al., 2010; Ghoul and Mitri, 2016). However, correlation alone does not imply causation, and method limitations such as the lack of temporal resolution and potential indirect effects must be considered when interpreting these networks (Needham et al., 2013; Kodera et al., 2022). Despite these constraints, co-occurrence networks are useful hypothesis-generating tools that can identify recurring ecological patterns and guide future experimental investigations.

The merged network revealed significant correlations and nodes shared across all three regions. Although there were significant differences between nearshore and offshore networks and different potential keystone microbes, some significant correlations were found in all three networks ( Figure 5 ). The cryptophyte P. prolonga was a potential keystone microbe in all three networks ( Figure 3 ; Table 2 ). P. prolonga is a common marine cryptophyte that is an important plastid donor in kleptoplastic interactions with ciliates and dinoflagellates (Altenburger et al., 2020; Cruz and Cartaxana, 2022). In our networks, P. prolonga was consistently associated with ciliates and dinoflagellates, including a grouping found in all three regions involving the ciliate order Oligotrichida and the dinoflagellate Gyrodinium ( Figures 5 , 6 ). Kleptopasty allows the host to obtain plastids from ingested photosynthetic prey (Tsuchiya et al., 2020; Cruz and Cartaxana, 2022). Marine oligotrich ciliates have been found to retain prey plastids which remain functional for days to weeks giving them the temporary ability to perform photosynthesis in addition to phagotrophy, making them non-constitutive mixotrophs (Stoecker et al., 1988; Jonsson, 1989; McManus et al., 2012). Ciliates within Oligotrichida have been identified in high latitude environments with the ability to perform mixotrophy (Stoecker and Lavrentyev, 2018; O’Hara, 2023; Strom et al., 2024), suggesting this strategy may be especially advantageous in cold, seasonally variable systems like the NGA. Mixotrophy enhances survival under low-prey or low-light conditions and supports vertical and horizontal transfer of organic matter (Stoecker et al., 2017). It can also lengthen food chains and improve trophic transfer efficiency by bridging microbial and metazoan consumers (Stoecker et al., 2017). Recent work highlights the ecological significance of mixotrophs in the NGA, where environmental gradients in light and nutrients may select for flexible nutritional modes (Strom et al., 2024). The consistent correlations between P. prolonga and Oligotrichida nodes suggest kleptoplasty is a prevalent nutritional strategy in the NGA, contributing to primary production and enhancing trophic connectivity in this high-latitude ecosystem.

Within the nearshore network microbes associated with parasitism and predation were more abundant and connected ( Figure 3 ). This includes potential keystone microbes and predatory BALOs ( Figure 7 ; Table 2 ). BALOs are parasitoid bacteria that feed upon and reproduce within gram-negative bacteria using bdelloplasts until the prey cell lyses releasing new BALO cells (Sockett, 2009; Negus et al., 2017; Laloux, 2019). BALOs have two life phases a motile attack phase where they can swim toward prey and a proliferative phase after attachment to prey (Sockett, 2009; Laloux, 2019). Subnetworks showed BALOs (including OM27) significantly correlated with flavobacteria such as NS7 marine group, Cryomorphaceae, and others in all three regions ( Figure 7 ). These flavobacteria are often particle-associated marine microbes (Mitulla et al., 2016; Milici et al., 2017). Since BALOs contain two different life phases, they may move toward particles where there is a higher abundance of potential prey (Lambert et al., 2006, 2011; Iida et al., 2009). Flow cytometry data showed higher heterotrophic bacteria and Synechococcus nearshore cell counts compared to offshore ( Supplementary Figure S2 ). Prey density has been linked to predator density as well as predator-prey interactions (Holling, 1959; Dick et al., 2014; Cuthbert et al., 2021). A minimum prey density is therefore required to support a predator population (Holling, 1959; Dick et al., 2014; Cuthbert et al., 2021). Higher cell counts could increase cell contact rates and therefore give parasitic and predatory microbes higher success in the nearshore network (Hu et al., 2013).

Interestingly, eukaryotic parasitoids, such as Cryothecomonas aesticalis, a nanoflagellate parasitizing diatoms (Drebes et al., 1996; Catlett et al., 2023), were the most connected microbes in the nearshore network ( Table 2 ), highlighting how complex microbial interactions like parasitism can dominate under dynamic, high-biomass conditions even when surface nutrients are seasonally depleted (Worden et al., 2015). Other highly connected microbes in the nearshore, such as Syndiniales Dino-Group I and MAST-1C, are also known or suspected parasites (Guillou et al., 2008; Massana et al., 2009; Lin et al., 2012; Cleary and Durbin, 2016; Clarke et al., 2019). These eukaryotic parasites were notably absent as the most connected microbes in the offshore network, where microbial parasites were not considered keystone species. The nearshore environment typically has higher micronutrient availability; however, in summer, surface macronutrient depletion may limit primary production, allowing excess micronutrients to persist. These conditions could still support abundant microbial communities particularly in subsurface layers or particle-rich zones that serve as hosts for parasites (Mitulla et al., 2016; Milici et al., 2017). Additionally, freshwater inputs nearshore often create a shallow pycnocline, which may concentrate potential hosts in the surface layer and enhance contact rates with parasites (Lucas et al., 1999). Greater environmental variability in salinity and temperature nearshore may also increase host susceptibility to infection, as physiological stress can compromise host defenses (Harvell et al., 2002). In addition, in macronutrient-depleted surface waters, some microbes may occupy deeper layers to access nutrients, where reduced light imposes further metabolic stress that could increase vulnerability to parasitism (Strom et al., 2006; Worden et al., 2015; Stoecker et al., 2017). In contrast, offshore waters are typically colder and iron-limited, with high nitrate but low chlorophyll concentrations characteristic of high-nutrient, low-chlorophyll (HNLC) systems (Aguilar-Islas et al., 2016; Coyle et al., 2019). Host density in these waters may be limited by iron availability, reducing opportunities for parasitic lifestyles (Boyd, 2007; Strom et al., 2024). However, episodic mixing events such as eddies or intrusions from the Alaska Current and Alaskan Stream may intermittently deliver iron to surface waters, potentially creating localized conditions favorable for parasitism (Ladd et al., 2007; Coyle et al., 2019). These ecological factors may help explain the stronger prominence of parasitic eukaryotes in the nearshore compared to offshore ecosystems.

Network analysis of microbial co-occurrences provides valuable insights into community structure, resilience, and stability by identifying keystone taxa, those microbes most crucial to maintaining ecological functions. Using measures like degree centrality, which reflects a taxon’s connectivity within the network, we can determine which species, if removed, would be most likely to disrupt community dynamics (Berry and Widder, 2014; Banerjee et al., 2018). The loss of these keystone microbes can lead to significant shifts in microbial community structure and function, as their absence can disproportionately affect ecological processes (Berry and Widder, 2014; Cagua et al., 2019). Network fragmentation, which occurs when key correlations between taxa are lost, further exacerbates this disruption, leading to communities that are less resilient to disturbances (Widder et al., 2014; Wang et al., 2023; Yang et al., 2024). Communities exhibiting higher functional redundancy and lower fragmentation are more likely to maintain stability and recover from perturbations, as they have a greater capacity to compensate for the loss of keystone taxa (Wang et al., 2023). Understanding these key interactions and network properties is essential when evaluating the stability and adaptability of microbial communities, especially in varying environmental conditions.

Potential keystone microbes were identified as the top ten taxa in each network with the highest degree centrality values ( Table 2 ; Figure 3 ). This cutoff follows previous microbial network studies and was selected to enable consistent cross-network comparison, though we acknowledge that degree centrality alone may not capture all aspects of keystone behavior. Our results suggest that potential keystone microbes shifted from predominantly eukaryotic taxa nearshore to an equal mix of bacteria and eukaryotes offshore ( Table 2 ; Figure 3 ). Offshore keystone microbes included many oligotrophic bacteria such as SAR11, SAR86, and SAR116 (Brown et al., 2012; Giovannoni, 2017; Hoarfrost et al., 2020; Roda-Garcia et al., 2021). While some of these taxa are also highly abundant, which may contribute to their high centrality, their consistent appearance as central nodes across samples suggests they may also play important ecological roles. Despite being abundant and widespread, many bacteria with oligotrophic lifestyles are adapted to low-nutrient conditions and remain poorly understood due to their uncultivated status (Faust and Raes, 2012; Zelezniak et al., 2015; Yan et al., 2023). These traits may confer an advantage in the iron-limited HNLC offshore waters of the NGA. Studies have shown that interactions among such streamlined taxa, including cross-feeding and metabolite exchange, can structure microbial communities and enhance ecosystem function, highlighting the importance of further investigating their ecological roles (Braakman et al., 2017).

The negative correlations between bacteria and eukaryotic keystone microbes found in our datasets ( Figure 3 ) are possibly driven by bacterivory, with eukaryotic microbes phagocytizing bacterial prey (Sakka et al., 2000; Christaki et al., 2002; Strom, 2002; Sarmento et al., 2013; Arandia-Gorostidi et al., 2022; Follett et al., 2022). For example, the cryptophyte Plagioselmis prolonga, a keystone microbe in the nearshore network, had exclusively negative correlations with bacteria ( Figure 5 ). Many cryptophytes are mixotrophs, supplementing autotrophy with phagotrophy of bacteria (O’Hara, 2023; Šimek et al., 2023), suggesting that P. prolonga’s role as a keystone microbe may be linked to its dual trophic strategy. Similarly, Phaeocystis, another keystone microbe present across all regions, had a high number of negative correlations with bacteria, possibly due to bacterial degradation of bloom remains, such as carbon released from cell lysis (Brussaard et al., 1995, 1996) or due to the ability of Phaeocystis to phagocytize bacteria (Koppelle et al., 2022). By linking microbial primary production and heterotrophic consumption, mixotrophs like P. prolonga and Phaeocystis may connect different parts of the microbial network by occupying different trophic modes. This metabolic flexibility likely contributes to their high centrality and supports their potential roles as keystone microbes in our dataset (Li et al., 2024). Similar patterns have been observed in other marine systems, where mixotrophy is associated with ecological versatility and central roles in microbial networks (Flynn et al., 2013; Mitra et al., 2014; Ward and Follows, 2016).

Network-wide stability metrics indicated that the offshore network exhibited nodes with significantly higher closeness centrality ( Figure 4B ) compared to the networks in the other two regions. Higher closeness centrality suggests that microbes in the offshore network are more interconnected, meaning disturbances may propagate more quickly, but the overall network is more resistant to fragmentation (Assenov et al., 2008; Berry and Widder, 2014; Cagua et al., 2019; MatChado et al., 2021). While we cannot directly infer specific interaction types from correlation-based networks, higher connectivity in other ecological systems has been linked to community stability and robustness (Jeong et al., 2000; Poisot et al., 2011; Hannigan et al., 2018; Qian and Akçay, 2020). Some studies have also found that highly connected networks are more common in systems where cooperative or mutually beneficial interactions prevail (Thébault and Fontaine, 2010; Montesinos-Navarro et al., 2017; Costas-Selas et al., 2024). Further investigation is needed to confirm whether similar dynamics are occurring in the NGA. In contrast, lower connectivity in the nearshore and shelf networks could indicate reduced stability, potentially driven by a greater prevalence of antagonistic interactions such as parasitism or competition. Although these interactions still contribute to network structure, they may result in more specialized or short-lived associations that fragment networks rather than broadly linking them (Weitz and Dushoff, 2008; Coyte et al., 2015). Antagonistic dynamics also tend to be more specific and variable than mutualistic ones, which could make strong correlations less consistent over time. This fits with ecological theory showing that mutualistic networks are often more cohesive and nested, while antagonistic ones can be more modular and fragmented (Thébault and Fontaine, 2010; Allesina and Tang, 2012). Alternatively, the lower connectivity we observed may reflect greater environmental variability in nearshore and shelf waters compared to offshore. These regions experience stronger fluctuations in temperature, salinity, and freshwater input due to river discharge, stratification, and glacial melt (Beamer et al., 2016; Ortega et al., 2025; Reister et al., 2024; Strom et al., 2024). Prior studies have linked that kind of variability to reduced stability and coherence in microbial networks (Shade et al., 2012; Tylianakis and Morris, 2017). Taken together, these patterns suggest that while the offshore network may be more resistant to perturbation due to its higher connectivity, the lower connectivity observed in nearshore and shelf communities could make them more vulnerable to environmental change or better able to reorganize in response to it through faster turnover or more flexible associations (Shade et al., 2012).

Defining a stable ecosystem, especially a marine one, is difficult since the system is inherently dynamic. Researchers have used fragmentation measurements of microbial co-occurrence networks as an indicator of ecological network sensitivity, with higher fragmentation suggesting greater instability (Widder et al., 2014; Banerjee et al., 2018; Yang et al., 2024). All three co-occurrence networks had fragmentation values greater than 0.5, suggesting microbial communities in all three regions of the NGA are more fragmented and potentially less stable than those observed in other marine ecosystems with lower fragmentation, such as freshwater lakes, the East China Sea, and seagrass meadows (Widder et al., 2014; Wang et al., 2023; Yang et al., 2024), though differences in taxonomic focus and how the networks were built make direct comparisons harder to interpret. Among the three regions, the nearshore network had the highest fragmentation value and the lowest connectivity ( Table 1 ), which could indicate greater vulnerability to disturbance. This region also experiences the highest environmental variability, such fragmentation may also reflect greater flexibility in microbial associations rather than instability per se. This region also contained a greater proportion of parasitic and predatory taxa, such as parasitic eukaryotes and BALOs, which have been shown to influence microbial community structure through top-down interactions and host-parasite dynamics (Thingstad and Lignell, 1997; Chow et al., 2013; Lima-Mendez et al., 2015; Worden et al., 2015). The presence of these taxa suggests strong top-down control, which may contribute to increased fragmentation. According to the kill-the-winner hypothesis, strain-specific predators like viruses can promote microbial diversity by targeting the most competitive taxa, preventing dominance and enabling coexistence (Thingstad, 2000; Winter et al., 2010; Thingstad et al., 2014). If parasitism disproportionately affects highly connected or central taxa, it may destabilize the network by weakening key interactions (Weitz and Dushoff, 2008; Chow et al., 2013; Lima-Mendez et al., 2015). Alternatively, top down pressure could allow less dominant taxa to emerge as new keystone species, potentially reshaping the structure and function of the microbial community. The high fragmentation in the nearshore network could reflect this where top-down control both supports diversity and increases the risk of instability.

Differences in network characteristics between nearshore and offshore regions may be driven by two key factors. First, microbial community composition and keystone taxa vary significantly between regions ( Figure 2 ; Table 2 ), influencing interaction types and network structure. Antagonistic interactions, such as parasitism or predation, may be more prevalent in the nearshore, where exploitative relations could dominate, whereas offshore microbial communities may experience a more complex mix of interactions. These interactions might not necessarily be mutualistic but could involve processes like micronutrient sharing or cross-feeding that may foster greater network cohesion, which we interpret here as one form of structural stability (Amin et al., 2012; Bertrand et al., 2015; Cooper and Smith, 2015; Arandia-Gorostidi et al., 2022). Second, geochemical differences between regions, driven by the mixing of water masses in the Northern Gulf of Alaska, impact nutrient availability and primary productivity (Strom et al., 2010; Aguilar-Islas et al., 2016; Coyle et al., 2019). Offshore microbial communities may be more stable, potentially due to selective pressures favoring micronutrient exchange processes such as cross-feeding or siderophore exchange, which could reinforce network stability (Amin et al., 2012; Morris et al., 2012; Zelezniak et al., 2015; D’Souza et al., 2018). At the same time, antagonistic interactions can also contribute to ecological stability by regulating population dynamics, suppressing dominance, and promoting turnover (Shade et al., 2012; Coyte et al., 2015; Tylianakis and Morris, 2017). However, the relative roles of different interaction types in promoting network stability remain debated in the ecological literature, and the link between interaction type and stability is likely context-dependent. In contrast, the nearshore environment, characterized by more variable micronutrient and freshwater inputs, may foster a greater prevalence of exploitative interactions, since environmental variability and resource pulses have been shown to push microbial networks toward more antagonistic or short-lived associations (Hawkes and Keitt, 2015; Mickalide and Kuehn, 2019). The differences suggest that moderate variability in nutrient availability and a mosaic environmental structure could promote microbial diversity and stability, but this remains speculative and contrasts with the observed fragmentation in the nearshore network. Ultimately, these findings raise questions about the potential for ecosystem transitions in the NGA. In 2019, a shift in microbial community composition was observed, including a notable increase in chlorophytes and oligotrophic bacteria (Cohen, 2022). The challenge moving forward is to identify early warning signs that might precede tipping points in ecosystem functioning. Identifying these signs could help predict future transitions and offer valuable insights for understanding the resilience of the system (Chisholm and Filotas, 2009; Scheffer et al., 2012).

This study found microbial communities, potential keystone microbes, and their correlations differ across the shelf/slope of the NGA. Potential keystone microbes that were the most connected in the micronutrient-rich nearshore network were all eukaryotes. However, as distance increased to offshore waters, there was a shift in keystone microbes from eukaryotes to smaller bacteria often found in oligotrophic waters. The offshore network exhibited higher connectivity and less fragmentation compared to the nearshore, suggesting it may be more stable either due to greater microbial interconnectivity or as a result of the less variable environmental conditions offshore. This stability may be further promoted by more cross-domain interactions and a more balanced representation between bacteria and eukaryotes as keystone microbes in the offshore network. While the nearshore network seemed to have a greater prevalence of possible antagonistic interactions, the offshore network likely involves a more complex mix of interactions, with processes like micronutrient sharing or cross-feeding potentially playing a role in fostering network cohesion. The difference in potential interactions suggests a shift from more exploitative relationships nearshore to a more cooperative system offshore where low micronutrient concentrations can be an environmental stressor. Additionally, this study identified two potential key microbial players found in all regions, the mixotrophic nanoflagellates Plagioselmis prolonga and Phaeocystis sp., suggesting these microbes strategies are well suited to environmental gradients across the shelf/slope domains, and may contribute disproportionately to microbial interaction network stability. Overall, different network characteristics and potential interactions were apparent by region, but further research is needed to determine how microbial community composition and environmental variability influence ecosystem stability. These findings contribute to the understanding of microbial ecology during summer conditions in the NGA, and could have important implications for understanding the resilience or vulnerability of the NGA in the face of climate change.