Spatial and temporal variation of Antarctic microbial interactions: a study around the west Antarctic Peninsula

Abstract

Background

The west Antarctic Peninsula (WAP) is a region of rapid environmental changes, with regional differences in climate warming along the north–south axis of the peninsula. Along the WAP, Palmer corresponds to a warmer region with lesser sea ice extent in the north compared to Rothera ~ 400 km to the south. Comprehensive and comparative, year-round assessments of the WAP microbial community dynamics in coastal surface waters at these two locations are imperative to understand the effects of regional climate warming variations on microbial community dynamics, but this is still lacking.

Results

We report on the seasonal diversity, taxonomic overview, as well as predicted inter-and intra-domain causal effects (interactions) of the bacterial and microbial eukaryotic communities close to the Palmer station and at the Rothera time-series site between July 2013 and April 2014. Our 16S- and 18S-rRNA gene amplicon sequencing data showed that across all seasons, both bacteria and microbial eukaryotic communities were considerably different between the two sites which could be attributed to seawater temperature, and sea ice coverage in combination with sea ice type differences. Overall, in terms of biotic drivers, causal-effect modelling suggests that bacteria were stronger drivers of ecosystem dynamics at Palmer, while microbial eukaryotes played a stronger role at Rothera. The parasitic taxa Syndiniales persevered at both sites across the seasons, with Palmer and Rothera harbouring different key groups. Up to 62.3% of the negative causal effects were driven by Syndiniales at Rothera compared to only 13.5% at Palmer, suggesting that parasitism drives community dynamics at Rothera more strongly than at Palmer. Conversely, SAR11 Clade II, which was less abundant but persistent year-round at both sites, was the dominant driver at Palmer, evidenced by many (28.2% and 37.4% of positive and negative effects respectively) strong causal effects. Article note: Kindly check first page article notes are correct.

Conclusions

Our research has shed light on the dynamics of microbial community composition and correlative interactions at two sampling locations that represent different climate regimes along the WAP.

Highlights

• Causal-effect modelling indicates that bacteria are stronger drivers of ecosystem dynamics at PAL, while microbial eukaryotes play a stronger role at RATS;

• Parasitism drives community dynamics at RATS but not PAL;

• SAR11 driven causal effects are abundant at PAL, though not by the most dominant clade.

Background

The region around the west Antarctic Peninsula (WAP) has been touted as one of the marine ecosystems experiencing the highest warming rate globally (0.46 °C ± 0.15 °C per decade) [1]. Increased temperature, rapidly melting ice sheets, and changing hydrography ultimately increases the vulnerability of marine organisms in their natural ecosystem and impacts regional biogeochemistry [2,3,4,5,6]. The implications of these climate change effects are significant since the WAP region is ecologically important – housing a productive food web [7], high rates of primary productivity [8] and significant CO₂ flux between the air and sea [9, 10]. These properties also exhibit strong seasonal variability in the WAP region due to large differences in irradiance and sea ice cover within the year, and interannual variability in climatically modulated processes (e.g. Amundsen Sea Low, Southern Annular Mode (SAM) and/or the El Niño-Southern Oscillation (ENSO) controlling mixing, stratification, deep-water nutrient supplies, sea ice cover, sea ice melt timing etc.) [11,12,13]. Annual sea ice extent and duration, which can have significant cascading effects on ocean heat flux, irradiance and consequently the marine food web [14], is declining following a north-to-south gradient at the WAP. The north experiences higher temperatures, lower ice extent and shorter duration of sea ice cover than the south [15], adding additional layers of complexity to predicting climate change effects on the WAP ecosystem [16].

The microbial community, consisting of bacteria, archaea, microbial eukaryotes (including but not limited to phytoplankton, heterotrophic protists, fungi) and viruses are key players regulating the productive WAP food web and biogeochemical cycling [17,18,19,20,21]. The WAP microbial community composition is highly dynamic, influenced by the physical and environmental forces that vary between seasons, years and decades due to climate change [22]. Furthermore, microbes form complex ecological interactions between each other, both within and between domains and with their environment [23]. These interactions can have either positive or negative effects on the species involved within the community, ultimately influencing the ecosystem’s microbial community composition in addition to influences by physical and environmental forcings, and seasonal patterns [24, 25]. Network inference from observational microbial abundance data can generate hypotheses about these interactions, and has delivered relevant ecological insights into microbial community structure of marine habitats elsewhere [26,27,28]. In this study, we used causal effect analysis [29,30,31], a network reconstruction approach originally developed to predict the impact of gene knock-out experiments [29]. The method was successfully used to predict interactions of bacteria and eukaryotes on marine plastic debris earlier [26]. Here, we integrated bacterial and micro-eukaryotic ASV relative abundances with measurements of environmental parameters and nutrients, thus allowing us to assess both biotic interactions and effects of physicochemical parameters on microorganisms simultaneously. Causal effect analysis has been shown to outperform other approaches for predicting cause-effect relationships including correlation-based approaches [30], by exploiting the local structure of the reconstructed graph to predict an interaction [29].

At the WAP, we now have a reasonably good grasp of the spatial and temporal (seasonal and interannual) patterns of physical climate dynamics, fluctuations in primary productivity rates and nutrient levels [31, 32], and some microbial community dynamics. This is in large part due to observations generated by long-term monitoring programs such as the Rothera Time Series program [33] and Palmer Long-Term Ecological Research program [34]. However, studies exploring the temporal/seasonal microbial composition and interaction patterns within the austral fall and winter seasons are still lacking, or do not consider seasonal successions continuously throughout the year. While productivity is low during the winter, understanding the composition and dynamics of the overwintering microbes and their alternative nutrient cycling pathways are critical since they influence the community and activity of the following seasons [35,36,37,38,39]. Additionally, other than the study by Luria et al. [40], most studies did not concurrently consider microbial community composition across domains (i.e. prokaryotes and eukaryotes from the same samples); thus inter-domain interactions have not been extensively investigated. Luria et al. [41] reported on the seasonal succession of bacteria in the free-living fraction at Palmer, on which this study builds. An adequate understanding of the spatial and temporal patterns of the microbial community composition at the WAP, and the microbial interactions occurring within the communities are critical, since they can directly impact the efficiency of carbon fixation, nutrient cycling, export and overall balance of ocean productivity.

In this study, we fill these knowledge gaps by tracking the seasonal community succession and interactions of the marine WAP bacteria and microbial eukaryotes using two sets of year-long datasets. This study leverages on coordinated (time-points, methods) sampling from the 2013 to 2014 winter to summer season at both Palmer and Rothera research stations, which allowed us to jointly consider the effects of localized environmental conditions and climatic variations (here, temperature and sea ice duration); on microbial community dynamics. While the WAP is warming at a faster rate than the global average [1], there is local variation: The northern part around Palmer station has warmed more and experiences higher reduction in sea ice cover than southern areas around Rothera [33]. Integrating microbial community data from both of these sites, we expect microbial community differences and differences in microbial interactions related to temperature. We use 16S (bacteria) and 18S (eukaryotes) rRNA gene metabarcoding data processed as amplicon sequence variants (ASVs), generated from 198 surface water samples collected from the two sites, and measurements of environmental parameters (temperature, salinity, chlorophyll a, nitrate + nitrite [NOx], phosphate, silicate, sea ice coverage and sea ice type) from the same samples. From these datasets, we (1) provide new insights on the temporal shifts of the bacterial and microbial eukaryotic community, and (2) evaluate how the climatic differences between Palmer and Rothera relate to the observed seasonal patterns. We then focus on (3) the intra- and inter domain microbial interactions using causal-effect modelling of the microbial abundance data, where we explore the seasonal and climatic gradient trends of the predicted key interactions.

Methods

Study sites, contextual data, sample collection and processing

The datasets used in this study were based on coastal seawater samples collected from two stations with long-term monitoring sites at the WAP (Fig. 1). The first is Rothera Research Station on Adelaide Island, where samples were taken from 15 m depth at the Rothera Time Series site (67.04°S, 68.08°W; water column depth ~ 500 m) in Ryder Bay (referred to as RATS from hereon in). The second is the Palmer Station on Anvers Island (64.77°S, 64.05°W), where samples were taken from the seawater intake, sampling from 6 m depth (referred to as PAL from hereon in).

Fig. 1

Map showing Palmer and Rothera station and the location of the west Antarctic peninsula relative to the Antarctic continent (inset). The Rothera and Palmer stations are marked to indicate the relative positions of and distances between the Rothera (RATS) and Palmer sampling sites. Maps were generated using ggOceanMaps (v 2.2.0) [61] with glacier polygons of glaciated areas and Antarctic ice shelves from naturalearthdata.com

For both sites, the analyses considered samples collected from all four seasons, between July 2013 (austral mid-winter) to April 2014 (early autumn). One to three samples were collected weekly per site (with lower intervals in Autumn/Winter – down to the lowest of one sample in three weeks), resulting in a total of 32 RATS and 75 PAL samples. Each sample was sequentially subjected to 3.0 µm followed by 0.2 µm polycarbonate membrane filtration, and filters from both size fractions were retained and analysed (referred to hereafter as RF3.0 and RF0.2). Detailed sample collection, environmental parameter measurements, DNA extraction, 16S and 18S rRNA gene library generation procedures are described in Luria et al. [41] and Rozema et al. [42]. Briefly, the V6 bacterial hypervariable region of the 16S rRNA gene was targeted using the 967F and 1064R primer pair [43]. For the 18S rRNA gene library generation, the V4 hypervariable region of the gene was targeted by PCR amplification using the V4F and V4RB primer pair and amplification conditions as described in Balzano et al. [44]. Positive and negative controls were included in all PCRs. Paired-end sequencing of the 16S and 18S rRNA gene libraries were performed on the HiSeq and MiSeq platforms, respectively as descried in Luria et al. [45] and Rozema et al. [42]. RATS sea-ice data was obtained from Venables et al. [33], while PAL sea-ice data was obtained from Environmental Data Initiative Portal [46] and Goodell et al. [47] (see Additional File 5).

Imputation of missing contextual data

For several of the PAL samples, some contextual data were not available due to equipment malfunction. Briefly, missing data for salinity, NOx, phosphate and silicate were imputed using the R package ‘mice_3.16’ [48] with the weighted predictive mean matching method (m = 10 imputations, maxit = 50 iterations per imputation). Dates and percentage of dates with missing data are detailed further in Supplementary Methods (Additional File 1). Imputed values are indicated in Table S1 (Additional File 2). Some parameters were transformed before causal effect analysis as described in Supplementary Methods (Additional File 1).

Sequence analyses

All raw, demultiplexed DNA sequences were downloaded from the Visualization and Analysis of Microbial Population Structures (VAMPS) portal (https://vamps2.mbl.edu/) from projects ‘LAZ_PAL_Bv6’, ‘LAZ_PAL_Ev4’, ‘LAZ_RAT_Bv6’ and ‘LAZ_RAT_Ev4’. Sequencing adapters and primers were removed using Cutadapt [v3.1; [49]]. The data were then processed using the CASCABEL amplicon sequence analyses pipeline [v.4.4; [50]] based on the standard configuration file for generating amplicon sequence variants (ASV) with paired-end reads, but with modifications to the following parameters applied to the 16S rRNA gene library: (1) at the ‘extract barcodes’ section, bc_length was set to --bc1_len 9 and --bc2_len 9"; (2) at the ‘dada2 trim and filter reads’ section, maxEE_FW was set to 2 and maxEE_Rv was set to 4; (3) at the ‘dada2 merge pairs’ section, minOverlap was set to 40 and maxMismatch was set to 10 to prioritize large overlaps of the forward and reverse read as expected from the amplicon length. At the ‘Assign taxonomy’ section, Silva (nr99, v138.1; Quast et al. [51]) and PR2 (v5.0; Guillou et al. [52]) databases were used for the 16S and 18S rRNA gene libraries respectively, with minBoot set to 60.

Sequences (ASVs) that classified as plastid (chloroplast or mitochondrial) sequences from taxonomic assignments with the Silva database were removed. Plastid sequences were further filtered out from the 16S rRNA gene libraries by taxonomy assignment with Global Alignment for Sequence Taxonomy (GAST; Huse et al. [53]) and PhytoREF [54]. After eliminating samples with insufficient sequencing depth (less than 40,000 and 10,000 reads for bacterial and eukaryotic amplicons, respectively), 30 RF0.2 and 28 RF3.0 RATS samples, as well as 69 RF0.2 and 72 RF3.0 PAL samples yielded both bacterial and eukaryotic sequences. The full list of samples for each size fraction and amplicon, and their corresponding contextual data are in Table S1 (Additional File 2).

Causal effect analysis

We retained ASVs that occurred at least 3 times in a minimum of 30% of the samples for the causal effect analysis, which resulted in 896 bacterial and 499 eukaryotic ASVs (RATS RF0.2), 616 bacterial and 325 eukaryotic ASVs (RATS RF3.0), 2014 bacterial and 383 eukaryotic ASVs (PAL RF0.2), as well as 2381 bacterial and 373 eukaryotic ASVs (PAL RF3.0) remaining for the subsequent analyses. For each sampling site and size fraction combination, the bacterial and eukaryotic ASV tables were merged and transformed by variance stabilizing transformation to render approximately Gaussian distributed ASV abundances. These were then merged with their respective contextual data (some transformed, see Supplementary Materials in Additional File 1 for details). We then generated 100 subsamples of the data, where ¼ of the samples (site/size fraction combination) were left out for each subsample run, standardized the features, and then ran the aIDA algorithm [55] with alpha = 0.5. For each subsampling run, aIDA reconstructs a partially directed acyclic graph from the feature abundances using the pc-algorithm which makes use of partial correlations and conditional independence tests [56]. Assuming sparsity, only the strongest correlations, explaining unique parts of the variation between features are represented by an edge in the graph. Subsequently, conditional independence tests allow to infer directionality for some but not all edges in the graph, since multiple graph skeletons can encode the same conditional independence assumptions. The local structure of the graph is then used to estimate causal effects between ASVs and environmental parameters and ASVs. This is done by means of linear regression with the cause feature and its direct upstream nodes (‘parent’ nodes) as explanatory variables and the target feature as dependent variable [29]. The regression coefficient of the cause represents the causal effect and is ‘corrected’ for spurious effects that should be attributed to nodes upstream in the graph. A causal effect of e.g. 1.5 indicates that if the causal feature increases by one standard unit of relative abundance, the affected feature will increase in 1.5 standard units of relative abundance. Subsequently, results from the subsampling runs were aggregated for the final causal effect estimate, as described in Taruttis et al. [55].

Statistical analysis

Multivariate analyses were run in R (v4.2.2) [57] with the packages ‘vegan_2.6-4’ [58], ‘ecole_0.9-2021’ [59] and ‘microViz_0.10.10’ [60]. ASV count tables were rarefied (Rarefy – ‘GUniFrac’) using the lowest total number of counts per sample per dataset as subsample before alpha diversity analyses (richness; Shannon index, inverse Simpson index; evenness)(hill_taxa—‘hillR’). Rarefaction curves and saturation points were drawn with the rarecurve function (‘vegan’). PAL and RATS microbial community ASV tables from both the RF0.2 and RF3.0 fractions were square-root transformed and Wisconsin double standardized (sqrt, wisconsin, dist_calc) to allow sample comparisons. From these data, we generated Bray–Curtis dissimilarity matrices. We then performed permutational MANOVA (PERMANOVA) with pairwise post-hoc tests with ‘PRIMER’ (v7, with PERMANOVA +) on Bray–Curtis dissimilarity matrices of the size fractions (RF0.2, RF3.0) individually using Site (PAL, RATS) and Season (spring, summer, autumn, winter) and their interaction term as explanatory variables. To test for homogeneity of dispersions among groups (as linked to their explanatory variables), we performed PERMDISP tests (‘PRIMER’). Similarities between community profiles were visualized by non-metric multidimensional scaling (nMDS) (metaMDS, ord_plot) on Bray–Curtis dissimilarity matrices.

Results and discussion

Seasonal water and sea-ice profile dynamics across PAL-RATS

Water temperature at RATS was lower on average than PAL by 0.87 °C, a pattern which persisted throughout all seasons (Fig. 2A) and reflects the south-north warming gradient at the WAP. As expected, chlorophyll a peaked in the summer, and its dynamics were very similar between both sites throughout the seasons with values peaking approximately 2 weeks later at PAL (Fig. 2C). Phosphate concentrations were similar at the two sites but displayed large fluctuation in the summer and autumn (Fig. 2D). In the winter and spring, NOx patterns were similar at both sites (27–35 µmol/L) and were consistently slightly lower (~ 5 µmol/L) at PAL (Fig. 2E). NOx levels dropped sharply in the summer displaying large variation at both sites, similar to trends also observed in 2015–2016 [62]. This resulted in comparable nitrogen to phosphorus ratios (N:P; Fig. 2G) between the two sites in winter, spring and autumn, following N:P observed over multiple years at RATS [62]. However, there were outlier values in the summer, where the ratio was very low (~ 1) for two time points at PAL (an indication of N limitation due to intense algal blooms) and high (> 20 for two time points; indicating P limitation) at RATS. These fluctuations in NOx, phosphate and their ratio might be due to phytoplankton blooms and viral lysis events, which may also be influenced by ocean temperature, climate and physical forcings (SAM/ENSO) [63,64,65].

Fig. 2

Seasonal physicochemical and nutrient profiles for the RATS and Palmer water samples from July 2013 to May 2014. For salinity, only RATS measurements were considered. See supplementary methods (Additional File 1) for more details on the physicochemical parameters not considered in the causal effects analysis of PAL samples

Maximum silicate levels were higher at RATS than PAL in winter and spring, while in the summer silicate decreased slightly at PAL and drastically at RATS from ~ 80 µmol/L to a low of ~ 15 µmol/L mid-summer (Fig. 2F). There were no signs of silicate limitation at both sites, concurring with the findings of previous studies [31, 65]. Overall, ice-retreat (timing, variability, sea-ice duration) and spring/summer ice types differed between RATS and PAL, which could influence phytoplankton and microbial community dynamics by altering light availability, water column stratification and nutrient availability. When calculated based on the definition of Stammerjohn et al. [13], day of retreat was 46 days later at RATS than PAL (Fig. 2H). Following sea-ice retreat at PAL (Dec 21), rapid retreat of mixed thin first-year and new/young ice (Fig. S1 in Additional File 1) left behind relatively high sea ice meltwater fractions in the surface waters [66]. 2013–2014 was a high ice year with retreat happening nearly a month later than the mean retreat day of Nov 29 (calculated based on 1992–2020 PAL in-situ records). The technical day of retreat at RATS was Feb 5, however, ice concentrations quickly returned to two-tenths or above for five consecutive days on March 14, 2014 at the start of the 2014 ice year. Brash ice persisted at RATS throughout the summer and fall at concentrations ranging from 10 – 80%.

Diversity and biogeography of PAL-RATS bacteria and eukaryotes

Following QC, the datasets yielded 136 million bacterial and 14.2 million eukaryotic sequencing reads, from which we identified 19,125 bacterial and 6,728 eukaryote ASVs, respectively. Rarefaction curves (Fig. S2 in Additional File 1) indicated both bacterial and eukaryotic ASV richness to be comprehensively represented at all sites and size fractions. Lower ASV richness (Fig. S3 in Additional File 1; Table S2 in Additional File 2) of both bacteria and eukaryotes was observed in the summer compared to other seasons at all site/size fractions, in agreement with general seasonal trends observed in Antarctic waters [41, 67, 68]. Bacterial ASV richness was considerably higher in PAL than RATS in both size fractions, with PAL RF3.0 having the most complex community. The difference in overall bacterial ASV richness was much larger between the two size fractions at PAL. The same cannot be said for the eukaryotes, where ASV richness was much more comparable across both site/size fractions. Further, decreased bacterial ASV richness at PAL RF3.0 persisted for longer periods (6 weeks) compared to the RF0.2 fraction (4 weeks), suggesting a prolonged influence of biologically significant events such as phytoplankton blooms on bacterial production and richness [41]. Bacterial ASV richness fluctuated more but less drastically at RATS in the summer, which could indicate bacterial production influenced by smaller and episodic (shorter duration) phytoplankton blooms. There was less seasonal variation observed for the evenness of the community species distribution. Overall, it was consistently lowest throughout the year in PAL RF0.2 and fluctuated substantially throughout the year for the rest of the site/size fractions for both bacteria and eukaryotes.

nMDS ordination showed that collectively, bacteria and eukaryotes at both PAL and RATS had seasonally differentiated communities across both size fractions (Fig. 3). The communities within a common season were also different between sites. PERMANOVA and PERMDISP analyses supported community differences in centroids between some seasons and differences between both centroids and dispersion for others for the RF0.2 size fraction. For the RF3.0 size fraction, season, site and interaction effects were significant both for centroid and dispersion differences, supporting that they differed significantly in at least one of them (Table S3 in Additional File 2). The RATS summer and autumn communities were much more scattered compared to PAL at both size fractions. When comparing between size fractions individually at the PAL and RATS sites, we observed there to be little to no community overlap of the sample clusters at both PAL and RATS (Fig. S4 in Additional File 1), suggesting communities to be distinct between the size fractions (at ASV level).

Fig. 3

Non-metric multidimensional scaling (nMDS) ordination plots of microbial communities (combined bacteria and eukaryotes) in the RF0.2 (A) and RF3.0 (B) size fractions from Palmer and Rothera. Samples collected in July/August are classified as winter samples, indicated by blue (B) circles and stars for samples from PAL and RATS, respectively, while Mar/Apr samples are classified as autumn samples and indicated by orange (O) circles/stars (PAL/RATS). Ellipses were drawn manually to group samples from the same or successive months

Community composition across seasons and sites

The bacterial and eukaryotic community taxonomic profiles of PAL and RATS of the two size fractions are summarized in Fig. 4. For both the RF0.2 and RF3.0, PAL and RATS showed distinct seasonal shifts in composition, though both sites and size fractions also had taxa that persisted across the seasons and were found across both sites. For example, in the RF0.2 fraction, substantial proportions of bacterium SAR11 (primarily Clade Ia and Clade II), Pseudomonadales, Cellvibionaceae and eukaryotes Phaeocystis antarctica (Prymnesiophyceae) and Dino Group I Clade 1 (Syndiniales) were found across all seasons at both sites (Fig. 4A, B). In the RF3.0 fraction, Dinophyceae (primarily Gyrodinium fusiform) and P. antarctica were observed year-round in PAL and RATS (Fig. 4D). A base level of Fungi was always detected across all seasons in both PAL and RATS in the RF0.2 fraction (0.02–15%), though the most abundant ASV (ev4.asv.23) could not be classified beyond the subdivision level (Fig. 4B, D), and the remaining ASVs classified as fungi were low abundance and not classifiable to fine taxonomic resolution.

Fig. 4

Taxonomic composition of Palmer and Rothera microbial communities. A RF0.2 fraction bacteria; B RF0.2 fraction microbial eukaryotes; C RF3.0 fraction bacteria; D RF3.0 microbial eukaryotes. ‘Other’ species (or their lowest classifiable taxonomic level thereof) include all taxa that were filtered based on min_sample_abundance = 0.1 in the tax_filter function of ‘microViz’. "Fall" and "Autumn" are used interchangeably throughout this manuscript

There were also taxa that were dominant exclusively in either one of the sites, and these could be divided into taxa that persisted year-round at that site, or only within specific seasons. Examples of site-exclusive taxa persisting year-round were the dinoflagellate G. fusiforme (RATS RF0.2 & RF3.0 all seasons, summer-autumn dominant), and single-celled algae Pedinellales (PAL RF0.2 & RF3.0 all seasons, dominant early summer) (Fig. 4B, D). Conversely, microbes dominant only within specific sites and season were bacterium Sulfitobacter sp. and stramenopile MAST-1A (PAL RF0.2; late summer-autumn and late spring-early summer), and Pseudoalteromonas piscicida (RATS RF3.0) (Fig. 4). Colwellia sp. and more specifically C. psychrerythraea was persistent year-round in PAL but only dominant winter-early spring in RATS RF3.0 (Fig. 4A, C). A similar pattern was observed for NS3a marine group Flavobacteria (peaking in RATS RF0.2 only in early summer) (Fig. 4A). In general, we found a similar seasonal succession of bacterial community composition in PAL RF0.2 as reported earlier by Luria et al. [41], though we did not detect Polaribacter sp. (a dominant genus reported in their study); this may be attributed to the variations in methods used for abundance table generation and taxonomy assignment. Further, the higher prevalence of Sulfitobacter sp. during summer in PAL could be linked to presence of specific taxa of diatoms in PAL not found in RATS (detailed below), which may support increased numbers of Sulfitobacter through mutualistic relationships such as those previously reported by Amin et al. [69].

The cryptophytes, consisting primarily of Geminigera cryophila were present at much higher relative abundances in PAL compared to RATS, particularly during the summer. In the RF0.2 and RF3.0 fraction, the maximum G. cryophila relative abundances were 62.1% and 26.2% respectively in PAL compared to 8.6% and 11.5% respectively in RATS (Fig. 4B, D; Fig. S6 in Additional File 1). Increased proportions of cryptophytes compared to other key phytoplankton groups (diatoms, haptophytes) typically found in coastal Antarctic surface waters have been observed in other studies [70, 71]. This could be related to the warmer seawater temperatures at PAL which are closer to the reported optimum growth temperature of cryptophytes [72]. Warmer conditions are also linked to increased sea-ice/glacial melt (hence decreased salinity), which was experienced in PAL that year (high sea-ice meltwater fraction). Decreased salinity has been previously reported to be associated with higher occurrences of cryptophytes [70]. Further, the rapid and drastic retreat of sea ice at PAL that year (i.e. drastic increase in irradiance to water column; Fig. 2H; Fig. S1 in Additional File 1) may favour cryptophytes, touted to have better adaptability to significant light-level fluctuations and fast photo-regulation response [73]. Another key phytoplankton group, Pelagomonadales, was interestingly much more prevalent in early spring at both sites (Fig. S6A in Additional File 1). Considering relative abundances of both size fractions, they were more abundant at PAL where maximum relative abundances in RF0.2 were ~ 40% over ~ 4 weeks compared to maxima of 30% over ~ 2 weeks at RATS. Pelagomonadales were also the most dominantly detected phytoplankton in the winter at PAL, observations which were consistent with the preferences of Pelagomonadales for higher temperatures (PAL) and lower irradiance levels (under sea ice or in a well-mixed water column; conditions more likely in winter/spring) [74].

Diatoms are one of the key phytoplankton classes with respect to carbon export in the WAP [75,76,77] and were supported by replete silicate levels at both sites. The taxa of diatoms detected at both PAL and RATS (Fig. S6 in Additional File 1) corroborated previous findings [71, 78, 79] in general. Though both sites harboured similar (relative) abundances of diatoms (Fig. S6B in Additional File 1), drastically higher silicate drawdown (Jan-2014; Fig. 2F) suggests higher diatom productivity in RATS compared to PAL. While relative abundances are mostly only indicative of diatom (eukaryotic) numbers, a plausible explanation for the larger silicate drawdown is that RATS harbours higher abundances of larger size-class diatoms (e.g. Corethron inerme, Proboscia alata; 100–200 µm size range, Fig. S6B in Additional File 1) throughout the productive season. In PAL, the extensive numbers of Fragilariopsis cylindrus, known to thrive in sea-ice [80], in early summer (end Nov-early Dec) corresponded with the high (95%) sea-ice coverage observed. The main initial ice retreat on Dec 22 at PAL could have seeded and led to the spike in numbers of Thalassiosira tumida, Stellarima microtrias and Porosira glacialis [81,82,83] which we found exclusively at PAL. Ice-type has previously been found to influence diatom community composition [71]; compared to the mostly ice-free conditions at PAL in the summer (Jan-May), the mobile brash ice fragments of fluctuating concentration at RATS (Fig. S1 in Additional File 1) may lead to more variable light and stratification conditions throughout the summer, ultimately resulting in the very different summer diatom communities observed at PAL and RATS.

The Syndiniales group is globally abundant and diverse, including in Antarctica [84,85,86], and an obligate parasitic clade [87] touted to substantially influence nutrient cycling and population dynamics [85, 88]. Syndiniales were overall more prevalent in the RF0.2 fraction, and a higher relative abundance was detected in RATS, constituting up to 40% of the RF0.2 eukaryotic community in the winter, as compared to PAL (max 16.5%; Fig. 4B, D; Fig S5A). Group-I-Clade-1 Syndiniales was persistent but was not the dominant Syndiniales group year-round at both sites (Fig. S5 in Additional File 1). In the winter-early spring, there was a higher diversity of different Syndiniales families detected, some of which were part of strong causal-effect relationships (detailed in the interactions section). The seasonal transitions were similar at both sites (peak abundance winter-early Spring), but dominant groups (other than Group-I-Clade-1) varied between sites (Group-II-Clade 10 and 11 at RATS; Group-II-Clade-1 and Clade-6 at PAL) (Fig. S5A in Additional File 1). The higher prevalence of Syndiniales in the winter corroborates the findings of Cleary and Durbin [85], a surprising finding because in other environments, parasite abundances commonly decline in the winter and colder waters [89, 90] due to infection suppression at low temperatures [89], lack of nutrients, sunlight and prey. This may also suggest the presence of polar ecotypes of Syndiniales with varying thermal preferences and lifestyle/survival strategies.

Microbial interactions based on causal effect estimation

A key interest of this study was to examine the inter- and intra-domain microbial interactions, as well as the interactions of the microbial community with pivotal environmental parameters. Understanding these interactions, their variations between sites that represent a climatic gradient and how they vary across different seasons can shed light on the microbial ecosystem structure and dynamics with better depth and accuracy [24]. We used a causal effect estimation method [29, 55], which differs from traditional correlation based network reconstruction approaches in that it allows for the inference of directed effects (i.e. there is a cause and effect node/ASV); hence, a change in relative abundance of the causal node is interpreted to cause a change in relative abundance of the effect node [29, 91]. The size of the causal effect indicates the induced change in standard units of the relative abundance of the effect node, thus larger effect sizes indicate a more prominent predicted effect. The nature of the interaction, be it a direct physical interaction or for example a metabolic dependency, cannot be inferred from relative abundances alone; though hypotheses about these interactions can be developed and further tested with other approaches.

Microbial-microbial effects: We first looked at the distribution of effect sizes at each site/size fraction (Fig. S7 in Additional File 1) and chose effect size cut-off points to focus on the strongest interactions. The effect size distribution and range differed by dataset, therefore we selected the tails individually (RATS RF0.2: interactions ≥ 0.95 and ≤  − 0.95; RATS RF3.0: interactions ≥ 1.00 and ≤  − 1.00; PAL RF0.2 & RF3.0: interactions ≥ 1.50 and ≤  − 1.50). The interactions described in the manuscript are based on effects that remained post-filtering at the respective cut off points with the total number of interactions per site/size fraction indicated in Fig. 5A. Overall, the maximum effect sizes observed for microbial-microbial effects were 7.92 and − 8.44. We summarized the top 20 positive and negative microbe-microbe interactions across both Palmer and Rothera at the RF0.2 and RF3.0 size fractions in Fig. 6 and Table S4 (Additional File 3). Strong inter- and intra-domain interactions were observed (Fig. 6), but the ratio of inter- and intra- domain interactions varied widely between the sites and size fractions (Fig. 5A). Bacteria→Bacteria (arrow indicates a causal effect from a microbial ASV on another microbial ASV) interactions were always the most prevalent at all sites and size fractions, except in the RATS RF3.0 fraction, where Eukaryote→Bacteria and Eukaryote→Eukaryote interactions dominated the positive and negative interactions, respectively. In all site/size fraction combinations other than RATS RF3.0, Eukaryote→Eukaryote interactions were less prevalent. This was most apparent in the PAL RF0.2 fraction where accumulatively, Eukaryote→Eukaryote interactions only made up 0.6% of both positive and negative interactions, suggesting that bacteria may play a larger role in structuring the PAL microbial ecosystem than at RATS. This is congruent with the findings of Sailley et al. [92] who found an increased importance of carbon flow through bacteria in the (warmer) northern WAP based on their inverse food-web model approach. The increased carbon flow could be linked to the sloppy feeding of krill on large phytoplankton (e.g. diatoms) which resulted in an additional source of DOC for bacterial production. They concurrently also observed decreased primary production (likely due to the decrease in large phytoplankton) and increased proportion of small phytoplankton in northern WAP waters, which is in-line with our observations at PAL (Fig. S6A in Additional File 1).

Fig. 5

A Percentage of different types of predicted positive (pos) and negative (neg) effects from the RF0.2 and RF3.0 size fractions at Rothera and Palmer. Effects considered here were those ≥ 0.95 and ≤  − 0.95 in RATS RF0.2; ≥ 1.00 and ≤  − 1.00 in RATS RF3.0; and ≥ 1.50 and ≤  − 1.50 in PAL RF0.2 & RF3.0. B Percentage of predicted positive (pos) and negative (neg) effects of different environmental parameters on microbial ASVs from the RF0.2 and RF3.0 size fractions at Rothera and Palmer. Effects ≥ 0.5 and ≤  − 0.5 were considered. For both plots, the total number of effects are indicated above the bar of each site-size-fraction combination

Fig. 6

Sankey diagrams of the 20 strongest positive and negative microbial-microbial interactions inferred from causal-effect (response) modelling at the different sampling sites and size fractions. Positive and negative interactions are represented by green and red flow edges, respectively. Microbes are grouped at class level by colours explained in respective figure legends below the positive and negative interaction diagram

At RATS RF0.2, we observed strong Eukaryote→Eukaryote interactions, while the total number of Eukaryote→Eukaryote interactions were small. For example, we observed large (positive) effects between the dinoflagellate (class Dinophyceae) Prorocentrum sp. (ev4.asv.817) and P. antarctica (haptophyte class Prymnesiophyceae) (ev4.asv.322 and 258) (Fig. 6, Table S4 in Additional File 3) (indicating an increase in Prorocentrum sp. abundance causes an increase in P. antarctica abundances). Prorocentrum was previously isolated from colonial P. antarctica and associated with a potential kleptoplastic (plastid sequestering)-mixotrophic feeding strategy of Prorocentrum sp. [93]. The relationship has implications on changing the proportions of carbon-exported versus cycled and remaining in the Antarctic microbial loop. Prorocentrum sp. as causal taxa was found in up to 6.5% of the interactions considered at RATS RF0.2 (Fig. 7), suggesting it has a relatively strong influence on the RATS microbial community dynamics. Relative abundance trends of ev4.asv.817 with ev4.asv.322 and 258 (Fig. S8A in Additional File 1) showed a seasonal association between Prorocentrum sp. and P. antarctica, the two species co-occurred only during the spring P. antarctica bloom.The decline in Prorocentrum sp. (ev4.asv.817) may be linked to another Bacteria→Eukaryote interaction with bv6.asv.1343, identified as a Pseudoalteronomas. sp. (Fig. S8B in Additional File 1), which is known to have algicidal effects on dinoflagellates including those from the Prorocentrum genus [94,95,96,97].

Fig. 7

Percentage of causal (A, B) and response (C, D) taxa in each of the size classes for the bacterial and eukaryotic community at Palmer and Rothera. The bar charts on the right of the heat map indicate the number of unique ASVs per species (A, B) or per class (C, D) detected within the communities. The top heatmap in A-B represent the classes of microbes which are unclassified at the species level. In the heatmaps, effects of all ASVs within the indicated taxon are considered

At site-size fraction combinations (RATS RF3.0 and 0.2) where causal nodes in the interactions were prevalently eukaryotes (i.e. Eukaryote→Eukaryote or Eukaryote→Bacteria interactions), causal eukaryotes in the strongest interactions consisted primarily of dinoflagellates Gyrodinium fusiforme, Syndiniales (Group I Clade 5 and Group II Clade 12), unicellular chlorophytes (Pterosperma sp.) and Haptophytes (Haptolina sp.) (Fig. 5, Table S4 in Additional File 3). As expected, strong predicted interactions involving Syndiniales were primarily negative given their parasitic nature; they were found to infect other dinoflagellates such as Heterocapsa sp. (ev4.asv.143, infected by ev4.asv.258, a Group II Clade 12 Syndiniales), a known Syndiniales host [98, 99]. This parasitic Syndiniales-Heterocapsa relationship may represent a control mechanism of dinoflagellate blooms [98] occurring in the late summer, where we observed elevated relative abundances of Heterocapsa nei/rotunda concurrent with a lower prevalence of Group II Syndiniales (Fig. 4B, D; Fig. S8C in Additional File 1). We also observed strong negative interactions between ev4.asv.258 Syndiniales with two bacterial ASVs (bv6.asv.1418 and bv6.asv.1038), though they were unfortunately unclassified beyond domain level. Within RATS RF0.2, Syndiniales cumulatively contributed as causal taxa in 62.3% of the negative effects/interactions, of which Group II Clade 10 and 11 Syndiniales alone was involved in 42.9% of the interactions (Fig. 7), highlighting the strong influence of Syndiniales (and by extension, parasitism) in driving the microbial dynamics at RATS. The same cannot be said for PAL, where Syndiniales was a causal taxa in only 4.7% of the negative effects, driven solely by Group II Clade 7 Syndiniales. Host specificity and infection dynamics of Syndiniales remain mostly obscure. Based on predicted interactions at RATS RF0.2 (where Group II Clade 10 and 11 Syndiniales was the primary causal taxa; Table S4-I, J in Additional File 3), 84% of the 99 interactions involved bacteria (primarily family Rhodobacteraceae, SAR11 and Sphingomonadaceae) while the remaining interactions involved stramenopiles (MAST-8D, Parmales and Pelagomonadales), dinoflagellates (Prorocentrum sp.) and haptophytes (Prymnesiophyceae Clade B3). These findings suggest Group II Clade 10/11 Syndiniales have more of a generalist/opportunistic infection strategy as described previously for micro- and macroeukaryotic hosts [88, 90, 100]. The predicted negative interactions between Syndiniales and bacteria likely represent indirect effects that might be related to environmental conditions favouring parasitic or free-living lifestyles alternately. Moreover, we cannot completely exclude that these interactions might be related to compositional artifacts that could not be completely removed by variance stabilizing transformation.

One bacterial taxon that contributed substantially as causal taxon at RATS (but not PAL) was Nitrospina sp., contributing 3.1% and 29.6% of positive interactions at RATS 0.2 and 3.0, respectively (Fig. 7). It was the causal taxon of the top 4 positive causal effects at RATS 0.2, with bv6.asv.80 predicted to interact with Shewanella benthica (Enterobacterales) and other Gammaproteobacteria, Chaetoceros neogracilis (Mediophyceae diatom), and Diaphanoeca sp. (Choanoflagellatea, a bacteriovore) (Fig. 6). Nitrite oxidizers Nitrospina spp. are a critical component supporting the marine nitrogen cycle, and their positive interactions with C. neogracilis likely reflects their role as nutrient source for phytoplankton [101, 102]. On the other hand, some species of Diaphanoeca have been reported to be bacteria filter feeders [103].

An intriguing finding from the PAL RF0.2 fraction was that a single ASV, bv6.asv.45, classified as SAR11 Clade II, was the causal ASV that both positively and negatively affected 40 different ASVs (24 bacterial, 16 eukaryotic) (Fig. 6, Table S4 in Additional File 3). Cumulatively, SAR11 Clade II (5 unique ASVs) was found to be the causal taxon in up to 28.2% and 37.4% of positive and negative effects respectively at this site and size fraction (Fig. 7), and > 90% of these interactions were involving bv6.asv.45. The bacterial response ASVs were mainly Alpha- and Gammaproteobacteria, Planctomycetes and Campylobacteria, while the eukaryotic response ASVs were mainly dinoflagellates, Syndiniales and marine stramenopiles, Telenomia (known heterotrophic predators feeding on bacteria, pico- and nano-phytoplankton) (Fig. 4, Fig. 6C–D). SAR11 is known to be one of the most abundant and ubiquitous free-living bacterioplankton, with very small and streamlined genomes [104]. The strongest predicted interactions were not driven by the more dominant Clade Ia SAR11 ASVs (Fig. 4A, C), which supports the concept that abundant taxa do not always drive significant interactions within an ecosystem [43, 105]. The copious and strong observed effects, in combination with the streamlined genome of SAR11 suggests cross-feeding/syntrophy to be an important feature in PAL RF0.2. A nutrient/substrate limitation (e.g. DOC) or environmental stressor may be triggering higher levels of cross-feeding to be predicted here and not at RATS. Other studies have also found SAR11 to have high co-occurrences [106, 107]; pointing to interdependent relationships with the other microbes within the PAL RF0.2 ecosystem to obtain substrates that they are unable to self-metabolize. One example is the positive relationship observed between bv6.asv.45 with Alteromonadaceae bacterium bv6.asv.167 (Paraglaciecola sp.) and bv6.asv.529 (P. chathamensis) (Table S4 in Additional File 3). The Alteromonadaceae family members are known polysaccharide degraders commonly associated with algal blooms [108, 109]. When we tracked the relative abundance trends of bv6.asv.45 and bv6.asv.167 (Fig S8D), we observed that bv6.asv.45 is prevalent throughout the year while bv6.asv.167 increased mid-spring and peaked mid-summer, declining again in fall. This indicates that this is a seasonally observed causal effect that is also co-dependent on activities in a higher trophic level (e.g. algal blooms) during the productive months. Given the year-round prevalence of SAR11, this likely explains its dependence on the many other interactions/exudates from other microbes for survival.

The Sva0996 marine bacterium group, an Actinobacteria (Acidimicrobiia), was observed to be a dominant causal taxa in 24.3% of the negative interactions at PAL RF3.0 (Fig. 7). Specifically, the strongest predicted interactions were driven by Actinobacteria bv6.asv.172 with Micavibionales, Amylibacter, Cryomorphaceae, other Gammaproteobacteria and Hypotrichia (a ciliated protist) (Fig. 6, Table S4 in Additional File 3). For the eukaryotes, a causal taxon that was found to be dominant in both PAL fractions is Pelagomonadales clade B1, a heterokont picoeukaryotic alga, contributing both positively and negatively to between 26.6 and 65.9% of the interactions. Ev4.asv.53, the Pelagomonadales member that was most involved in strong causal effects, was frequently associated with bacteria such as Fulvivirga sp., Granulosicoccus sp. and several bacterial ASVs not classified beyond the domain level. While Pelagomonas was also prevalent at RATS (Fig. 3B), it was involved in substantially fewer (and weaker) interactions. Fulvivirga sp. is known to secrete carbohydrate (polysaccharide) degrading enzymes [110], and its strong positive relationship with Pelagomonadales indicates a role in phytoplankton degradation, consumption of phytoplankton derived organic matter and driver of organic matter remineralisation within PAL. The specifics of the negative relationship between Granulosicoccus sp. and Pelagomonadales is less clear, though Pelagomonadales have been postulated to potentially have a mixotrophic lifestyle involving bacterivory [111, 112].

Environmental parameter-microbial effects: the distribution of effect sizes involving environmental parameters were separately evaluated at each site/size fraction, and cut-off points of ≥ 0.5 for positive effects and ≤  − 0.5 for negative effects were applied to filter environmental-microbial effects (Fig. S9 in Additional File 1). The total number of interactions considered post-filtering are indicated in Fig. 5B. Overall, the maximum effect sizes observed for environmental-microbial effects were 1.22 and − 1.37. Temperature had some of the strongest causal effect relationships with microbial community members (Fig. 8) across all sites, particularly at PAL (both size fractions) where they made up more than 43% of the observed effects (Fig. S10 in Additional File 1). This supports the theory that the warming gradient between the two sites may play an important role in structuring the microbial community differences observed between the two sites. Gammaproteobacteria, Bacteroidia, Nitrospina, and Dinoflagellates (Dinophyceae) were among the taxa most likely (i.e. strong causal effects observed) and frequently influenced by temperature. For example, 40.7% of the temperature-linked positive effects at RATS RF3.0 involved Bacteroidia. Temperature was also predicted to have causal effects on Alphaproteobacteria (up to 20.1%, positive and negative), Syndiniales (up to 49.5%, primarily negative) and Filosa-Thecofilosea (40.9%, primarily positive), and Cryptophyta (17.9%, positive). The high percentages of positive temperature-cryptophyte causal effects supports the higher occurrence of cryptophytes (G. cryophila) at PAL (Fig. S6A in Additional File 1) which is also backed by strong positive correlation between temperature-cryptophytes observed in other studies [71].

Fig. 8

Sankey diagrams of the 20 strongest positive and negative environmental variable → microbial interactions at the different sampling sites and size fractions. Positive (A–D) and negative (E–H) interactions are represented by green and red flow edges, respectively. Microbes are grouped at the class level by colours explained in respective figure legends below the positive and negative interaction diagrams. Chla - Chlorophyll a, NOx - Nitrate + Nitrite, PO₄ - Phosphate, Sal - Salinity, SiO - Silicate, Temp - Temperature

We observed the strong negative effects at PAL RF0.2 fraction to be primarily phosphate driven; the 20 strongest predicted interactions were with bacteria that were from the classes Campylobacteria (Arcobacteraceae), Bacteroidia (Vicingus sp., Reichenbachiella sp., Flavobacteroiales), Gammaproteobacteria (Enterobacterales, Thalassolituus oleivorans and other Pseudomonadales) and Verrucomicrobiae (Luteolibacter sp.). A majority (40.7%) of phosphate driven effects were involving Gammaproteobacteria (Fig. 8, Fig. S10 in Additional File 1), followed by Bacteroidetes (12%) and Campylobacteria (7.2%). Chlorophyll a also contributed considerably to the strongest environmental effects at all sites, except at PAL RF3.0 where it made up less than 10% of the interactions. At all RATS size fractions, Chlorophyll a had strong positive influence on the Bacteroidia class (Fig. 8); other taxa that were similarly influenced include Paraglaciecola sp. (Gammaproteobacteria) and Group I Syndiniales (Table S5 in Additional File 4). Nitrogen to phosphorus ratio was not a dominant driver of predicted effects (interactions) at all sites and was only found to have a strong positive effect on MAST-2D Stramenopiles (Gyrista) at PAL RF0.2.

Salinity and silicate drove many of the 20 strongest interactions at RATS (Fig. 8, Table S5 in Additional File 4), the main taxa influenced being Gammaproteobacteria, Bacteroidia, Alphaproteobacteria, Syndiniales, Dinophyceae and Bacillariophyceae. In the RF3.0 fraction, Nitrospina, Verrucomicrobiae and Actinobacteria also seem to be negatively affected by both salinity and silicate (Fig. 8C-D, I-J). The negative relationship observed between silicate and Syndiniales in up to 49.1% of the interactions at RATS RF0.2 fraction (Fig. S10J in Additional File 1) suggests that Syndiniales may (indirectly) be dependent on diatoms for their parasitic diet, though silicate has also been previously found to be an important driver of Syndiniales composition in previous correlation network analyses [90].

Concluding remarks

We report on the first comparative year-long seasonal dynamics of bacteria and microbial eukaryotes from two WAP sites with a spatial climate warming gradient. Our datasets reinforce earlier observations that while the two sites RATS and PAL are within close spatial distance, they have considerably varied environmental conditions due to the climatic gradient. This effect cascades down to the bacterial and microbial eukaryotic communities, where we observed variations in seasonal patterns of (alpha) diversity and community composition between the two sites. Causal effect estimation showed diverging trends of the predicted microbe-microbe interactions between the two sites, with bacteria such as SAR11 being more frequently found to drive many strong causal effects at PAL, while eukaryotes played a larger role in structuring the RATS ecosystem, featuring a strong parasitic component. Future research avenues should consider focusing on multi-year comparisons of such seasonal dynamics to characterize interannual variability of these patterns and include investigations into other key microbial domains that are known to influence biogeochemical cycling such as archaea and viruses.