Increased contribution of parasites in microbial eukaryotic communities of different Aegean Sea coastal systems

Collection of samples and DNA extraction

Samples were collected from seven gulfs in the Aegean Sea, Greece, set within a polygonal area of 72,600 km² (Table 1; Fig. S1), characterized by a range of environmental conditions due to differences in hydrology, geomorphology, substrate, terrestrial runoff and local anthropogenic pressures (Spatharis et al., 2019). Within each site, five stations were sampled at 1 m (surface), at the Secchi disc depth (Secchi) and at the two times Secchi disc depth (2xSecchi) or maximum depth, depending on site depth (Table S1), collecting 1 l seawater with a Niskin type sampler. In some cases duplicates or triplicates were collected in order to check sampling and sequencing consistency (Table S1).

Sampling took place during July 2014 and was carried out within 19 days (5–24 July) to minimize the effect of temporal variation on assemblage composition. July was selected for sampling as previous studies have shown that physico-chemical variables and phytoplankton composition are relatively stable during this summer period, at least for a subset of the coastal areas included in the present study (Spatharis et al., 2007; Naselli-Flores et al., 2003); this is in contrast with winter months, during which episodic rainfall events, and strong wind mixing of the water column and sediment resuspension—which at this large scale may not simultaneously affect all sites—add noise that could distort food web dynamics. Water samples were filtered using low vacuum filtration (<150 mmHg) on 0.2 µm isopore filters (Sartorius Stedim Biotech, Germany).

Environmental parameters

At each station and depth, several environmental parameters were recorded. Specifically, salinity and temperature were recorded onsite, while 3 l seawater samples were collected with a Niskin type sampler for later nutrient measurement (nitrate, nitrite phosphate, silicate, organic nitrogen, organic phosphorus). Organic nitrogen and phosphorus were strongly correlated with dissolved inorganic nitrogen (DIN) and dissolved inorganic phosphorus (DIP) and were therefore excluded from further analysis. The covariates used in the analysis were thus salinity, temperature, DIN, DIP and silicate.

DNA extraction and sequencing

DNA extraction from filters was performed with the MoBio Power Soil kit (MoBio Inc. Carlsbad, CA, USA) following its standard protocol with minor modifications for filters processing.

Sequencing of the V2–V3 region of the 18S rRNA gene was performed upon amplification using the primer pair 18S-82F (5′-GAAACTGCGAATGGCTC-3′) (López-García et al., 2003) and Euk-516r (5′-ACCAGACTTGCCCTCC-3′) for Eukaryotes (Amann et al., 1990). Construction of libraries was performed by ‘Genes Diffusion’ company (Lille, France) and amplicons were finally sequenced with Illumina MiSeq PE 2x300 (CNRS-UMR8199; Illumina, Lille, France).

18S rRNA gene amplicon analysis

Processing of the resulting sequences, i.e., sequence assembly and quality control, was performed with the MOTHUR software (v 1.35) (Schloss et al., 2009). Only sequences with ≥480 bp, no ambiguous bases and homopolymers shorter than 8 bp were considered for further analysis. These sequences were aligned using the SILVA SSU database (release 119) (Pruesse et al., 2007). Chimeras were removed using the Uchime Software (Edgar et al., 2011). All sequences were binned into operational taxonomic units (OTUs) and were clustered (average neighbour algorithm) at 97% sequence similarity. Single singletons, that appeared only once in the whole dataset, were removed using MOTHUR (v 1.35). Coverage values were calculated with MOTHUR (v 1.35) as well as diversity indices. The batch of sequences from this study has been submitted in NCBI Short Read Archive under accession number PRJNA515026. Taxonomic classification was assigned using BLAST (Altschul et al., 1990) on the Protist Ribosomal Reference (PR2 version 4.14.0) curated Database (built on GenBank; June 2021), containing 197,602 sequences (Guillou et al., 2013).

Similarity profiles, multivariate analysis and differentially abundant features

The Bray–Curtis similarities coefficients were calculated, based on non-transformed number of reads in order to identify relationships between the samples, using R (R Core Team, 2020). NMDS ordination plots were prepared in R (package vegan; Oksanen et al., 2020) using Bray-curtis similarities.

Differentially abundant OTUs between gulfs, were identified with DESeq2 package version 3.0.2 (Anders & Huber, 2010; Love, Huber & Anders, 2014) which performs differential analysis of count data, estimating dispersions and fold changes, thus enabling a more quantitative analysis focused on the strength rather than the mere presence of differential expression.

Abundant and rare OTUs, generalists and specialists and trophic groups

OTUs were classified as abundant or rare based on total relative abundances of on non-transformed number of reads and on the thresholds used in previous studies (Logares et al., 2014) using >1% for abundant OTUs and <0.2% for rare OTUs. Habitat specialization was calculated as described in Székely & Langenheder (2014) using Levins’ niche width index (B) (Levins, 1968), where $\mathrm{B}=1/{\sum }_{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{p}}_{\mathrm{ij}}^2$ (pij: the proportion of OTU j in sample I; N: the total number of samples). Thus B, although is not taking the environmental conditions in a local community into account, is describing the extent of niche specialization based on the distribution of OTUs abundances.

Specialization categories were set along arbitrary cut-off values (B > 15, B:10–15, B:1–10 B = 1) with B = 1 corresponding to extreme specialists (one sample) and B > 15 corresponding to top generalists (>66% of available habitats). Niche width was calculated based on the presence of each OTU at all sites and depths (112 samples) as each sample was considered to be a different habitat in terms of physico-chemical conditions. OTUs were sorted into major trophic groups, such as autotrophs, nanoheterotrophs (picograzers), nanoplankton grazers (nanograzers), microplankton grazers (micrograzers), mixotrophs and parasites following the classification shown in Genitsaris et al. (2016).

Correlation and network analysis

For the network analysis, we focused on the OTUs indicated by DeSeq as they were considered ‘key’ species for community structure regulation. The relationships/interactions amongst these OTUs were characterized by computing the maximal information coefficient (MIC) between each OTU pair (Reshef et al., 2016), calculated using MICtools v1.1.4 (Albanese et al., 2018). MIC captures associations between data and provides a score that represents the strength of the relationship between data pairs. MIC values > 0.5, corresponding to a p-value < 0.05, were used for the networks. For networks visualization Cytoscape 3.0 (Smoot et al., 2011) was used.

Protist community structure

Overall 4,389,394 sequences were obtained resulting in 5,005 operational taxonomic units (OTUs) in all 112 samples. Diversity indices indicated that the highest richness as well as Shannon diversity were observed in LK, while in S the respective values were the lowest ones followed by T and KV (Table S2). Shannon diversity index was stable over depth in most gulfs apart from S, T and KV where diversity increased in deeper layers (Table S2).

Protistan community composition (PCC) similarities within gulfs decreased over depth apart from K, S and KV (Fig. S2 a). Regarding between gulf similarities anova results on hellinger transformed OTUs abundances indicated several significant differences (p < 0.05; Table S3A) at the surface level while at Secchi depth the only difference was observed between T and K (p = 0.041) and at 2xSecchi no differences were observed (Table S3A). Bray-Curtis similarities between gulfs ranged from 0.25 to 0.30 with the highest values being observed in surface (Fig. S2B). This range was very low, although differences were significant between surface and the other depths (t-test Monte-Carlo permutation; p = 0.0001 for both Secchi and 2xSecchi). Overall these results implied more distinct between gulfs, but homogenous within gulfs, communities at the surface level while communities between gulfs were less distinguishable at higher depths.

Similar results were observed after analysis with NMDS showing spatial grouping, mostly for the surface and Secchi level samples (Fig. 1). Ordination stress was quite high (0.12, 0.14 and 0.18 for different depths respectively). The only environmental parameters that were important for these ordinations (p < 0.001) were silicate and ammonia (p = 0.00099) for both surface and secchi samples, NO2 only for surface (p = 0.00099) and phosphates and Chla (p = 0.00099) for Secchi depth. Finally at 2xsecchi depth only temperature appeared as significant for the ordination (Fig. 1).

Taxonomic classification for each gulf and depth revealed different patterns. Overall nine major taxonomic groups comprised ∼90% of taxonomic diversity in all samples (Fig. 2). Bacillariophyta, Dinophyceae and Syndiniales were the most abundant groups in most cases, exhibiting fluctuating abundances across gulfs and depths. However Chrysophyceae prevailed in S samples where Bacillariophyta decreased and Haptophyta exhibited their highest relative abundances in K (Fig. 2).

Effect of depth in trophic groups; abundant OTUs and niche breadth

Overall autotrophic algae (shown in green colors in Fig. 2) abundance dropped from surface to deeper layers. It has to be noted though that in KV, M and LK samples 2xSecchi depth coincides with maximum depth; thus, these samples were collected from bottom water (Table S1). Overall parasites (mainly Syndiniales) relative abundances were constantly >10% in all gulfs (Fig. 2). Correlation analysis and the calculation of MIC scores for each depth separately, showed that in surface strong correlations were observed between all trophic groups apart from autotrophs, while in Secchi depth correlations were observed only between grazers and parasites and no correlations were observed in 2xSecchi (Table S4).

The whole dataset included 25 abundant (>1%) and 60 common (>0.2%) OTUs, accounting for 0.5% and 1.19% of the total OTUs number respectively (Table S5). However, in terms of relative abundances these OTUs accounted for >60% of protistan community composition. These abundant species belonged to all trophic groups, exhibiting different abundance patterns amongst gulfs and depths (i.e., prevalence of Chrysophyceae in S samples) (Fig. 3).

Amongst the 5,005 OTUs detected, 2,314 represented extreme specialists (present only in one site and depth), 2,038 represented specialists (present in <30% of all sites and depths; 1 < B < 10), 227 represented generalists (10 < B < 15) and 426 were top generalists (B > 15) reaching ∼50% of OTUs number/sample (Fig. S3A; Table S5). In most samples top generalists and generalists accounted for >90% of OTUs relative abundances (Fig. S3B), apart from secchi and 2xsecchi samples in KV, S and T and 2xSecchi in G. Overall the number of extreme specialists and specialists increased with depth with differences being significant between surface and 2xsecchi (p = 0.0126; F = 5.06).

The vast majority of abundant OTUs consisted of top generalists (17/25) and generalists (5/25), while a similar pattern was observed for common OTUs (50/60; 6/60) (Table S5). Only 17/426 top-generalists were abundant and 47/426 were common, while the rest were rare species. Similarly for generalists 216/227 species were rare (Table S5).

Differentially abundant OTUs across depths exhibit different trophic and niche breadth patterns

The number of differentially abundant OTUs across the seascape (using the binomial test and false discovery rate < 0.05), varied between depths, with 44, 23 and 51 OTUs at surface, Secchi, and 2XSecchi depths, respectively (Table S6). For surface samples these OTUs belonged to all trophic groups (Table S6) and more specifically included increased relative abundances of phototrophic algae in T, KV and G, of nanograzers in G and S, of mixotrophs in S and of parasites in KV and M (Fig. 4). For secchi depth, differences between sites were attributed to fluctuations of OTUs, belonging to mixotrophic and parasitic groups (Fig. 4) that both peaked in KV followed by S, while fluctuations of nanograzers were also important. Finally, for 2xsecchi, the majority of differences were attributed to OTUs that mostly belonged to heterotrophic and mixotrophic groups, although in specific sites such as LK, M, S and T increases of algae (Chrysophyceae) were responsible for gulf separation, while parasites peaked in KV (Fig. 4, Table S6). We have to highlight that in all depths >70% of differentially abundant OTUs belonged to rare OTUs (Table S6) showing the importance of rare members of the total protistan community for spatial differentiations in PCC.

Regarding these OTUs niche breadth, at surface 28/44 OTUs were top-generalists and generalists, while at Secchi and 2xSecchi depths the respective ratios were only 7/23 and 9/51 (Table S6), corroborating with the total number of specialists increasing with depth (Fig. S3A).

Amongst these generalists and top-generalists differentially abundant species only 4/28, 1/7 and 3/9 species were also abundant at surface, Secchi and 2Xsechhi, respectively (Table S6, Fig. 3). These OTUs were related to Leptocylindrus (OTU00031) and to Chaetoceros (OTU00036) in surface, to Chrysophyceae (OTU00022) in 2X Secchi, while the uncultured Syndiniales (OTU00056) was important for all depths and an uncultured Dinophyceae (OTU00004) was important in surface and 2XSecchi (Table 2). Apart from OTU00056, two more species, OTU00288 and OTU00737 were important for differences between gulfs in all depths. OTU00288 was a generalist although rare mixotrophic haptophyte, and OTU00737 was a rare specialist nanograzer belonging to Spirotrichea.

Correlations between differentially abundant OTUs; network analysis for different depths

Network analysis for the 44 surface species that were responsible for between gulfs differences indicated 155 significant relationships between 42 species while two species indicated no correlation (Table S7). Amongst these species 11 species exhibited 105 significant interactions amongst them as well as with 19 other species and were responsible for the majority of interactions observed (Fig. 5A,; Tables 2; S7).

These eleven species belonged to algae (Chlorophyta, Bacillariophyta) that mostly interacted with parasites followed by other algae and mixotrophs, parasites (Syndiniales) apart from algae also interacted with mixotrophs, and mixotrophs (Haptophyta, Dinophyceae), nanoheterotrophs and nanograzers that interacted almost equally with all groups (Table S6). The highest MIC scores were observed between algae, mixotrophs and some parasites while the correlations with other heterotrophs were minor (Tables 2; S7). Most importantly these 11 OTUs belonged to generalists and top-generalists and only the nanonoheterotroph OTU00263 and the parasitic OTU00177 (Syndiniales) were specialists (Table S6). On the other hand, the Secchi network had much lower centrality (that is few central nodes that are close to all the other nodes) with only 15 species (out of 23) exhibiting 23 relationships with MIC scores >0.05, thus indicating weaker interactions between species. However species exhibiting the higher number of interactions included mainly parasites, nano-grazers and the rare generalist mixotrophic haptophyte (OTU00288; 100% identity to AB058358 Haptolina) that was present in all three networks. Phototrophs (Bacillariophyta) had low centralities interacting mostly between them and with nanograzers (Ciliophora) (Table S6). The top generalist but rare nanograzer OTU00409 (Spirotrichea) strongly correlated with both Syndiniales, and mixotrophs while the parasitic Blastodinium (OTU00052) interacted with mixotrophic Dinophyceae and connected to the rest of the network through OTU00056. Finally, at 2xsecchi two networks occurred with the first including only six species (mainly parasites; Syndiniales and Blastodinium) and the second including only 10 interactions between nine species and being built on the associations between parasites and algae (Chrysophyceae; Bacillariophyta) that were generalists or top generalist but not abundant (Fig. 5C; Tables 2; S7).