Functional Seasonality of Free-Living and Particle-Associated Prokaryotic Communities in the Coastal Adriatic Sea

Steiner, Paul A.; Geijo, Javier; Fadeev, Eduard; Obiol, Aleix; Sintes, Eva; Rattei, Thomas; Herndl, Gerhard J.

doi:10.3389/fmicb.2020.584222

ORIGINAL RESEARCH article

Front. Microbiol., 16 November 2020

Sec. Aquatic Microbiology

Volume 11 - 2020 | https://doi.org/10.3389/fmicb.2020.584222

Functional Seasonality of Free-Living and Particle-Associated Prokaryotic Communities in the Coastal Adriatic Sea

$\r\nPaul A. Steiner*$ Paul A. Steiner^1*

Javier Geijo²

Eduard Fadeev¹

Aleix Obiol³

Eva Sintes⁴

Thomas Rattei²

Gerhard J. Herndl^1,5

¹Department of Functional and Evolutionary Ecology, University of Vienna, Vienna, Austria
²Department of Microbiology and Ecosystem Science, Division of Computational Systems Biology, University of Vienna, Vienna, Austria
³Institut de Ciències del Mar, Institut de Ci ncies del Mar – Consejo Superior de Investigaciones Cient ficas (ICM-CSIC), Barcelona, Spain
⁴Instituto Español de Oceanografia, Centre Oceanogràfic de les Balears, Palma, Spain
⁵Royal Netherlands Institute for Sea Research (NIOZ), Department of Marine Microbiology and Biogeochemistry, Utrecht University, Utrecht, Netherlands

Marine snow is an important habitat for microbes, characterized by chemical and physical properties contrasting those of the ambient water. The higher nutrient concentrations in marine snow lead to compositional differences between the ambient water and the marine snow-associated prokaryotic community. Whether these compositional differences vary due to seasonal environmental changes, however, remains unclear. Thus, we investigated the seasonal patterns of the free-living and marine snow-associated microbial community composition and their functional potential in the northern Adriatic Sea. Our data revealed seasonal patterns in both, the free-living and marine snow-associated prokaryotes. The two assemblages were more similar to each other in spring and fall than in winter and summer. The taxonomic distinctness resulted in a contrasting functional potential. Motility and adaptations to low temperature in winter and partly anaerobic metabolism in summer characterized the marine snow-associated prokaryotes. Free-living prokaryotes were enriched in genes indicative for functions related to phosphorus limitation in winter and in genes tentatively supplementing heterotrophic growth with proteorhodopsins and CO-oxidation in summer. Taken together, the results suggest a strong influence of environmental parameters on both free-living and marine snow-associated prokaryotic communities in spring and fall leading to higher similarity between the communities, while the marine snow habitat in winter and summer leads to a specific prokaryotic community in marine snow in these two seasons.

Introduction

Marine snow, described as detrital particles larger than 500 μm, plays an important role in the export of organic carbon to the deep sea and the sequestration via the biological carbon pump (Cho and Azam, 1988; Ducklow et al., 2001; Herndl and Reinthaler, 2013). Chemical characterization, microbial community composition, and activity of marine snow in a large variety of oceanic regions have been extensively studied (Alldredge and Silver, 1988; Simon et al., 2002; Turner, 2002, 2015). However, the majority of studies represent snapshots of the biotic and abiotic characteristics of marine snow captured with sediment traps (De La Rocha and Passow, 2007; Buesseler et al., 2008), marine snow catchers (Riley et al., 2012; van der Jagt et al., 2018), SCUBA divers (Kaltenböck and Herndl, 1992; Müller-Niklas et al., 1994; Rath et al., 1998) or filtration systems (Salazar et al., 2015; Mestre et al., 2017). Only a few studies investigated the dynamics of marine snow and its associated microbial community over a few weeks (Kaltenböck and Herndl, 1992; Müller-Niklas et al., 1994). Recently, the importance of increasing the spatial and temporal resolution of research studies to understand the dynamics in the community composition and functioning of ecosystems has been emphasized (Chow et al., 2013; Fuhrman et al., 2015; Ward et al., 2017). Marine snow-associated microbial communities are considered to be relatively insensitive to changes in environmental parameters of the bulk seawater. With increasing size and volume of the individual marine snow, the physico-chemical conditions determine the composition and activity of the marine snow associated microbiome (Vojvoda et al., 2014; Yung et al., 2016).

Free-living microbial communities are exposed to rapidly changing environmental parameters, resulting in seasonal changes in taxa and closely related ecotypes (Eren et al., 2013; Fuhrman et al., 2015; Yung et al., 2015; Ward et al., 2017; Steiner et al., 2019). Additionally, organic matter leaching from marine snow (Stocker, 2012) links and shapes the free-living community composition (Vojvoda et al., 2014). In contrast, the marine snow-associated microbial community is primarily affected by the origin, composition and developmental stage of the marine snow, and secondarily by the seasonal changes of environmental parameters (Vojvoda et al., 2014; Yung et al., 2016; Duret et al., 2019).

Marine snow represents a specific microhabitat with organic and inorganic nutrient concentrations potentially orders of magnitude higher than in the ambient water (Alldredge and Silver, 1988; Kaltenböck and Herndl, 1992; Müller-Niklas et al., 1994). In addition, it also provides plenty of surfaces for microbial attachment. Therefore, it is not surprising that the prokaryotic community associated with marine snow is distinctly different from that of the ambient water under such conditions (Kaltenböck and Herndl, 1992; Rath et al., 1998; Simon et al., 2002; Simon et al., 2014). However, these reports have been performed during extensive marine snow formation mainly in the summer months of temperate seas. Whether these differences between the free-living and marine snow-associated prokaryotic community composition are maintained throughout the seasonal cycle in temperate seas, however, remains unknown.

The northern Adriatic Sea is a temperate, coastal system with a high abundance of marine snow receiving considerable attention since the 1980’s (Herndl and Peduzzi, 1988; Herndl, 1992; Kaltenböck and Herndl, 1992; Rath et al., 1998). The extent of marine snow formation has been suggested to be due to changes in temperature and nutrient concentrations as well as to the specific water mass circulation patterns in the northern Adriatic Sea (Cozzi et al., 2004). We hypothesized that marine snow offers more stable environmental conditions resulting in a more stable community over the two seasons than the ambient water for the free-living microbial community. Thus, we analyzed the dynamics of the composition of the free-living and particle-associated microbial community based on 16S rRNA gene amplicon sequences over a seasonal cycle in the northern Adriatic Sea. Furthermore, the functional potential of marine snow-associated and free-living prokaryotic community at two contrasting seasons (summer 2015 and winter 2016) was assessed using metagenomics analysis. Finally, the eukaryotic and prokaryotic community composition was related to the environmental parameters and the implications of the prokaryotic metabolic potential for biogeochemical cycles are discussed.

Materials and Methods

Sampling

Samples were collected at twelve days covering a 1.5-year period from April 2015 to July 2016 in the coastal northern Adriatic Sea off Rovinj, Croatia (45.08347°N, 13.60518°E). Marine snow was collected with 100 mL sterile syringes at 15 ± 2 m depth by SCUBA diving and processed as described elsewhere (Steiner et al., 2019). In brief, the marine snow collected in several syringes was pooled into a 0.1 M HCl-rinsed glass bottle on deck and stored at in situ temperature in the dark. The marine snow included the pore water and a small amount of the water surrounding marine snow. The ambient water was collected at 15 m depth with 5 L Niskin bottles. The samples were transferred to the Center for Marine Research of the Ruder Bošković Institute at Rovinj within 30 min after sample collection. Temperature, salinity and in vivo fluorescence representing chlorophyll a concentrations were determined monthly with a SBE25 conductivity-temperature-depth probe (Sea-Bird Electronics, WA, United States) at the long-term monitoring station RV001 as described in Ivancic et al. (2018).

DNA Extraction

Marine snow (150 mL to 500 mL) and 1 L of ambient water were filtered onto 0.2 μm polyethersulfone filters (47 mm diameter, Supor, PALL Gelman) using an aspirator pump (Cole-Parmer). Filters were placed in cryovials (Biozym), flash-frozen in liquid nitrogen and stored at −80°C. DNA was extracted using a standard phenol-chloroform method as described in Steiner et al. (2019). Hereinafter, ambient water and marine snow are abbreviated as ambient water (AW) and marine snow (MS), respectively. Metagenome samples are labeled indicating the season, winter (W) or summer (S), followed by AW or MS, and the duplicate sample number 1 or 2. Duplicates collected in winter (W) correspond to samples collected on 1 and 2 February 2016 and summer (S) duplicates correspond to samples collected on 24 and 25 June 2015. All sequence data are publicly available at the European Nucleotide Archive (ENA) at EMBL-EBI under accession number PRJEB38662¹. Data was deposited using brokerage service of the German Federation for Biological Data (GFBio; Diepenbroek et al., 2014) in compliance with the Minimal Information about any (x) Sequence (MIxS) standard (Yilmaz et al., 2011).

Prokaryotic Community Composition Using Amplicon Sequencing

Twelve MS and 11 AW samples were collected for amplicon sequencing of the 16S rRNA gene at the same location as the metagenomics samples. Water samples were collected seasonally over 1.5 years and processed as described elsewhere (Steiner et al., 2019). The primers 341_ill forward (TCGTCGGCAGATGTGTATAAGAGACAGCCTACGGGNGG CWGCAG) and 802_ill reverse (GTCTCGTGGGCTCGGAGA TGTGTATAAGAGACAGGACTACHVGGGTATCTAATCC) targeted the V3–V4 hypervariable region (∼460 bp) of the 16S rRNA of bacteria and archaea. In brief, PCR amplification was performed with the primers 341_ill forward and 802_ill reverse containing adaptors and a KAPAHiFi Mastermix (Peqlab) using the following program: initial denaturation at 94°C for 3 min, followed by 20 cycles of 94°C for 30 s, annealing at 56°C for 30 s and extension at 72°C for 90 s. A final extension step was carried out at 72°C for 7 min, followed by cooling at 4°C. PCR products were purified using Agencourt AMPure XP magnetic particles (Beckman Coulter) and quantified with a Quant-IT PicoGreen^® Assay (Invitrogen). Subsequently, Nextera PCR was performed with the same thermocycling conditions as described above for 10 additional cycles, followed by pooling and sequencing on an Illumina MiSeq system (Microsynth AG, Balgach, Switzerland) using v2 chemistry. Amplicon sequences were demultiplexed and trimmed by Microsynth AG (Balgach, Switzerland). All samples were sequenced with a coverage of 1.0 and the number of reads ranged from 7,256 to 6,4801 per sample.

Bioinformatics analysis followed the ‘DADA2 1.12’ (Callahan et al., 2016) pipeline tutorial with specific adaptations when filtering and trimming: truncLen = c(225, 210), maxN = 0, minQ = 2, maxEE = c(3, 3), truncQ = 0. Singletons, mitochondria and chloroplast sequences were removed. Statistical analysis of the amplicon sequences was conducted using R² in RStudio v.1.2.5033 (RStudio-Team,, 2015) with the package phyloseq v.1.30.0 (McMurdie and Holmes, 2013). The obtained amplicon sequence variants (ASVs) were variance stabilized (McMurdie and Holmes, 2014) and a differential gene expression analysis was conducted using the R package DESeq2 v.1.26.0 (Love et al., 2014). Differences in variance stabilized abundance of ASVs were considered significant if p < 0.05. Only ASVs with a log₂ fold change higher than 1 in sequence abundance were further analyzed.

Metagenomics Analysis of Prokaryotes

Metagenomic DNA was quantified fluorometrically using the Quant-iT^TM PicoGreen^® Assay (Invitrogen) on a plate reader (Infinite 200). The library preparation of metagenomic DNA was done with the Westburg kit and sequencing was performed on an Illumina HiSeq 2,500 (2 × 150 paired-end platform) at the Vienna BioCenter Core Facilities GmbH (VBCF). Metagenome sequences were quality checked using fastQC v.0.11.8. Sequences were trimmed with trimmomatic v.0.39 (Bolger et al., 2014) and a sliding window size of 3 bases and a quality of 30. Assembly was done with megahit v.1.1.2 (Li et al., 2015).

The program phyloFlash v.3.3 (Gruber-Vodicka et al., 2019) was used to analyze the prokaryotic community composition based on 16S rRNA obtained from metagenomes with the following settings: -lib LIB -emirge -poscov -treemap. The correlation between the 16S rRNA gene amplicon data matrix and the extracted 16S rRNA genes from the metagenome was tested using Mantel test implemented in the R package vegan v.2.5.6. Prokaryotic OTUs of eight metagenomes were correlated to metataxonomic ASVs determined in the corresponding MS and AW samples collected on 24 and 25 June 2015 and 1 and 3 February 2016.

Further taxonomic and functional analyses were performed with MEGAN Community Edition v.6.17.0 (Huson et al., 2007) with default settings (min score: 50, max expected: 0.01, min percent identity: 0, top percent: 10, min support percent: 0.05, min support: 0, LCA algorithm: naive, percent cover: 100, read assignment mode: readCount). Reads were blasted against the NCBI nr database using DIAMOND with short reads settings: –b 12.0 –k 1 –f 100 –e 0.00001 –p 20. The taxonomic affiliation was determined by aligning the reads to the database prot_acc2tax-Jul2019 × 1.abin. The database acc2eggnog-Oct2016X.abin was used for functional assignment in MEGAN6. Reads putatively assigned to prokaryotes were extracted from the MEGAN6 file using the ‘extract to new document’ function. The prokaryotic reads were visualized using taxonomy and EGGNOG viewer (Powell et al., 2014) of MEGAN6 in absolute comparison mode. Taxonomic and gene assignments to clusters of orthologous groups (COGs) were extracted using the ‘extract to new document’ function to determine the taxonomic affiliation of the COGs. Subsequently, prokaryotes were sorted out from eukaryotic and other reads. The ‘export’ function was used to extract absolute numbers of reads assigned to COGs. Absolute read counts were variance stabilized (McMurdie and Holmes, 2014) and differential gene expression analysis was conducted using the R package DESeq2 v.1.26.0 (Love et al., 2014). COGs were selected for further analysis if differences in their variance-stabilized abundance were significant (p < 0.05) and the log₂ fold change of sequence abundance was > 1. COGs with assignment to transport function were extracted from the variance-stabilized data set and used for further analyses.

Marker gene analysis was performed based on the Kyoto Encyclopedia of Genes and Genomes Orthology (KOs) functional classification of metabolisms (Supplementary Table S1). The presence of 69 selected key marker genes (Acinas et al., 2019) was checked against the aligned metagenomic reads. Predicted genes of aligned reads were determined with MetaProdigal v.2.6.3 (Hyatt et al., 2012) and quantified using FPKM (Fragments Per Kilobase Million) normalization.

Metagenomics Analysis of Eukaryotes

The eukaryotic community composition was analyzed as described elsewhere Obiol et al. (2020). Briefly, a BLAST search against a subset of the eukaryotesV4 database with > 90% similarity and > 70% alignment coverage was performed to retrieve fragments of the 18S rRNA V4 region from the metagenomes. The extracted fragments were mapped to the full eukaryotesV4 database using USEARCH v.11 local alignment at 97% identity. All hits with maximal score for every fragment were selected. All top hits of each read were merged into one, keeping the minimum consensus taxonomic level for taxonomic classification. Reads which resulted in no consensus taxonomic level remained unclassified.

Analyses of the Structure of the Microbial Community

Diversity indices were calculated using the package vegan v.2.5.6 in R on rarefied data to the lowest number of ASVs. Diversity indices, one-way ANOSIM, PERMANOVA, SIMPER and principal component analysis (PCA) were calculated in PAST3 (Hammer et al., 2001) using Euclidean similarity index. Graphs were created with the program SigmaPlot 14.0 and Venny 2.1 (Oliveros, 2015).

Results and Discussion

16S rRNA Amplicon Sequence Analysis

Taxonomic Analysis

On average 891 ± 232 ASVs were hit by 46797 ± 13521 sequences in the AW and 437 ± 219 ASVs were hit by 31548 ± 17230 sequences in MS (Supplementary Figure S1). In total 1711 ASVs had a prevalence of more than 10%, implying the presence of these ASVs in at least two samples.

Diversity Indices

Species richness was on average 1.5 ± 0.6 times higher in the AW than MS in spring, summer and fall (Supplementary Table S2 and Supplementary Figure S2). In the winter, however, species richness in the AW was 5.4 ± 6.9 times higher than in MS (Supplementary Table S2 and Supplementary Figure S2). The observed low species richness of MS-associated communities as compared to AW is in agreement with previous findings and was linked to the high nutrient concentrations these microhabitats offer (DeLong et al., 1993; Acinas et al., 1999). The noticeably lower species richness in MS in winter compared to the other seasons and to AW concurred with lower Shannon diversity (3.3 ± 0.2) as compared to the average diversity in MS (4.9 ± 0.4; Supplementary Figure S2). This is likely due to the dominance of Synechococcus ASVs in the winter MS communities (Supplementary Figure S3).

Prokaryotic Community Analysis

The core community of each month, resembled by the number of prokaryotic classes found in both, the AW and MS, ranged between 20 and 23 in spring, summer and fall. In winter, the core community was made up by a lower number of classes (16 classes) compared to the other seasons (Supplementary Figure S4). However, the number of bacterial classes present exclusively in the AW was relatively high in winter (14 classes as compared to, on average, 4.7 ± 3.4 classes in spring, summer, and fall; Supplementary Figure S4) and comprised 1.4% of all sequences. The largest part thereof was made up by Euryarchaeota, Dadabacteria, and Patescibacteria. These classes contain ecotypes with seasonal patterns, surface water and terrestrial or groundwater clades (Hugoni et al., 2013; Herrmann et al., 2019; Graham and Tully, 2020). The fraction of bacterial classes occurring exclusively in either the AW or MS was on average 0.05 ± 0.04% in spring, summer, fall and in MS in winter, indicating a minor role of habitat specialists in terms of sequence proportions.

Summer-spring communities were distinct from winter-fall prokaryotic communities, with the free-living communities showing larger seasonal dissimilarities than the MS-associated prokaryotic communities (Figure 1). This suggests a more pronounced impact of seasonality on the free-living than on MS-associated prokaryotes possibly due to the more stable microenvironmental conditions in MS than in the AW (Vojvoda et al., 2014; Yung et al., 2016).

FIGURE 1

Figure 1. Principal component analysis showing the seasonal differences in the marine snow (MS)-associated and ambient water (AW) prokaryotic community composition based on amplicon sequencing (symbols represent data collected on 12 and 11 days, respectively) labeled according to the season. The percentage of variation explained by each axis is indicated in parenthesis.

Prokaryotic Taxa Enriched in Marine Snow and Ambient Water

Overall, the free-living and MS-associated communities differed significantly (one-way ANOSIM p < 0.05; Figure 1). Yet, the low number of samples precluded to statistically test whether the free-living and MS-associated communities differed in particular seasons. However, ASVs with significantly different sequence abundances between seasons were obtained. The AW sample collected at 30 July 2016 was identified as an outlier and excluded from the differential sequence abundance analysis (DESeq2) due to high compositional similarity to the summer MS community (Figure 1).

In winter, 589 ASV were enriched in the AW as compared to MS, and 19 ASVs in the MS as compared to the AW. Enriched ASVs in the AW represented 91% of the total AW community in winter, enriched ASVs in MS accounted for 26% of the total MS community. Twelve of the 589 AW-enriched ASVs had a mean abundance > 1% and were assigned to Synechococcales, SAR11, Thiomicrospirales, SAR86, SAR406, and Actinomarinales (Figure 2A). These taxa comprise groups typically found in the free-living community (Mestre et al., 2017; Preston et al., 2019), some consisting of many ecotypes, including ecotypes characteristic for the winter in temperate marine waters (Giovannoni, 2017; Hoarfrost et al., 2020). Seven MS-ASVs exhibited a mean relative abundance > 1% in winter (Figure 2A). These MS-ASVs were affiliated to orders including diverse phenotypes, such as Caulobacterales, capable of attachment and motility (Poindexter, 1964), Obscuribacterales, an order of non-phototrophic fermentative and motile cyanobacteria (Soo et al., 2014), the symbiotic nitrogen fixing Rhizobiales (Summers et al., 2013), pigmented photoautotrophs of the order Sphingomonadales (Siddaramappa et al., 2018), and Propionibacteriales, usually found on human skin (Barka et al., 2016) and in liver and kidney tissue of marine fish (Meron et al., 2020; Figure 2A). The high abundance of potential symbionts and pathogens in MS might be related to the influx of waste water from a fish cannery about 1.5 km away the sampling location (Paliaga et al., 2017). No ASVs were found enriched in spring AW and MS, indicating freshly generated and hence, recently populated MS, or a strong effect of environmental parameters on both AW and MS communities. The slight increase in seawater temperature and photosynthetically active radiation (Trisolino et al., 2018) triggered an increase in primary production as indicated by the increasing chlorophyll a concentrations (Supplementary Figure S5) favoring copiotrophic opportunistic prokaryotes such as Rhodobacterales (ASV4), Flavobacteriales (ASV5, 11 and 140) and Opitutales (ASV9; Wagner-Döbler and Biebl, 2006; Chafee et al., 2018; Jain et al., 2020). In summer, 50 ASVs were enriched in the AW as compared to MS and 68 ASVs in the MS over AW. AW-enriched ASVs had a mean relative abundance < 1% during the summer, however, together accounting for 7% of the total AW community. All MS-enriched ASVs together accounted for 31% of the total MS-associated community in summer. Noticeably, ASVs of the orders Vibrionales, Pirellulales, and Caulobacterales (Figure 2B) were prominently enriched in MS, in agreement with previous findings of members of these taxa attached to particles, particularly in coastal systems in the summer (DeLong et al., 1993; Mestre et al., 2017; Duret et al., 2019). The only two enriched ASVs in MS in fall were putatively assigned to Gemmatales and Alteromonadales, together accounting for 3% of the total MS-associated community (Figure 2C). Both taxa belong to previously reported particle-associated classes (Planctomycetacia and Gammaproteobacteria; Mestre et al., 2017). The low number of enriched ASVs in fall might be explained by similar reasons as described above for spring. Typically, in fall, water column mixing introduces nutrients from the nutrient-rich bottom waters into the nutrient depleted surface waters, thereby stimulating the microbial communities (Supplementary Figure S5; Ivancic et al., 2018).

FIGURE 2

Figure 2. Amplicon sequence variants (ASVs) with relative abundance > 1% and a log₂ fold change > 1 significantly (p < 0.05) enriched in the ambient water (AW) and in marine snow (MS) in (A) winter, (B) summer, and (C) fall. Bubble size represents the relative abundance and colors indicate the assigned taxonomy at the order level.

Enriched Taxa in Different Seasons

Significant enrichment of ASVs in a specific season as compared to the previous and the following season (e.g., fall vs. winter and winter vs. spring) was tested in the AW and MS communities separately. In total, 964 ASVs were enriched in one season compared to the previous or the following season and 58 out of the 964 ASVs exhibited an average relative abundance > 1% (Supplementary Figure S6). The pattern of ASVs enriched in a specific season compared to the previous or following season was similar in the AW and MS (Mantel test: r = 0.7, p < 0.05). This indicates that seasonal conditions shape the taxonomic distribution of both, free-living and MS-associated communities.

In the following, only ASVs with a relative abundance > 1% and those with significant enrichment in a particular season as compared to the previous and the following season are considered. Hence, only truly season-specific ASVs were selected (Figure 3). In the winter, the free-living community was largely characterized by one abundant Synechococcus ASV3, Marinimicrobia (SAR406 clade) ASV65 and Actinomarinales ASV77 (Figures 3A,D). Temperature, nutrients and light regime have been shown to structure Synechococcus into a large number of ecotypes, potentially including a free-living winter ecotype (Ahlgren and Rocap, 2006; Sohm et al., 2016). Marinimicrobia are ubiquitous and typically found in high abundance in the mesopelagic but also in the euphotic zone in the Adriatic Sea in winter (Yilmaz et al., 2016; Mucko et al., 2018). Actinomarinales persist throughout the year in the open Mediterranean Sea and are found below the season pycnocline during the summer (Haro-Moreno et al., 2018). Hence, the high abundance and specific occurrence of Marinimicrobia and Actinomarinales in the winter in the shallow waters of the northern Adriatic Sea are presumably due to water column mixing (Supplementary Figure S5) and lateral advection (Haro-Moreno et al., 2018).

FIGURE 3

Figure 3. Amplicon sequence variants (ASVs) with a relative abundance > 1% and a log₂ fold change > 1 significantly (p < 0.05) enriched when compared to both, the previous and the following seasons in (A–D) the ambient water (AW) and (E–H) in marine snow (MS). Bubble size represents the relative abundance and colors indicate the assigned taxonomy at the order level.

No ASVs were specific for the winter in the MS-associated community (Figures 3E,H). However, compared to the fall, Synechococcus ASV2 and ASV3, both highly abundant, were enriched and dominated the MS-associated and free-living prokaryotic community in winter (Supplementary Figure S6). The high relative abundance of these two Synechococcus ASVs (Supplementary Figure S6) in the AW as well as in MS in winter might be due to a variety of reasons: (i) inefficient separation between MS-attached and free-living prokaryotes, particularly of abundant taxa due to sampling biases introduced by using syringes when MS is small in size; (ii) possible frequent change between the free-living and MS-associated lifestyle of Synechococcus cells as described for Pseudomonas aeruginosa (Son et al., 2015) via RTX proteins. RTX proteins are related to cell motility in cyanobacteria (Linhartová et al., 2010) and were highly abundant in winter MS (see below). (iii) Integration of copepod fecal pellets containing Synechococcus DNA into MS. (iv) Additionally, Synechococcus could initiate particle formation via release of extracellular polymeric substances (Kaltenböck and Herndl, 1992; Deng et al., 2016).

The SAR86 ASV39, Rhodobacterales ASV4, Opitutales ASV9 and two Flavobacteriales ASV34 and ASV35 were enriched in the AW in spring as compared to winter and summer (Figures 3A,B). Rhodobacterales, Flavobacteriales and Opitutae are flexible and generalist heterotrophs (Choo et al., 2007; Newton et al., 2010; Cohan, 2016; Chen et al., 2020). The SAR86 clade comprises aerobic heterotrophs with the potential of energy acquisition via proteorhodopsin (Dupont et al., 2012). Flavobacteriales ASV35 and Puniceispirillales ASV27 were also enriched in MS in spring as compared to winter and summer (Figures 3E,F). The enrichments of the Flavobacteriales ASV35 specifically in spring and of largely the same prokaryotic groups in the AW and MS in spring as compared to winter or summer (Supplementary Figure S6) indicate the similarity of AW and MS communities in spring. This suggests freshly produced MS originating from the spring phytoplankton bloom or a strong effect of environmental parameters on the prokaryotic communities in spring, as mentioned above. Typically, the increasing temperature and solar radiation in spring promote a phytoplankton bloom in the Mediterranean Sea, and hence opportunistic copiotroph prokaryotes such as Flavobacteria (Chafee et al., 2018) experience favorable conditions for growth.

In the summer, Flavobacteriales ASV71 and ASV17 were enriched in the AW community as compared to spring and fall (Figures 3B,C). Rhodobacterales ASV25 was enriched in the AW and MS-associated community (Figure 3). Vibrionales ASV24, ASV7 and Pirellulales ASV18 were enriched in MS and were previously reported as primarily MS-associated prokaryotes (Smith et al., 2013; Salazar et al., 2015; Steiner et al., 2019). The seasonal re-occurrence of Vibrio in the Adriatic and Mediterranean Sea has been linked to water stratification, elevated dissolved organic carbon concentrations and production of transparent exopolymeric particles (Zaccone et al., 2002; Tinta et al., 2015).

Similar to spring, the fall was characterized by the enrichment of the same prokaryotic groups in both the AW and MS in contrast to prokaryotic community composition in summer and winter. Actinomarinales ASV102 (Figure 3) and the abundant prokaryotic groups of the orders Rhodobacterales, Flavobacteriales, SAR11 and Pirellulales were enriched in both, the free-living and MS-associated prokaryotic community in fall (Supplementary Figure S6). Fall-specific free-living ASVs were Pirellulales ASV231 and Thiomicrospirales ASV142 (Figures 3C,D). Thiomicrospirales has been reported from oxygen depleted coastal sites and deep oxygen minimum zones linked to nitrate and the sulfur cycle (Aldunate et al., 2018; Muck et al., 2019). Pirellulales is a heterotrophic order within the class Planctomycetacia (Dedysh et al., 2020). MS-associated ASVs in fall were Planctomycetales ASV314, Cellvibrionales ASV179, Pirellulales ASV265 and Flavobacteriales ASV140 (Figures 3G,H). Planctomycetales produce secondary metabolites such as bacteriocins and ectoines (Wiegand et al., 2020) and are typically found in MS (Mestre et al., 2017; Duret et al., 2019). Cellvibrionales consist of oligotrophs as well as copiotrophs and have been reported from coastal systems and distinct nutrient-rich niches (Spring et al., 2015).

Overall, enrichment analysis revealed that ‘transition seasons,’ i.e., spring and fall, exhibited a high abundance of enriched ASVs compared to winter and summer (1,182 in total) and low abundance of habitat-specific ASVs (two in total). Whereas ‘peak seasons,’ i.e., winter and summer, showed the opposite pattern, with a comparably low number of enriched ASVs compared to spring and fall (827 in total) and a large number of habitat-specific ASVs (728 in total). High habitat specificity in AW vs. MS particularly in winter might be caused by large difference between the two habitats in terms of nutrient and possibly, oxygen availability (e.g., low oxygen microzones in MS). In the summer, MS is generally larger in size than in winter due to low turbulence conditions in the stratified water column facilitating aggregation of colloidal organic matter (Kaltenböck and Herndl, 1992; MacIntyre et al., 1996). Differences in the substrate availability between AW and MS are higher in summer and winter than in spring and fall, when phytoplankton blooms are providing labile substrate to the AW prokaryotic community (Aubry et al., 2004; Bernardi Aubry et al., 2006).

Metagenome Analysis

Taxonomic Analysis

The majority of reads from the metagenomes were assigned to the domain Bacteria, ranging from 43 to 56% (Supplementary Table S3), followed by ‘not assigned’ reads (on average 41 ± 4% for all samples). Reads mapping to archaea were more abundant in winter, both in AW (3–4%) and MS (2%; Supplementary Table S3), in accordance with previous reports on seasonality and lifestyles of archaea in surface waters (Galand et al., 2010; Smith et al., 2013; Haro-Moreno et al., 2017). Eukaryotes contributed substantially less with a maximum contribution of 7–13% in MS in summer (Supplementary Table S3) in the coastal Adriatic than previously reported for MS-associated microbes in coastal (up to 44%) and deep waters (27–40% of all small subunit rRNA genes). Most reads assigned to viruses were detected in the AW in winter (3–4%) and in MS in summer (2%; Supplementary Table S3). Despite reports of lower importance of virus induced mortality as compared to grazing by protists in the Mediterranean Sea, viruses can drastically alter prokaryotic community composition, particularly in eutrophic coastal areas (Siokou-Frangou et al., 2010).

Free-Living and MS-Associated Eukaryotic Community Composition

The largest fraction of eukaryotic reads retrieved from the metagenomes was assigned to copepods (on average 49% ± 27%, Supplementary Table S4 and Figure 4), with Acartia clausii (on average 15 ± 15%) as the most abundant eukaryote. Hydrozoa comprised a substantial fraction of assigned reads in AW in the summer (35 ± 6%). Certain Hydrozoa reproduce during the summer period, reaching high abundances in the Mediterranean Sea (Boero and Fresi, 1986; Di Camillo et al., 2012), suggesting that Hydrozoa sequences in summer originated from free-living hydrozoan larvae. Sequences assigned to Polychaeta (on average 24%) characterized the winter MS community and Platyhelminthes the summer MS community (on average 12%, Supplementary Table S4 and Figure 4). These findings support the notion that MS represents an important microenvironment for eukaryotes, particularly larvae of Polychaeta and juvenile turbellarians (Platyhelminthes; Bochdansky and Herndl, 1992a; Bochdansky et al., 2017). Larvae use MS as a food source and to enhance dispersal since a considerable fraction of MS is only slowly sinking or neutrally buoyant hence being laterally transported over large distances (Bochdansky and Herndl, 1992b). Furthermore, a high number of ‘not assigned’ metazoa (NA_Metazoa), Gastropoda, Bivalvia, and Tunicata were also present in the MS eukaryotic community (Supplementary Table S4 and Figure 4). The high diversity and larger proportion (89 ± 6%) of metazoan reads in MS compared to AW (73 ± 12%; Supplementary Tables S2, S3) suggest an important contribution of metazoans to the remineralization of marine snow in coastal systems (Bochdansky and Herndl, 1992b). Conversely, protists were relatively more abundant in the AW (27 ± 12%) than in MS (11 ± 6%). AW protists were mostly assigned to Alveolata (including dinoflagellates), Rhizaria (including radiolaria), Stramenopiles (including diatoms), and Archaeplastida (including red and green algae), all taxa commonly found in the northern Adriatic Sea (Monti et al., 2012; Piredda et al., 2016; Steiner et al., 2019).

FIGURE 4

Figure 4. Averaged relative abundance of reads assigned to eukaryotes. The fraction of protists has color filling and the fraction of metazoa is shown as white bars with colored shading. Eukaryotic groups with lower relative abundance than 1% are grouped under “other protists” and “other metazoan.”

Functional Assignment of Prokaryotic Communities

The prokaryotic community composition revealed by 16S rRNA amplicon sequences was similar to the prokaryotic community composition obtained by the metagenomes (Mantel test: r = 0.55, p < 0.05, n = 8, Supplementary Figure S3). This allows linking the taxonomic analysis to functional analysis based on metagenomes. The similarity was tested at the taxonomic level ‘order.’ In total, 146 ± 91 orders were detected in each sample. Thereof, 22–59 orders were common in the amplicon and metagenome data in each sample covering the main fraction of sequences (84–96%). This indicates that the similarity between the communities is determined by few abundant orders. Diversity patterns were similar as in the amplicon sequence analysis, however, less pronounced. Species richness in MS (252 ± 10.5) and Shannon diversity in MS in winter (2.7) were lower than in summer and AW samples. The 4,181 identified COGs comprised on average 11 ± 4.4% of all reads with functional assignment.

Enriched COGs in the Free-Living as Compared to MS-Associated Community

The AW community was defined by 189 significantly enriched COGs as compared to the MS-associated community, accounting for 5.2% of all COGs assigned to AW prokaryotes. A particularly high relative abundance (4%) was observed for COGs involved in metabolism, indicating that metabolic capabilities led to a differentiation of the AW from the MS-associated communities. The most abundant COGs (> 0.1%) of the supergroup ‘metabolism’ were related to ‘C – Energy production and conversion,’ including enzymes participating in terpenoid backbone biosynthesis (COG1304), precursors of carotenoids, the building blocks of the proteorhodopsin chromophore retinal (Béja et al., 2001; DeLong and Beja, 2010; Figure 5A). The latter suggests that photoheterotrophy is a survival strategy of the AW community, related to oligotrophy and carbon limitation (Moran and Miller, 2007). The vast majority of reads (65%) aligned to COGs and enriched in the AW community over MS were assigned to Candidatus Pelagibacter sp. HTCC7211, the alpha proteobacterium HIMB114 and Candidatus Puniceispirillum marinum, all abundant proteorhodopsin-containing Alphaproteobacteria (Giovannoni et al., 2005; Lee et al., 2019; Supplementary Table S5). Other genes enriched in AW communities encoded for oxidoreductases and transferases involved in the acquisition of energy from organic compounds via the citric acid cycle (COG1894, COG0074, and COG0045; Figure 5A). An acyl-CoA transferase (COG1804) involved in the catabolism of compatible solutes, such as carnitine, a betaine structurally similar to dimethylsulfoniopropionate (DMSP; Elssner et al., 2001; Todd et al., 2007; Curson et al., 2011) was abundant as well (Figure 5A). Compatible solutes are either generated by or imported into bacterial cells. Three of the most abundant AW specific COGs (category ‘E - Amino acid transport and metabolism’) were components of the ATP-binding cassette (ABC)-type proline/glycine betaine transport system (COG4176, COG2113 and COG4175), indicating potential utilization of compatible solutes as substrate or adaptations to salt stress (Figure 5A). A phytanoyl-CoA dioxygenase (COG5285) involved in the biosynthesis of polyketides (mitomycin antibiotics/fumonisin) indicated the production of antibiotics by the AW prokaryotes, possibly as a survival strategy (Figure 5A; Karimi et al., 2019).

FIGURE 5

Figure 5. Significantly (p < 0.05) enriched clusters of orthologous groups (COGs) in (A) the ambient water (AW) versus marine snow (MS), (B) in winter (WAW) versus summer (SAW) ambient water communities; and (C) in winter (WMS) versus summer (SMS) marine snow. Only COGs with a log₂ fold change > 1 and a relative abundance > 0.1% (indicated by bubble size) are shown. COG numbers and corresponding COG categories are indicated: E – Amino acid transport and metabolism, H – Coenzyme transport and metabolism, P – Inorganic ion transport and metabolism, C – Energy production and conversion, I – Lipid transport and metabolism, G – Carbohydrate transport and metabolism, F – Nucleotide transport and metabolism and Q – Secondary metabolites biosynthesis, transport and catabolism.

Enriched COGs in the MS-Associated as Compared to the Free-Living Community

The MS-associated community was characterized by 219 enriched COGs, accounting for 2.9% of all COGs assigned to MS-associated prokaryotes. Half of the MS-enriched COGs (1.5%) were assigned to functions involved in metabolism. Particularly, a choline-glycine betaine transporter (COG1292) and a Ca²⁺ binding protein RTX toxin-related domain (COG2931) were abundant (Figure 5A). RTX proteins are involved in cell defense acting as cytotoxins in Vibrio or contribute to motility in Cyanobacteria (Linhartová et al., 2010). Furthermore, RTX proteins may form bacterial surface layers as protection and have been observed in pathogenic and non-pathogenic heterotrophic bacteria and cyanobacteria (Linhartová et al., 2010). Most reads (55%) aligned to COGs enriched in MS were assigned to diverse Synechococcus species (CC9311, CC9902, BL107, and WH 7805) and Vibrio caribbeanicus (Supplementary Table S6).

Seasonally Enriched COGs in the Free-Living Community

The functional potential of the AW prokaryotic community in the summer as compared to winter was characterized by 396 enriched COGs constituting 4% of the total summer COGs of the AW prokaryotes, of which approximately half (2.3%) were related to metabolism. The most abundant gene encoded the malic enzyme (COG0281), important in various metabolic pathways including the citric acid cycle (Figure 5B; Takahashi-Íñiguez et al., 2016). The main fraction of all reads aligned to enriched COGs in the summer was assigned to Pelagibacterales (53%) and unclassified Alphaproteobacteria (24%; Supplementary Table S7).

In the AW prokaryotes in winter, 531 COGs (accounting for 8.2% of COGs of the AW community) were enriched compared to the summer, and were mainly associated to Synechococcales (59%) and Pelagibacterales (24%; Supplementary Table S8). More than half (5.9%) of all winter-specific COGs were related to metabolism. Specific functions were assigned to oxidative phosphorylation (COG1009), the biosynthesis of amino acids and other secondary metabolites (COG0147), a cobalamin biosynthesis protein (COG1429) involved in porphyrin and chlorophyll metabolism, an exopolyphosphatase (COG0248), a polyphosphate kinase (COG0855) and a phosphate/phosphonate transport system (COG3639) indicating phosphorus stress (Figure 5B).

Seasonally Enriched COGs in the MS-Associated Community

The MS-associated community showed pronounced differences in the relative abundance of enriched COGs when comparing summer and winter. Enriched COGs in summer represented 2% and in winter 47.3%. However, a similar number of enriched COGs (853 in summer and 806 in winter) were found in both seasons. The reads aligned to summer-specific COGs were mainly assigned to Vibrionales (28%), Pelagibacterales (22%) and Rhodobacterales (13%), and mostly related to metabolism (1.2% of all enriched summer MS COGs, Supplementary Table S9). However, no single COG contributed more than 0.1% to the total COGs in MS in the summer. In winter, almost all enriched COGs were assigned to Synechococcales (92%), reflecting the main difference in the taxonomic assignment of COGs between summer and winter MS-associated communities (compare Supplementary Tables S9, S10). The most abundant COGs (> 0.1%) of the supergroup ‘metabolism’ (57 COGs) spanned eight COG categories (Figure 5C). Of particularly high relative abundance were COGs involved in nitrogen assimilation and storage (COG0436, COG0070 and COG0458), pointing to an adaptation to changing nitrogen availability (Eisenberg et al., 2000; Zhang et al., 2018). Fatty acid biosynthesis genes (COG0304 and COG0439) were enriched in the winter MS-associated community, particularly in Synechococcales members (Supplementary Table S11). Unsaturated membrane lipids in cyanobacteria protect the photosynthetic machinery from photoinhibition associated with low temperatures (Nishida and Murata, 1996). A high relative abundance of enzymes involved in photosynthesis including carotenoid biosynthesis (COG1233), terpenoid backbone biosynthesis (COG1154 and COG0142) and porphyrin and chlorophyll metabolism (COG1429 and COG0155) underlines the impact of Synechococcales on the functional potential of the MS-associated community in the winter (Figure 5C).

Functions of Prokaryotic Communities – Transporter Analysis

Transporters provide insights into the dynamics of nutrient and organic matter uptake of microbial communities (Bergauer et al., 2018). COGs associated with substrate transport comprised 6.6% of the MS- and AW-COGs combined. ABC transporters were the majority (44.4 ± 0.6%) of transporters, in agreement with results from diverse open ocean regions (Bergauer et al., 2018). Highest dissimilarity (determined via SIMPER analysis) between summer MS, winter MS, summer AW and winter AW was reflected in COGs (ENOG4111IMY and ENOG4111FRH) involved in transport of the compatible solute ectoine in winter AW (Supplementary Table S12), which acts as an osmoprotectant and is synthesized in conditions of salinity or temperature stress (Kuhlmann et al., 2011). Another ABC transporter (COG3638) found in high abundance particularly in winter in both AW and MS is associated to phosphonate/phosphite/phosphate transport and linked to phosphorus limitation (Gebhard and Cook, 2008; Supplementary Table S12). Many ABC transporters responsible for the dissimilarity between summer and winter were associated to iron transport and mostly present in summer in both AW and MS (COG4181 and COG3127) and in MS in the summer (COG4594; Supplementary Table S12). This includes high affinity Fe³⁺-citrate siderophores and uptake systems to pirate siderophores of other microbes (Gaballa and Helmann, 2007; Banerjee et al., 2016), pointing to iron limitation and high competition for iron in MS. Other ABC transporters characteristic for the MS-associated community in summer were important in the assembly of cytochrome bd terminal oxidase which is responsible for micro-aerophilic respiration and is overexpressed under nitric oxide stress (Jünemann, 1997; Mason et al., 2009).

The second most abundant transporter subfamily was the major facilitator superfamily (MFS) accounting for 21 ± 0.2% averaged over all seasons and habitats (AW and MS). The MFS is the largest family of secondary transport carriers (Saier et al., 1999; Reddy et al., 2012). Unfortunately, due to poor annotation a detailed functional description was not possible. However, a strong seasonal separation was observed, with a range of COGs seasonally enriched in summer or winter independent of the habitat, and specific COGs linked to summer MS or to summer AW (Supplementary Table S13). Other abundant superfamilies were the tripartite ATP-independent periplasmic transporter (TRAP-T) family comprising 5.7 ± 0.8% of transporter COGs averaged over all seasons and habitats, and the resistance-nodulation-cell division (RND) superfamily with 4.2 ± 0.3% of transporter COGs, also averaged over all seasons and habitats. The presence of TRAP transporters implicates preferential utilization of organic compounds and an energy saving lifestyle, which is advantageous under oligotrophic conditions (Kelly and Thomas, 2001; Bergauer et al., 2018). In this study, the few TRAP transporters showed low dissimilarity between samples and were mainly associated to a high affinity transport system encoded by the dct locus (Forward et al., 1997; Supplementary Table S14). This system transports the C4-dicarboxylates malate, succinate, and fumarate (Forward et al., 1997). The RND family includes transporters for multidrug efflux in Gram-negative bacteria (Nikaido and Takatsuka, 2009). Almost all COGs identified as members of the RND family were assigned to the membrane fusion protein HlyD (Supplementary Tables S15, S16), a specific transporter for the RTX hemolytic toxin HlyA (Pimenta et al., 2005). HlyA is a pore-forming toxin released by many pathogens (including Vibrio spp.) into the medium or directly into the host cell (Costa et al., 2015). However, it has been suggested that HlyD is involved in sync1217 secretion, a RTX toxin domain protein that creates an outer layer membrane barrier to toxins in Synechococcus (Stuart et al., 2013). In this study, 66% of all reads of all seasons and habitats combined targeting COGs of RND family transporters were aligned to COGs of the HlyD protein of Synechococcus (Supplementary Table S16). The large variety of HlyD assigned COGs and the heterogenous distribution through seasons and habitats indicate an important role of defense mechanisms for prokaryotic communities in this coastal system.

Functions of Prokaryotic Communities – Marker Gene Analysis

Functional prokaryotic community analysis based on 69 selected key marker genes revealed the presence of 39 KOs corresponding to diverse genes involved in carbon, nitrogen, sulfur, methane, and carbon monoxide metabolism. Genes involved in H₂ oxidation were not present in the dataset. Interestingly, AW prokaryotes were characterized by a high relative abundance of anaerobic methanogens (mttB, Figure 6) which has been typically described for marine sediments and linked to archaea (Ferry and Lessner, 2008; Sun et al., 2019). However, anaerobic methanogens have also been found in oxygenated coastal marine environments and have been hypothesized to inhabit suboxic and anoxic microenvironments within organic particles such as MS (Alldredge and Cohen, 1987; Ditchfield et al., 2012). MS collected from coastal waters as well as laboratory made artificial MS showed small patches (100 to 300 μm diameter) of reduced oxygen concentrations in marine snow (Shanks and Reeder, 1993). These microhabitats were subsequently suggested as sites of anaerobic processes (DeLong, 1992; Shanks and Reeder, 1993; Bianchi et al., 2018) including methanogenic archaea (Ditchfield et al., 2012). Our results indicate a larger potential for methanogenesis in the AW than in MS. The AW also includes small particles, which might harbor methanogens. Another possible explanation could be frequent attachment and detachment from MS by methanogenic prokaryotes, being active in microzones of MS and detaching as MS becomes dispersed (Son et al., 2015). The aprA and aprB genes were very abundant in the AW, particularly in the winter (Figure 6). The aprA and aprB genes are involved in anaerobic dissimilatory sulfate reduction by sulfate reducing organisms (SRO) such as Syntrophobacterales, Thermodesulfobacterium, Thermodesulfovibrio, Archaeoglobus, and some deltaproteobacterial lineages (Meyer and Kuever, 2007b). SRO have been found to be frequently associated with methanogens (Grein et al., 2013) which were abundant in winter AW as mentioned above. However, apr genes are also found in anoxygenic phototrophic and chemolithoautotrophic sulfur oxidizing bacteria (SOB; Frigaard and Dahl, 2008), as well as chemolithoheterotrophic SOB (Meyer and Kuever, 2007a), including the abundant Pelagibacterales (Giovannoni et al., 2005) and thus, are widespread throughout the oxygenated oceanic water column (De Corte et al., submitted).

FIGURE 6

Figure 6. Fragments per kilo base million (FPKM) normalized abundance of the 39 KOs (Kyoto Encyclopedia of Genes and Genomes Orthology) identified and corresponding to genes involved in the carbon, nitrogen, sulfur, methane, and carbon monoxide (CO) metabolism. Gene names are listed and described in Supplementary Table S1. WAW1 and WAW2 indicate winter ambient water duplicates, SAW1 and SAW2 indicate summer ambient water duplicates, WMS1 and WMS2 correspond to winter marine snow duplicates, SMS1 and SMS2 indicate summer marine snow duplicates.

Anaerobic metabolisms such as dissimilatory nitrate reduction (nrfA, napA, and napB) were only present in summer MS (Figure 6), indicating that large marine snow present in summer (as compared to winter as observed by SCUBA divers) was colonized by prokaryotes specialized in attachment, biofilm formation and anaerobic nitrate respiration. Most bacterial lineages can perform dissimilatory nitrite reduction to ammonium (Kuypers et al., 2018). Other genes involved in denitrification and dissimilatory nitrate reduction (nirK, narI, and narV), the assimilatory nitrate reduction gene narB and genes encoding nitrate reductases – nitrite oxidoreductases (narH, narY, nxrB, narG, narZ, and nxrA) and methane/ammonia monooxygenases (pmoA-amoA) were generally more frequent in winter than in summer in both habitats (Figure 6), although in low abundances. Nitrite and nitrate reductases are widespread among bacteria and archaea and occur in low oxygen and anoxic environments where nitrite and nitrate is present (Kuypers et al., 2018). NirK is commonly used as a marker for denitrification, however, it is not exclusively found in denitrifiers (Kuypers et al., 2018). Ammonia monooxygenases occur in beta- and gammaproteobacterial classes including Nitrosomonas, Nitrosococcus, and Nitrospirae (Chain et al., 2003; Klotz et al., 2006; Daims et al., 2015) and the archaeal phylum Thaumarchaeota (Brochier-Armanet et al., 2012). None of them occurred in large abundances in this study (Supplementary Figure S3). Summer AW prokaryotes exhibited the highest relative abundance of genes involved in CO-oxidation (coxL, coxS, coxM, and cutM; Figure 6). These genes occur mostly in Roseobacter to supplement heterotrophic growth with inorganic carbon (Moran and Miller, 2007). In this study, however, most reads (46%) assigned to CO-oxidation COGs (COG1529, COG2080, COG1319) affiliated to the Alphaproteobacterium HIMB114 and Candidatus Puniceispirillum marinum (Supplementary Table S17) indicating supplementing heterotrophic growth with energy derived from carbon monoxide (Moran and Miller, 2007; Oh et al., 2010). Winter MS exhibited the highest abundance of genes from the Calvin cycle (rbcS, rbcL, and prkB; Figure 6) in accordance with the high abundance of cyanobacteria in MS (Supplementary Figure S3).

Summary and Conclusion

We analyzed the community composition and functional potential of the AW and MS-associated communities over a seasonal cycle. The winter and summer were characterized by habitat specific prokaryotes, while the ‘transition seasons,’ i.e., spring and fall were defined by season-specific prokaryotes, similar in the AW and MS-associated communities.

Differences in the nutrient availability between AW and MS are larger in winter and summer than in spring and fall favoring pronounced differences in the microbial community composition in the winter and summer between AW and MS. A selective advantage of prokaryotes associated to MS possibly includes the utilization of RTX toxin-related proteins for protection and motility. Adaptations to low temperature in winter might involve protection of the photosynthetic machinery from photoinhibition by low temperature. Synechococcus was not significantly enriched in MS, however, most of the enriched potential functions were assigned to Synechococcus in MS in the winter. The prominent appearance of Synechococcus associated to MS suggests that MS represents a hotspot for autotrophy in winter. In contrast, the development of low oxygen microzones within MS during water column stratification and low turbulence conditions in summer offers suitable conditions for anaerobic metabolism such as nitrite respiration.

Adaptations to a free-living life-style in winter included coping with phosphorus limitation and energy acquisition from labile organic matter such as compatible solutes. In the summer, prokaryotes supplement heterotrophic growth with solar radiation energy harvested with proteorhodopsins (Gómez-Consarnau et al., 2019) or with CO-oxidation (Moran and Miller, 2007).

The ‘transition seasons’ spring and fall were largely characterized by similar composition and functional potential of the AW and MS-associated prokaryotic communities in contrast to winter and summer. The fact that the same ASVs were enriched in the AW and MS indicates that environmental parameters of the ‘transition seasons’ influence both, the AW and MS-associated communities in the same way leading to an overall convergence of the microbial community composition and functional potential between AW and MS during these ‘transition seasons.’

Overall, a pronounced seasonality in both, the AW and the MS-associated communities was detectable. The MS-associated microbial community, however, was less sensitive to seasonal changes than the AW microbial community. Thus, the MS-associated community was more stable over the seasonal cycle than the AW microbial community. The large diversity and proportion of metazoans associated with MS might substantially affect the substrate composition of marine snow and thus, prokaryotic community composition and function, particularly in the summer. In contrast, the relatively low abundance of protists (including eukaryotic phytoplankton) together with the high abundance of Synechococcus in winter MS suggests a major prokaryotic contribution to marine snow in this season.

Data Availability Statement

The datasets generated in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ebi.ac.uk/ena, PRJEB38662.

Author Contributions

PS, JG, TR, ES, and GH conceived the research. PS conducted DNA extraction and prepared 16SrRNA and metagenome sequencing. JG, AO, and PS analyzed the data under supervision and following the advices of TR and EF. PS wrote the manuscript. All authors revised the draft version and approved the final version.

Funding

This study was funded by the Wittgenstein prize (Austrian Science Fund, FWF, Z194-B1) to GH and the projects ARTEMIS (P28781-B21) and DK microbial nitrogen cycling (W1257-B20) to GH. This work was in partial fulfillment of the requirements for a Ph.D. degree from the University of Vienna to PS.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We thank the research group from the Ruder Bošković Institute at Rovinj for providing environmental data, their help and great support in organizing field trips and lab work. Special thanks go to Marino Korlević, Paolo Paliaga, Ingrid Ivanèić, Mirjana Najdek, and the captain of the R/V Burin. Additionally, we thank Elisabeth L. Clifford, Daniele De Corte, and Maria Pinto for support in sampling, in the lab and in bioinformatics analysis. We also thank Miguel Guerreiro for support on bioinformatics analysis and Abhishek Srivastava for support on microbial pathway analysis. Part of the analysis of this study was done at the Institut de Ciències del Mar (CSIC) in Barcelona during a secondment. PS particularly thanks Josep M. Gasol for help and care for any issue during the secondment.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fmicb.2020.584222/full#supplementary-material

Supplementary Figure 1 | Rarefaction curves of amplicon sequences normalized to the sample with smallest number of sequences (08 February 2016 marine snow).

Supplementary Figure 2 | Shannon diversity index (A) and species richness (number of taxa; B) of rarified prokaryotic communities in different sampling days based 16S rRNA amplicon sequence variant (ASV) analysis. Symbols indicate ambient water (AW) and marine snow (MS) habitats.

Supplementary Figure 3 | Contribution of the 30 most abundant (A) MS-associated (MS) and (B) ambient water (AW) prokaryotic orders to the communities based on variance-stabilized 16S rRNA amplicon sequences and the 20 most abundant (C) MS-associated (MS) and (D) AW prokaryotic orders based on variance stabilized 16S rRNA sequences extracted from the metagenomes. Less abundant prokaryotic orders are grouped under “Others.” Samples from April to November were taken in 2015 and samples from February to July were taken in 2016.

Supplementary Figure 4 | Venn diagrams showing the number of prokaryotic classes occurring in the ambient water (AW; blue), the marine snow (MS; yellow) and in both (overlapping circles) each season.

Supplementary Figure 5 | Temperature, salinity and chlorophyll a profiles at station RV001 in (A) spring (14/04/15 black bar and 06/04/16 grey bar), (B) summer (26/06/15 black bar and 19/07/16 gray bar), (C) fall (26/10/15 black bar and 09/11/15 gray bar) and (D) winter (20/01/16).

Supplementary Figure 6 | Amplicon sequence variants (ASVs) with a relative abundance > 1% and a log₂ fold change > 1 significantly (p < 0.05) enriched in one of two consecutive seasons in (A) the ambient water and (B) in marine snow (MS). Bubble size represents the relative abundance and colors indicate the assigned taxonomy at the order level.

Supplementary Table 1 | Date, sample type and diversity indices of samples rarified to the lowest number of amplicon sequence variants (ASVs).

Supplementary Table 2 | SIMPER analysis based on relative abundance of reads assigned to eukaryotes. WAW, winter ambient water; SAW, summer ambient water; WMS, winter marine snow; SMS, summer marine snow.

Supplementary Table 3 | Abundance and taxonomic assignment of trimmed reads, and diversity indexes of the metagenomes.

Supplementary Table 4 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in the ambient water (AW).

Supplementary Table 5 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in marine snow (MS).

Supplementary Table 6 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in summer ambient water (AW).

Supplementary Table 7 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in winter ambient water (AW).

Supplementary Table 8 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in summer marine snow (MS).

Supplementary Table 9 | Taxonomy, number of reads and relative abundance of reads aligned to significantly enriched species in winter marine snow (MS).

Supplementary Table 10 | Taxonomy, number of reads and relative abundance of reads aligned to COG0304 and COG0439 in winter marine snow (MS).

Supplementary Table 11 | SIMPER analysis based on variance stabilized abundance of reads aligned to ATP-binding cassette (ABC) transporter clusters of orthologous groups (COGs). WAW, winter ambient water; SAW, summer ambient water; WMS, winter marine snow; SMS, summer marine snow.

Supplementary Table 12 | SIMPER analysis based on variance stabilized abundance of reads aligned to major facilitator superfamily (MFS) transporter clusters of orthologous groups (COGs). WAW, winter ambient water; SAW, summer ambient water; WMS, winter marine snow; SMS, summer marine snow.

Supplementary Table 13 | SIMPER analysis based on variance stabilized abundance of reads aligned to tripartite ATP-independent periplasmic transporter (TRAP-T) family clusters of orthologous groups (COGs). WAW, winter ambient water; SAW, summer ambient water; WMS, winter marine snow; SMS, summer marine snow.

Supplementary Table 14 | SIMPER analysis based on variance stabilized abundance of reads aligned to resistance-nodulation-cell division (RND) superfamily clusters of orthologous groups (COGs). WAW, winter ambient water; SAW, summer ambient water; WMS, winter marine snow; SMS, summer marine snow.

Supplementary Table 15 | List of clusters of orthologous groups (COGs) associated with the hylD gene and a table showing the taxonomic affiliation (Taxa), average (AV), sum and relative abundance (%) of reads aligned to clusters of orthologous groups (COGs) associated with the gene hylD of all samples.

Supplementary Table 16 | Taxonomy (Taxa), number of reads (Reads) and relative abundance (%) of reads aligned to COGs associated to CO-oxidation (COG1529, COG2080, COG1319) in winter MS.

Supplementary Table 17 | Marker genes used normalized abundance. The presence (+) in at least one sample of the dataset is indicated. WAW1 and WAW2, winter ambient water duplicates; SAW1 and SAW2, summer ambient water duplicates; WMS1 and WMS2, winter marine snow duplicates; SMS1 and SMS2, summer marine snow duplicates.

Footnotes

References

Acinas, S. G., Antón, J., and Rodríguez-Valera, F. (1999). Diversity of free-living and attached bacteria in offshore western Mediterranean waters as depicted by analysis of genes encoding 16S rRNA. Appl. Environ. Microbiol. 65, 514–522. doi: 10.1128/aem.65.2.514-522.1999

PubMed Abstract | CrossRef Full Text | Google Scholar

Acinas, S. G., Sánchez, P., Salazar, G., Cornejo-Castillo, F. M., Sebastián, M., Logares, R., et al. (2019). Metabolic architecture of the deep ocean microbiome. bioRxiv [Preprint] doi: 10.1101/635680