Six one-day cruises were performed onboard the R/V Mediterranean Explorer over a 2-year period: two in early winter (December 1st, 2014; November 17th, 2015), two in spring (March 24th, 2015; March 30th, 2016), and two in summer (July 14th, 2015; July 25th, 2016). For each cruise, four stations were sampled (labeled in ascending order with distance from shore) along a coastal-to-offshore transect from the Herzliya marina, Israel (Figure 1). At each station, we measured the vertical profile of in situ parameters up to a depth of 500 m using a conductivity, temperature, and depth (CTD) probe (SeaBird CTD profiler SBE 19plus V2), a fluorescence probe (Seapoint fluorometer) to measure chlorophyll a and turbidity with a Seapoint turbidity sensor. The fluorometer was calibrated with bottle chlorophyll measurements before each cruise. Both the fluorometer and the turbidity sensor were mounted on the CTD. Data were extracted using the Seasoft V2 software suite and plotted using Ocean Data View 4.7.4 (Schlitzer, 2015).¹

Seawater for nutrient, phytoplankton pigments, cell counts, DNA, and RNA analysis was collected from 10 m depth using a 12-bottle rosette sampler equipped with 8 L Niskin bottles. To better identify gradients along the transect, CTD and nutrient data were collected for two additional stations: one located between stations 1 and 2, and the other between stations 2 and 3. Station coordinates, bottom depths, and distance between stations and the shore are summarized in Table 1.

Satellite-based maps of surface chlorophyll concentration are derived from the Copernicus Marine Environment Monitoring Service (CMEMS).² We used the level-3 Mediterranean Sea reprocessed surface chlorophyll concentration product, consisting of merged SeaWiFS, MODIS-Aqua, MERIS, and VIIRS satellite data (OCEANCOLOUR_MED_CHL_L3_REP_OBSERVATIONS_009_073). Using multi-satellite data allows continuous tracking of fine-scale chlorophyll filaments as they are advected and deformed by the currents (Lehahn et al., 2017). Surface chlorophyll concentrations are estimated via the MedOC4 (Volpe et al., 2007) and the AD4 (D’Alimonte and Zibordi, 2003) algorithms for case-1 and case-2 waters, respectively. Spatial and temporal resolution is 1 km and 1 day, respectively. To reduce the area masked by clouds, chlorophyll maps from 30.11.2014, 20.11.2015, 01.04.2016, to 26.07.2016 were used instead of 01.12.2014, 17.11.2015, 30.03.2016, and 25.07.2016, respectively. Furthermore, maps for 30.11.2014, 24.3.2015, 14.7.2015, and 26.7.2016 are composed of 5, 3, 3, and 3 consecutive images, respectively.

Seawater for nutrient analysis was collected in 15 ml Falcon tubes pre-rinsed with sample seawater. For each nutrient, duplicate non-filtered samples were frozen onboard directly after collection and kept at −20°C until analysis. Silicate, nitrate + nitrite, and soluble reactive phosphate content were measured within 6 weeks calorimetrically on a LACHAT Instruments Quik-Chem 8,500 Flow Injection Analyzer (Grasshoff et al., 2009) at the service unit of the Interuniversity Institute for Marine Sciences in Eilat, Israel. The internal precision (the average deviation from the mean of duplicate samples) of the measurements and the effective limit of detection for these frozen samples were 0.04 μmol/L for phosphate and 0.05 μmol/L for both silicate and nitrate + nitrite, respectively. Accuracy was obtained by calibration against a suite of solutions prepared by dilution of commercial, high-concentration standards (1,000 mg/L, Merck). The calibration curve was run several times during every flow injection analyzer session.

During the July 2016 cruise, phosphate, nitrate, and in addition ammonia were also measured in non-frozen samples. These samples were filtered through 0.4 μm filters, transferred into 15 ml Falcon tubes, stored at 4°C, and analyzed within 24 h at the Morris Kahn Marine Laboratory, Sdot Yam, Israel, using a SEAL AA-3 autoanalyzer. Samples for silicate analysis from these samples were acidified and kept frozen until 5 days before analysis. Silicate, nitrate + nitrite, and soluble reactive phosphate were measured calorimetrically following the manufacturer’s operating protocols (SEAL version 2015). Silicate was measured by a modification of the silico-molybdate blue method (Grasshoff et al., 1983), nitrate + nitrite determined by a purple diazo dye after Cd reduction (Armstrong et al., 1967), while soluble reactive phosphate was determined using a long flow cell (50 cm) and a modified Murphy and Riley (1962) molybdophosphate blue method. Ammonia was determined using a fluorometric method based on Kérouel and Aminot (1997). The precision of these methods at the used range were 33, 21, 2, and 6 nM, respectively, and the limits of detection were 55, 63, 2, and 5 nM (defined as three times the precision of the seawater blank).

For flow cytometry sample collection, 1.5 ml triplicate seawater samples were fixed with glutaraldehyde (0.125% final concentration), incubated in the dark for 10 min, stored in liquid nitrogen onboard, and kept at −80°C in the lab until analysis. For analysis, samples were thawed at room temperature and run on a BD FACSCanto™ II Flow Cytometry Analyzer Systems (BD Biosciences) with 2 μm diameter fluorescent beads (Polysciences, Warminster, PA, United States) added as an internal standard. Samples were run first without staining to use the natural fluorescence of the cells (chlorophyll and phycoerythrin pigments) for identification of three cell types: (i) Prochlorococcus, (ii) picoeukaryotic phytoplankton, which were identified based on the PerCP channel and forward scatter with picoeukaroytes being larger and having higher fluorescence signals than Prochlorococcus; and (iii) Synechococcus cells, which were identified based on their emission in the phycoerythrin (PE) channel. Samples were then stained with SYBR Green I (Molecular Probes/Thermo Fisher) to a final concentration according to the manufacturer’s instructions and the total microbial population was counted. Data were acquired and processed with FlowJo software. Flow rates were determined several times during each running session, and the average value for a specific run was used for calculating cells per ml.

Phytoplankton community structure was identified by pigment analysis. Four to 11 L seawater were filtered on Glass fiber filters (25 mm GF/F, Whatman, nominal pore size 0.7 μm). Excess water was removed from the filters by placing their underside on a kimwipe and then transferred to cryovials, stored in liquid nitrogen until arrival at the laboratory, and kept at −80°C until extraction. The collected cells were extracted in 1 ml 100% methanol for 2.5 h at room temperature. Extractions were immediately clarified with syringe filters (Acrodisc CR, 13 mm, 0.2 μm PTFE membranes, Pall Life Sciences) and transferred to UPLC vials. Samples were run on a UPLC and pigments identified based on retention time and spectrum absorbance. Several known standards [DHI Water and Environment Institute, (Hørsholm, Denmark)] were used to ease identification and calculate pigment concentrations (see Supplementary Information for details).

Seawater samples were collected from the Niskin bottles into 10 or 20 L polycarbonate carboys. Tubes for water transfer and the carboys were pre-rinsed three times with sample water. Samples for DNA (5–11.5 L) were pre-filtered through 11 and 5 μm nylon filters (Millipore), cells collected on 0.22 μm Sterivex filters (Millipore) and kept in storage buffer (40 mM EDTA, 50 mM Tris pH 8.3, and 0.75 M sucrose). Samples for RNA (0.5–4.5 L, according to volume filtered within 15 min from Niskin bottles being on-deck) were collected without prefiltration, directly on 0.2 μm filters (Supor-200 Membrane Disk Filters, 25 mm; Pall Corporation) to speed the filtration process and obtain enough biomass within a 15 min filtration time. While the lack of prefiltration might affect the perceived community composition, particle-associated bacteria are thought to be a relatively minor fraction of the community (see Supplementary Figure S2 in Mestre et al., 2017). Filters for RNA extraction were preserved in RNA Save (Biological Industries). DNA and RNA samples were stored in liquid nitrogen onboard and kept at −80°C in the laboratory until extraction. Nucleic acids were extracted at the BioRap unit, Faculty of Medicine, Technion, Israel using a semi-automated protocol, which includes manually performed chemical and mechanical cell lysis before the automated steps (see Supplementary Information for details).

Prior to reverse transcription, all RNA samples were tested by PCR for the presence of contaminating DNA. Total RNA was reverse transcribed using the iScript cDNA synthesis kit (Bio-Rad) according to the manufacturer’s instructions. A two-stage “targeted amplicon sequencing” protocol (e.g., Green et al., 2015) was performed to PCR amplify the 16S rRNA gene from cDNA to DNA (see Supplementary Information for a detailed description). The primers used in the first PCR stage consisted of the 16S primer set 515F-Y and 926R (Parada et al., 2016) that targets the variable V4-5 region with common sequence tags added at the 5′ end as described previously (e.g., Moonsamy et al., 2013). The first PCR stage was performed in triplicate, and the three PCR reactions were pooled after validation on 1% agarose gels. Subsequently, a second PCR amplification was performed to prepare libraries. These were pooled after a quality control and sequenced (2 × 250 paired-end reads) using an Illumina MiSeq sequencer. Library preparation and pooling were performed at the DNA Services (DNAS) facility, Research Resources Center (RRC), University of Illinois at Chicago (UIC). MiSeq sequencing was performed at the W.M. Keck Center for Comparative and Functional Genomics at the University of Illinois at Urbana-Champaign (UIUC).

Processing of the sequenced data followed the same pipeline as outlined in Roth Rosenberg et al. (2021). Briefly, after quality control of the obtained pair of fastq files for each sample, paired-end sequencing reads were merged using Flash (version 1.2.11) with the following parameters: max-overlap: 95; min-overlap: 85, with all other parameters set to default (Magoč and Salzberg, 2011). Merged reads were denoised and pre-processed using Dada2 (version 1.1.6; Callahan et al., 2016). Reads were trimmed to 400 nucleotides, while all other parameters were set to default values. Chimeric PCR products were removed using the “tableMethod” parameter set to “consensus.” Taxonomy was assigned to Exact Sequence Variants (ESVs) with the “classify.seqs” command in MOTHUR, the SILVA database (version 128), and an 80% identity cutoff (Schloss et al., 2009). ESVs identified non-prokaryotic, chloroplasts or mitochondria origin were removed. In total, the 56 samples had 2,678,521 quality sequences (range 21,488–75,136 sequences per sample, average: 478301). To avoid bias related to differences in library size, all libraries were rarified to 21,000 reads per sample, using the “rrarify” command implemented in the R Package Vegan (Oksanen et al., 2017). Read data were deposited in the NCBI SRA database under the project number PRJNA548664. Station 4 data are labeled N1200, a change made to ease the reading of the article.

Four negative controls from the first PCR were added randomly into the sequencing plates to monitor overall potential for cross contaminations in both PCR and sequencing. These negative control samples averaged 929 reads, compared to samples averaging 50,503 reads. Accordingly, contaminant DNA should poorly compete with sample DNA for amplification.

To assess within station variability and robustness of the observed trends, duplicates for DNA and RNA samples were collected and extracted for samples from station 4 for the summer 2015 cruise, the spring 2015 and 2016 cruises, and for station 1 for the spring 2016 cruise. We chose these samples to have a representation from different seasons, as well as from the most coastal and most offshore stations. Based on their Bray–Curtis dissimilarity, replicates were significantly closer to each other than to other stations from the same cruise (paired Wilcoxon signed-rank test, RNA: z = 3.9395 and p < 0.001; DNA: z = 4.1899 and p < 0.001).

The R package vegan version 2.5-7 (Oksanen et al., 2017) was used to calculate alpha diversity parameters (Chao1 richness, Shannon H′ diversity, and inverse Simpson), beta-diversity by pairwise calculation of the Bray–Curtis dissimilarity, and for non-metric multidimensional scaling analysis. Vegan was also used for non-metric multidimensional scaling analysis based on the Bray–Curtis dissimilarity matrix with k = 2, 100 runs and 100 iterations. Significance of the grouping factors molecule type (RNA vs. DNA), season (early winter, spring, and summer), and sampling station (stations 1–4) was assessed with the Vegan “adonis” command was used based on Bray-Curtis dissimilarities with 999 permutations.

Correspondence between abiotic and biotic measurements and variation in community composition was examined by variation partitioning analysis and distance-based redundancy analysis. First, the proposed explanatory variables were divided into three matrices: (I) physical: distance from shore, temperature, salinity, and turbidity; (II) nutrients: phosphate, NO₃ + NO₂, and silicate; and (III) biotic: total fluorescence as proxy for chlorophyll, picoeukaryote counts, and total cell counts. Variation partitioning calculates the canonical coefficient of determination (adj. R²), representing the forecasting potential of multiple regression between the response variables (i.e., the ESVs) and explanatory variables (i.e., environmental measurements). Variation partitioning was performed using Vegan “varpart” command. Permutation tests, to validate the significance of the coefficients (adj. R²), were executed both directly and conditionally (excluding variation that is related to interactions between the different explanatory variable matrices) by calculating distance-based redundancy analysis (dbRDA) using the “dbRDA” command followed by “anova” for the model and the tested variables.

Seasonal and spatial preferences within a season of ESVs were analyzed with the “multipatt” command with 9,999 permutations of the R package indicspecies version 1.7.7 (De Cáceres et al., 2010).

In addition to statistical tests described above, we used paired Wilcoxon signed-rank tests to test if Bray–Curtis dissimilarity between replicates from the same sample location was lower than between samples from different locations from the same cruise. A Kruskal–Wallis test was used to find significant differences in alpha diversity indices between seasons, followed by Dunn’s post-hoc test with Bonferroni correction of value of p. Two-tailed Mann–Whitney U tests were used to determine for each season if Bray–Curtis dissimilarity between station 1 and the other three stations was higher than between the three further offshore stations. Spearman correlations were performed to test for correlation between station distance and Bray–Curtis dissimilarity between samples from the same cruise. Spearman correlations were also used to test for correlation between DNA and RNA abundances at the family level. All these tests were performed in PAST version 4.05 (Hammer et al., 2001).

Satellite (Figure 1) and CTD data (Supplementary Figure S1) indicated that our sampling times occurred in three distinct stages of the annual mixing and stratification cycle in the EMS (Kress and Herut, 2001; Ben Ezra et al., 2021): (i) the stratified water column in summer, (ii) early mixing and phytoplankton bloom in early winter, and (iii) the declining surface water phytoplankton bloom in spring after its usual peak in January/February (Reich et al., 2021; Roth Rosenberg et al., 2021). Pigment data indicated that these stages differed in their phytoplankton community. Early winter samples were dominated by haptophytes of either Phaeocystaceae or members of the Prymnesiaceae and Isochrysidaceae families. Spring samples were characterized by peak abundances of Prochlorococcus (Supplementary Figure S5D; Supplementary Table S6), and summer samples had the lowest Synechococcus abundances at the sampling depth at stations 2–4 (Supplementary Figure S5C; Supplementary Table S6; Supplementary Information includes further details on the pigment analysis and results).

The ultra-oligotrophic status of the Levantine basin with extremely low nutrient concentrations in its surface waters (Krom et al., 2005; Ben Ezra et al., 2021; Reich et al., 2021) was consistently observed in our samples (Supplementary Tables S2 and S3). Phosphate concentrations were near or below our practical detection limit (40 nM) in almost all cruises. Nitrate + nitrite concentrations for these samples ranged from <50 to 240 nM at the most offshore station and from <50 to 690 nM at the most coastal station. At the most coastal station, winter samples were clearly higher in nitrate + nitrite concentrations than spring and summer samples, with no clear difference between the latter two. At the further offshore stations, no systematic difference in nitrate + nitrite could be found between winter and spring samples, but summer samples were always at or below our detection limit. In early winter nitrite + nitrate concentrations decreased from station 1 toward the offshore stations, whereas in spring, they increased from station 1 to station 3 before decreasing again toward station 4, the most offshore station. Additional analysis of fresh (unfrozen) samples provided a lower detection limit (Ben Ezra et al., 2021) and was carried out in the summer 2016 cruise (Supplementary Table S3), one of the most nutrient depleted sampling periods. These data show that, although most of our nutrient measurements fell below the detection limit, the sampling transect was able to capture a nutrient gradient impacted by the proximity to the shore.

Other environmental parameters varied with different degrees of change. For instance, surface seawater temperature varied significantly with the seasons, with the lowest temperatures observed in spring (17.7°C–18.6°C) and highest in summer (26.2°C–28.7°C). However, the seasonal changes in salinity were small (<0.7 PSU from lowest to highest value; Supplementary Table S2), suggesting minimal water mass intrusions from other locations.

Overall, environmental parameters changed more between seasons than between station locations within a cruise. The average difference between minimal and maximal temperature across seasons at a sampling site was 10.7°C compared to only 0.3°C between sites within a given cruise. The same could be observed for salinity (0.59 vs. 0.14 PSU), fluorescence (0.36 vs. 0.08), and turbidity (4.4 vs. 0.54 FTU). Within the nutrient data, these differences were smaller. Nitrate + nitrite differed on average 0.44 μmol/L between the lowest and highest value at a given site between seasons compared to an average difference of 0.3 μmol/L between the highest and lowest value within a cruise. For silicate, this was 0.62 μmol/L compared to 0.54 μmol/L. We did not compute these differences for phosphate as values were frequently below the detection limit.

Along the transect the collected environmental data differed the most between station 1, the shallowest and most coastal sampling site (16 km from shore), and the other three more offshore stations (25–50 km from shore). Station 1 often had higher nutrient concentrations (e.g., silicate in spring and nitrite + nitrate in early winter; Supplementary Tables S2 and S3) and chlorophyll concentrations (Figure 1) throughout the year. The station closest to the shore also differed from the rest in temperature (higher in early winter and spring), salinity (lower in summer), fluorescence (higher in early winter), and turbidity (higher in early winter and summer 2016, lower in summer 2015; Supplementary Table S5; Supplementary Figure S1).

The best-represented microbial groups were recurrently identified throughout the study. However, their proportions and levels of activity (retrieved as community 16Sr RNA) showed some variation at different sampling locations and times (Figure 2). Bacteria dominated the microbial community at the 10 m sampling depth throughout the year at all stations. After subsampling, they made up 99.7% of all reads (3,240 ESVs) compared to 0.3% of archaeal reads (39 ESVs). Most bacterial reads (82.6%) belonged to one of three classes: Alphaproteobacteria (47.0% of all reads), Gammaproteobacteria (18.5%), and Cyanobacteria (17.1%). The Alphaproteobacteria were dominated by the SAR11 clade (57.4% of the Alphaproteobacteria), the Cyanobacteriia by Prochlorococcus and Synechococcus (together 94.1% of the Cyanobacteriia), and the Gammaproteobacteria by the SAR86 clade (53.3% of the Gammaproteobacteria). Only four other classes contributed more than 1% of total reads: Bacteroidia (6.0%), Verrucomicrobia (3.3%), Marinimicrobia (SAR406 clade) (2.4%), and Acidimicrobiia (1.4%). The archaea were dominated by marine group II Thermoplasmatota (68.4% of all archaeal reads), followed by Thermoplasmatota unidentified at the family level (21.0%) and Nitrosopumilaceae Thaumarchaeota (10.6%). Figure 2 shows the relative abundance of the 18 families representing >1% of either the total DNA (Figure 2A) or RNA (Figure 2B) reads.

Comparison of the microbial community structure using 16S rRNA allowed us to further identify the most active microbial groups in these communities. For 15 of the 18 families shown in Figure 2, we found a significant positive correlation (Spearman correlation, r > 0.6, p < 0.05) between relative abundance in DNA and corresponding RNA samples (Supplementary Table S7), indicating a strong link between presence and activity. However, the average RNA:DNA ratios differed between family-level taxonomic groups (e.g., 0.36 ± 0.09 for SAR11 clade I and 2.48 ± 0.91 for Cyanobiaceae) and the composition of resident (DNA) and active (RNA) microbial communities clearly differed (Figure 2), as supported by non-metric multidimensional scaling analysis (NMDS; Figure 3) and Adonis test (F = 12.9, R² = 0.219, p < 0.001). Given this difference, we performed all subsequent analyses of microbial community composition (DNA) and activity (RNA) separately.

Sampling season had a clear impact on both the microbial community structure (DNA) and activity (RNA; Figure 3). The effect was significant and explained 54.9 and 55.9% of the variability in RNA and DNA samples, respectively. In contrast, the effect of station location across the coastal to offshore transect on the observed variability was small (8.7% and 7.6% in RNA and DNA samples, respectively) and not significant (Adonis test, DNA: season R² = 0.55944, station R² = 0.07566; RNA: season R² = 0.54948, station R² = 0.08692; p < 0.001 for season and p > 0.05 for station in both analyses, see details in Supplementary Table S8).

In addition to microbial community structure, the overall microbial diversity also changed significantly with season (Supplementary Figure S3). For DNA samples, the three alpha diversity indices (Shannon H, inverse Simpson, Chao1) indicated the lowest diversity in summer (Kruskal–Wallis test, p < 0.001 followed by Bonferroni corrected Dunn’s post-hoc test p < 0.05), when the system was most oligotrophic, with no difference between winter and spring samples. All three indices also varied significantly seasonally in the RNA samples (Kruskal–Wallis test, p < 0.01, followed by Bonferroni corrected Dunn’s post-hoc test p < 0.05), but the pattern was less clear.

Given the particularly high abundance of members of the SAR11 group in the EMS in our samples (average of 36.9% of DNA reads) and in previous studies (e.g., Dubinsky et al., 2017), we further analyzed the seasonality of its different clades. The SAR11 clade Ia was the dominant SAR11 clade contributing on average 61.3% of all SAR11 DNA reads (range 41.0%–90.2%). It was more abundant during the summer months than in spring and winter of the two sampling years (Supplementary Figure S4). In contrast, an opposite trend was observed for SAR11 clades Ib, II, and IV, for which relative abundances were lowest during the summer, especially at the more oligotrophic stations 2–4 (Supplementary Figure S4). These data support the notion that different SAR11 clades potentially occupy distinct ecological niches in this ultra-oligotrophic environment.

Other relatively abundant groups with clear seasonal patterns included the AEGEAN-169 marine group of the Alphaproteobacteria family Rhodospirillaceae as well as the SAR86 and SAR116 clades (all most abundant in summer; Figure 2). Within the family Cyanobiaceae, Prochlorococcus showed a clear peak in spring, whereas no clear pattern was seen for Synechococcus (Supplementary Figure S5).

To better understand the influence of different parameters on microbial structure and activity, we examined three environmental data matrices: (i) physical data (distance from shore, temperature, salinity, and turbidity), (ii) biological data (fluorescence as a proxy for chlorophyll, picoeukaryotic phytoplankton, and total microbial cell number), and (iii) nutrient data (phosphate, nitrate + nitrite, and silicate). Together, the three matrices explained about 67% of the observed variability in both the microbial structure (DNA samples) and the microbial activity (RNA samples). Partition analysis indicated that in both analyses, the physical data matrix explained most of the variability (44.5% and 49.4% for DNA and RNA, respectively), followed by the biological data matrix (30.2% and 28.3% for DNA and RNA, respectively) and the nutrient data matrix (17.1% and 13.0% for DNA and RNA, respectively). The physical data matrix also had the largest percent of variability not explained by any of the other two matrixes (DNA: 22.6%; RNA: 24.2%), whereas the nutrient matrix had the lowest (DNA: 0.4%; RNA: 1.1%; Supplementary Figure S6). This suggested that particular environmental factors within these data groups might be responsible for the different levels of prediction. To identify those key parameters, we then performed a distance-based redundancy analysis. As in the partition analysis, the physical matrix was able to explain the largest amount of variance, followed by the biological matrix and the nutrient matrix. Yet, the percent variation explained by the nutrient matrix was not significant after correcting for the other two matrices. Overall, environmental factors that tended to change with season, such as temperature and salinity from the physical matrix, and fluorescence and total cell count from the biological matrix were factors with a significant level of prediction in the DNA and RNA analysis (both in the unconditioned and conditioned analysis, p < 0.05, Supplementary Figure S6).

We used an indicator species analysis to identify ESVs with seasonal preferences. Significant association with one or two seasons was found in 299 ESVs in the DNA samples representing 77.5% of all DNA reads (Supplementary File 1—sheet 2). In the RNA samples, 386 ESVs representing 67.3% of all RNA reads had significant associations with one or two seasons (Supplementary File 1—sheet 3). When we compared the 202 ESVs that had significant seasonal preferences in both the DNA and RNA analyses, 181 ESVs had the same seasonal preferences in both analyses and another 20 ESVs had an additional season in one of the analyses.

The indicator analyses identified several taxonomic groups favoring specific seasons (see Supplementary File 1—sheets 2 for DNA and sheet 3 for RNA indicator ESVs): for example, all indicator Alteromonadaceae ESVs were specific for summer; several taxonomic groups were (mostly) indicator for spring (e.g., most ESVs of Prochlorochoccus, SAR11 clade Ib and II, all ESVs of Actinomarinaceae, and SAR11 clade IV) and some for early winter (e.g., most indicator ESVs of the Chloroflexi clade SAR202, the Deltaproteobacteria clade SAR324, and the Gammaproteobacteria clade OM182). However, most taxonomic groups had indicator species from several seasons: for example, most ESVs from SAR11 clade Ia were either associated with summer or spring at the DNA level. The Marinimicrobia clade SAR406 had ESVs indicator for spring and early winter, but not for summer. Likewise, the significant Synechococcus ESVs associated with spring, spring and early winter, or early winter, but not summer. ESVs from the SAR86 clade associated with summer, spring, and early winter at the DNA level, with most of the abundant ESVs being indicators for summer.

To analyze which ESVs dominate the community and whether dominant and rare ESVs show different patterns, we classified ESVs into four categories broadly following the categories in Auladell et al. (2021): dominant (present in ≥75% of the samples with ≥1% relative abundance), rare (always <1%), intermediate persistent (present in ≥75% of samples), and intermediate fluctuating (present in <75% of samples). Despite the large influence of season on the microbial community, the microbial community structure was dominated by only a few bacterial ESVs throughout the year (Table 1). The eight dominant ESVs in the DNA data made up together on average 38.5% of the community and were especially dominant in summer when they averaged 48.4%. Four of these ESVs belonged to SAR11 clade Ia and one each belonged to SAR11 clade II, AEGEAN-marine group, Ca. Actinomarina and Synechococcus. Especially noteworthy was ESV000001, a SAR11 clade Ia strain, which was the most abundant ESV in all but three DNA samples and made up on average 11.6% of the community.

In the RNA data, only five ESVs fell into the dominant group: three SAR11 clade Ia, a Synechococcus ESV—all of which were also dominant in the DNA data—and a SAR86 ESV. Together these ESVs represented on average 19% of the RNA community data (Table 1). The Synechococcus ESV (ESV000003) was the most active one across most samples and represented on average 8.7% of the RNA reads per sample. The only other ESV contributing always more than 1% to the RNA reads was ESV0000001 of SAR11 clade Ia, which was the most abundant ESV in the DNA data. ESV0000002, a Prochlorococcus strain, that was in the intermediate persistent group, had the highest overall abundance (32.6% of the RNA reads), was among the top 5 most abundant ESVs in 12 of the 24 samples (average RNA relative abundance 8.69%), and had a clear seasonality. Together these three strains represented on average 22% and maximal 40.9% of the RNA reads.

Distance from shore was not a significant factor when analyzing the effect of physical parameters on the community (Supplementary Figure S6). We also found no significant correlation of distance between stations and Bray–Curtis dissimilarity within cruises (neither when all cruises were analyzed together or cruises were split by season; Spearman correlation, p > 0.05). However, samples from the most coastal station 1 grouped apart from the samples of the more offshore stations 2–4 in the NMDS plot in early winter and summer, but not in spring (Figure 3).

As station 1 also differed in environmental parameters from the three more offshore stations (see environmental settings above), we tested if this difference could be reflected in the microbial community structure (DNA) and activity (RNA). In both early winter and summer within cruises, Bray-Curtis dissimilarities were significantly higher between station 1 and the other three stations than among the three offshore stations, while no difference was observed in spring (Figure 4).

Indicator species analyses performed for each season separately found 92 ESVs in the DNA data exhibiting a spatial preference (either station 1 or station 2–4), of which 13 showed significant association in multiple seasons with consistent spatial preference (details for each ESV are given in Supplementary File 1 sheets 4). Similar results were found in the active community as 129 ESVs showed a significant spatial preference in the RNA samples (Supplementary File 1 sheet 5).

Overall, several taxonomic groups showed a spatial preference among their indicator ESVs (Supplementary File 1 sheets 4 and 5). Indicator ESVs of archaea, SAR11 clade IV, SAR202 Chloroflexiota, Bdellovibrionaceae clade OM27, Marinimicrobiota SAR406, Planctomycota, and most Verrucomicrobiota (except ESVs of Lentimonas and R76-B128) were associated with the coastal station 1. In contrast, several ESVs of the SAR11 clade I and most of SAR11 clade II, SAR86, and SAR116 were indicators for the offshore stations 2–4 (Supplementary File 1 sheets 4 and 5). Within Cyanobacteria Cyanobium ESVs were associated with station 1, while Prochlorococcus and Synechococcus ESVs were associated with the offshore stations 2–4.