Surface waters were collected on 8th Aug 2013 in a 6 km inshore to offshore transect off the coast of Barcelona, NW Mediterranean. Five stations were sampled along the transect (Figure 1). Samples were sieved through a 200-μm mesh and transported to the laboratory within 2 h. Basic physical data was measured in situ with a conductivity, temperature and depth probe and surface salinity was analyzed with an AUTOSAL salinometer. The concentration of inorganic nutrients was determined spectrophotometrically by using an Alliance Evolution II autoanalyzer according to standard procedures (Grasshoff et al., 1983). Chlorophyll a (Chl a) concentration was measured from acetone extracts by fluorometry from the total fraction (< 200 μm) and the fractions less than 20 and 3 μm. To collect microbial biomass, about 5 l of surface seawater was sequentially filtered through a 3- and a 0.2-μm pore-size polycarbonate filters (Poretics, GE Osmotics, Delft, Netherlands) using a peristaltic pump. The filters were stored in cryogenic vials containing 1.7 ml of lysis buffer (50 mM Tris-HCl pH 8.3, 40 mM EDTA pH 8.0 and 0.75 M sucrose) at −80°C until further processing.

The 0.2-μm filters were treated with lysozyme, proteinase K and sodium dodecyl sulfate, and the nucleic acids were extracted with phenol and concentrated in an Amicon 100 (Millipore), as described in Massana et al. (1997). DNA was quantified spectrophotometrically (Nanodrop, Thermo Scientific), and two subsamples from each extraction were sent for sequencing. Sequencing was performed by the Research and Testing Laboratory (Lubbock, TX, USA; http://www.researchandtesting.com/). Primers 341F (5′-CCTACGGGNGGCWGCAG-3′), and 805R (5′-GACTACHVGGGTATCTAATCC-3′) were used to amplify bacterial 16 S rRNA gene (Herlemann et al., 2011) and primers TAReukFWD1 5′-CCAGCASCYGCGGTAATTCC-3′ and TAReukREV3 5′-ACTTTCGTTCTTGATYRA-3′ were used to amplify the V4 region of the eukaryotic 18 S rRNA gene (Stoeck et al., 2010). Pyrosequencing was performed using the bTEFAP method by 454 GL FLX technology as described previously (Dowd et al., 2008). Illumina MiSeq 2 × 250 flow cells were used for Illumina sequencing following protocols described elsewhere (Cúcio et al., 2016). Approximately 30,000 raw sequences per sample were obtained.

High-Performance computing analyses were run at the Marine Bioinformatics Service of the Institut de Ciències del Mar (ICM-CSIC) in Barcelona. Reads from the two sequencing methodologies underwent method-specific quality filtering before being pooled. Bacterial-454 data was filtered by quality using QIIME (Quantitative Insights Into Microbial Ecology, Caporaso et al., 2010) as described in Sánchez et al. (2013). Briefly, sequences from the 454 run were assigned a sample IDs using a mapping file and the barcode assigned to each sample. After sample IDs were assigned, bacterial sequences were removed from the subsequent analyses if they were shorter than 150 bp or longer than 500 bp, had an average quality score < 25 calculated in sliding windows of 50 bp, contained more than two ambiguous characters or had an uncorrectable barcode. Eukaryotic-454 reads were quality checked and demultiplexed with QIIME following the same parameters described in Pernice et al. (2015). Shortly, sequences shorter than 150 bp or longer than 600 bp, with more than three mismatches in the primer, or having homopolymers longer than 8 bp were removed. Phred quality was analyzed in 50 bp running windows. Illumina sequences from bacteria and picoeukaryotes were quality filtered following a custom made pipeline (https://github.com/ramalok). Briefly, BayesHammer error correction of sequence reads was performed with SPAdes software (Nurk et al., 2013). Sequences were assembled with PEAR (http://pear.php.net/) and quality filtered in UPARSE (fastq_maxee value = 1). Clean bacterial-454 and bacterial-Illumina sequences were pooled and processed together; eukaryotic-454 and eukaryotic-Illumina sequences were also pooled together. Since 454 and Illumina sequences may have different length, bacterial sequences were truncated at equal depth (400 bp). However, for picoeukaryotes we did not truncate sequences, since large natural variability in the length of 18 S rRNA from different taxa occur. Sequences of both datasets were clustered into operational taxonomic units (OTUs) at 97% cutoff using the UPARSE algorithm implemented in USEARCH (Edgar, 2013). Both de novo chimera check and by comparison to reference database (SILVA) were done using the UCHIME algorithm (Edgar et al., 2011). Chimeric sequences and singleton OTUs (those represented by a single sequence) were removed. Taxonomic assignment of bacterial OTUs was performed using the BLAST classifier and the version 119 of the SILVA SSURef non-redundant database. OTUs assigned to chloroplasts were removed for subsequent analyses. For picoeukaryotes, OTUs were taxonomically classified by using BLAST against two reference databases: PR² (Guillou et al., 2013) and a marine microeukaryote database (MASS9013, Pernice et al., 2013). After taxonomic assignment, metazoan OTUs were removed. Sequence data has been submitted to the Genbank Sequence Read Archive under accession number SRP079955.

Statistical analyses were performed using the R statistical software (R Developement Core Team, 2015) and the packages vegan, labdsv, venneuler, hmisc, and corrgram. Alpha- and betadiversity analyses were performed using an OTU abundance table that was previously subsampled down to the minimum number of reads in order to avoid artifacts due to an uneven sequencing effort among samples. For alphadiversity analyses, we calculated the Chao1 index as a measure of richness and Shannon and Simpson indices as diversity metrics. Differences in microbial composition (betadiversity) were assessed using hierarchical clustering of Bray-Curtis dissimilarity matrices and the Unweighted Pair Group Method with Arithmetic Mean algorithm (UPGMA). To search for “indicator species” we used the IndVal (INDicator VALues; Dufrêne and Legendre, 1997) analysis, which identifies indicator species based on OTU fidelity and relative abundance. Only OTUs with significant p-values (< 0.05), and >0.3 IndVal values were considered. To assess links between diversity and environmental data we performed linear regressions and pairwise correlations (Pearson's correlation coefficient). The results were thresholded at p < 0.05. Analysis of variance was run to test for differences among diversity data and categories (sequencing method, station) with Tukey-Kramer post hoc comparisons at the 5% significance.

Sequencing of all bacterial and most picoeukaryotic samples was successful yet two picoeukaryotic replicates (1a, 454, and Illumina) resulted in a low number of reads and were discarded from further analyses. The bacterial dataset resulted in 277,212 high quality reads that clustered in a total of 658 OTUs at 97% similarity. From those only 34.7% of OTUs were shared between samples sequenced by either 454 or Illumina (Figure 2). However, the unique OTUs in each methodology correspond to rare members; the proportion of shared OTUs (350 out of 658) represented 99.4% of the reads. We found a good correlation between the relative abundance of each OTU sequenced by both methodologies (R = 0.87, p < 0.001). Likewise, when grouping OTUs into the main bacterial taxa, a good agreement between contributions obtained by 454 or Illumina was found (R = 0.81, p < 0.001). In both cases, most bacterial sequences were related to the phyla Proteobacteria (average of all bacterial dataset, 72%), Bacteroidetes (20%), and Cyanobacteria (5%). Within the Proteobacteria, the most prevalent classes were the Alpha-(50%) and the Gammaproteobacteria (20%), whereas the Beta-, Delta-, and Epsilon- were present at low relative abundances (grouped as “Other Proteobacteria,” Figure 3). Within the Alphaproteobacteria, the OTUs showing higher relative abundances were affiliated to the Rhodobacterales, Rhodospirillales, Rickettsiales, and the SAR11 clade. The Bacteroidetes were largely represented by members of the Flavobacteriia. The Actinobacteria represented on average 1% of the total reads. Several other groups such as the Acidobacteria, Firmicutes, Gracilibacteria, Parcubacteria, Planctomycetes, and the Verrucomicrobia were also detected but at low read abundances (< 1%) and were grouped as “Other Bacteria” for plotting purposes. Analysis of the variance resulted in no significant differences in the contribution of the major taxa retrieved by each sequencing methodology.

A similar pattern was observed for picoeukaryotes (0.2–3 μm size fraction). The 556,143 clean reads were clustered into an OTU table at 97% similarity that contained 768 OTUs; from those only 37.1% were shared between the two methodologies, but these represented the vast majority of reads (96.4%; Figure 2). OTUs recovered with only one of the sequencing methodologies represented very rare members. In fact, as for bacteria, we found very good correlations when comparing the relative abundance of the different taxa both at the OTU level or clustering them at the taxonomic group level (R = 0.84, p < 0.001, and R = 0.91, p < 0.001, respectively; see Figure 3). The picoeukaryotic OTUs were classified into 70 class-level groups. The taxonomic affiliation was dominated by four groups that accounted on average for >55% of the total number of reads within the picoeukaryotic dataset: Mamiellophyceae (19% of the reads, dominated by Micromonas OTUs [97% of Mamiellophyceae]), Dinophyceae (17%), MALV-II (10%), and Cryptophyceae (10%). Other less abundant groups included MALV-I, Chlorarachnida, Picozoa, Prasinophyceae, Dictyochophytes, Chlorodendrophyceae, MAST-3, and Pelagophytes. The remaining 58 taxonomic groups presented very low relative abundances (< 1.1%) and were grouped as “Other Eukaryotes.” No statistically significant differences in the relative abundance retrieved by 454 or Illumina for the difference groups were found.

In order to further explore whether the sequencing methodology had an influence on the bacterial and picoeukaryotic diversity, we calculated various widely used indices of alphadiversity: the Chao 1 index for richness, and the Shannon and Simpson indices for diversity estimation (Hill, 1973; Magurran, 1988; Chao and Lee, 1992; Figure 4). Analysis of variance showed no significant differences between sequencing platforms for any of the indices tested, neither for Bacteria nor for Picoeukaryotes (P > 0.05). Additionally, to infer the variation of the microbial assemblages along the gradients, that is, beta diversity, the Bray–Curtis dissimilarity index was used on community composition. Dissimilarity matrices were constructed based on the relative abundance of each OTU. The distance between samples and replicates was visualized using hierarchical clustering. The results show that, in general, replication was good within each sequencing platform, but replicates sequenced using the same methodology were more similar among each other, which indicates that the sequencing chemistry has a certain influence on the community composition observed (Figure 5). For Bacteria, the samples grouped according to station regardless of the sequencing platform, except for Stations 3 and 4, which grouped by method, indicating that the platform introduces errors and artifacts to a certain extent at the OTU level. A similar trend was observed for the picoeukaryotic dataset, in which samples grouped by station, and in general were more similar among replicates subjected to the same methodology. Yet, in one case, the replicate obtained by Illumina (Illu−4a) was fairly different to the rest of the replicates from the same station. The number of OTUs in Illu-4a sample was much lower than in the other three replicates of Station 4 (one from Illumina and two from 454), indicating some biases in amplification or sequencing of this specific sample.

In order to obtain direct descriptors of the bacterial and picoeukaryotic diversity of plankton assemblages, we compared the diversity retrieved along the inshore-to-offshore gradient (Figure 1, Table 1). We observed significant differences in alphadiversity between stations (Figure 6). In particular, significant differences for Chao1 and Simpson indices were found for the bacterial dataset. The Chao1 showed a clear inshore-to-offshore decrease, whereas the Simpson index showed higher values in the transition zone from the coastal to the offshore station (Figure 6). The Shannon index showed a similar trend to the Simpson index but the differences detected were not significant. Interestingly, eukaryotic picoplankton showed a different trend. Whereas the values of Chao1 were quite constant along the gradient, the Shannon and Simpson diversity indices increased from coast to offshore. In fact, statistical analyses (ANOVA) confirmed significant differences, particularly between stations 1 and 5.

In addition, clear inshore-to-offshore changes in community composition were found both for bacteria and eukaryotic picoplankton. The larger differences were detected between Station 1 (coastal) and Station 5 (offshore), whereas a transition in community composition was observed at intermediate stations (Figure 3). In the case of Bacteria, some phylogenetic groups (Phylum, Class, and Order levels) showed a clear increase in their abundance from coast to offshore. These include the phylum Actinobacteria and the orders Rhodospirillales, Rickettsiales, and SAR11 within the class Alphaproteobacteria. An opposite trend was observed for the phylum Bacteroidetes, order Rhodobacterales (Alphaproteobacteria), and orders Alteromonadales and Oceanospirillales of the Gammaproteobacteria. Phylum Cyanobacteria were small contributors to community composition in the coastal station and peaked at Station 3 coinciding with the highest value of Chl a. The greater differences were observed for the order Alteromonadales which represented >25% of the reads in the coastal station and decreased to almost nil in the offshore station. Conversely, the SAR11 clade increased from 1 to >20% of the reads along the transect. Analyses of variance confirmed significant differences between stations for all the above-mentioned groups (details now shown).

The picoeukaryotic community also changed along the gradient being likewise Station 1 the most different from Station 5. The lineage Mamiellophyceae showed similar high relative abundance (>20%) in all stations except in Station 5, where they were virtually absent. The relative abundance of Cryptophyceae increased from Station 1 to Station 3 and then decreased toward offshore stations. Dinophyceae, Dictyochophytes, marine alveolates (MALV-II and MALV-I), and Stramenopiles showed increasing contributions along the transect. Contrarily, Ciliophora were important contributors only in the coastal station. Other groups presented quite constant contributions in all stations (Picozoa, Prasinophyceae, Dictyochophytes; Figure 3).

Potential “indicator species” were explored by calculating the indicator value (IndVal; Dufrêne and Legendre, 1997; Podani and Csányi, 2010) which identifies indicator species based on species (or OTU) fidelity and relative abundance, both for bacterioplankton and eukaryotic picoplankton. The IndVal of a species is a popular measure to express species importance in community ecology. Its potential to measure species explanatory power and to reflect environmental quality has been explored in biodiversity surveys (Gevrey et al., 2010; Keith et al., 2012; Lumbreras et al., 2016). We classified the stations into three categories, i.e., coastal (Station 1), transition (Stations 2, 3, and 4), and offshore (Station 5) and searched for indicator OTUs. We found 114 bacterial OTUs with significant IndVal values, potentially useful as indicator species. However, we considered only those OTUs showing (i) IndVal values >0.3, as this is the value that has been proposed to be a good threshold for habitat specialization (Dufrêne and Legendre, 1997), and (ii) overall relative abundance >1% since the potential as indicator species of rare OTUs is questionable taking into account the differences found between sequencing methods for the rare OTUs and the known biases of the PCR-based methodologies (Polz and Cavanaugh, 1998; Acinas et al., 2005). After selection, the list was reduced to 23 bacterial OTUs. We found OTUs with explanatory power for all three categories. The OTU with higher IndVal value was affiliated to a species of Gammaproteobacteria (Marinobacterium) and was indicative of coastal waters. On the contrary, alphaproteobacterial members of the SAR11 clade were explanatory for offshore waters and mainly Bacteroidetes for the transition zone. Within picoeukaryotes, a total of 164 OTUs with significant values were found but after filtering the table using the same criteria only 13 OTUs were retained. Most of them were explanatory for Station 5 in offshore waters. However, the indicator OTU presenting a higher contribution, OTU1, was classified as Micromonas pusilla, and was indicator for coastal waters in agreement with previous reports that have shown the preference of Micromonas species for coastal waters (Not et al., 2005, 2008). Overall, IndVal values for picoeukaryotic OTUs were lower than for bacterial OTUs. In both cases, the highest values were associated to rare species (details not shown) that were discarded based on abundance data. The selected IndVal scores and associated OTUs are listed in Table 2.

In addition to exploring potential “indicator species,” we explored the microbial profiles as possible descriptors of environmental status. That is, analyzing the relative abundance of the most abundant phylogenetic groups in each sample in relation to the degree of impact. The transect analyzed off the coast of Barcelona reflects a decreasing gradient of human impact from inshore (Station 1) to offshore (Station 5) which is somewhat reflected in the concentration of inorganic nutrients (see Table 1). The analysis of changes in community composition along the gradient together with the OTUs showing highest IndVal scores suggested the exploration of the ratios between taxa as potential indices of ecosystem health status. Interestingly, we found strong positive and negative correlations between the relative abundance of different bacterial groups as well as the ratio between taxa and the concentration of nutrients. The strongest correlation detected was a positive correlation between the relative abundance of Alteromonadales and all nutrients measured (phosphate, nitrite, nitrate, ammonium, silicate, R > 0.96, p < 0.0001). Likewise, the ratios Alphaproteobacteria/Gammaproteobacteria, Alteromonas/SAR11, and Alteromonas + Oceanospirillales/SAR11 were strongly correlated to nutrient concentration (R > 0.90, p < 0.0001). For picoeukaryotes, we found significant correlations between the relative abundance of certain taxa and the nutrient load, yet these correlations were in general weaker than for bacteria. The strongest positive correlations were found for Ciliophora and all nutrients (R = 0.85–0.88, p < 0.0001). Significant negative correlation between Chlorarachnida and nitrite (R = 0.79) and nitrate (R = 0.75) as well as between Dinophyceae and phosphate (R = 0.71) were also observed.

Up to date, several studies have investigated the potential biases on the estimations of richness and evenness in microbial communities associated with the primer selection and the PCR step in amplicon-based studies (Acinas et al., 2005; Hong et al., 2009; Engelbrektson et al., 2010; Parada et al., 2015). However, since cross-platforms studies are rare, currently it is unclear whether the inherent differences in chemistry and sequencing protocols will affect the quality of the sequences and the estimates of genetic diversity and community structure. Furthermore, despite variability is known to be introduced during sample manipulation, PCR amplification and sequencing, the numerous studies on microbial diversity using HTS lack analysis of replicates (Prosser, 2010). For these reasons, we compared the two most frequently used HTS platforms, the Roche 454 FLX Titanium, and the Illumina MiSeq, on a set of DNA samples obtained from an inshore-to-offshore transect in the coast of Barcelona. Additionally, we explored the reproducibility of the results by sequencing replicates. Overall, the platforms provided a comparable view of the marine picoplankton communities but some differences were found when comparing the datasets at the OTU level.

Different HTS platforms produce millions of short sequence reads, which vary in length. It is known that sequence length can impact diversity estimates (Claesson et al., 2010). Nowadays, pair-end Illumina can produce up to 300 bp nucleotide reads, and thus is feasible to do a careful compassion with 454 using the same primer set, providing the same amplicon length, and thus distinguish the performance of both methodologies based only in potential differences in the chemistry of the sequencing. Here, we found that the sequencing methodology does not significantly influence estimates of alphadiversity. No significant differences in Chao1, Shannon, and Simpson indices were found between platforms. A recent study comparing the Illumina and 454 platforms to study bacterial diversity via 16 S rRNA gene amplicons in sediments and soda lakes also found that both methodologies performed in a similar manner and that the general trends in alphadiversity were conserved with the exception of evenness estimates where correspondence between methods was low (Sinclair et al., 2015). It is known that the OTU clustering method can influence the estimates of diversity (Edgar, 2013; Flynn et al., 2015; Sinclair et al., 2015). We used the UPARSE algorithm, which offers an improved accuracy compared to other methods, resulting in fewer OTUs likely closer to the expected number of species in a community (Edgar et al., 2011). Using this methodology may have reduced the influence of sequencing and amplification artifacts and resulted in comparable estimates of diversity by the two sequencing methodologies.

We did observe some differences for betadiversity, that is the variation of the microbial assemblages along the transect, despite the trends identified were in general similar for both methodologies. Replication was good within each sequencing platform but in general replicates were more similar among each other depending on methodology, revealing thus a certain influence of its chemistry. We found that the bacterial communities in Stations 3 and 4 were more similar depending on the method indicating that the platform introduces biases. Oceanographic conditions were quite similar between these two stations (Table 1), and therefore microbial communities could be expected to be fairly similar. Sampling artifacts associated with random sampling (Zhou et al., 2008), PCR biases (Polz and Cavanaugh, 1998; Acinas et al., 2005) or errors directly related to the performance of the technology per se (Berry et al., 2011; Schirmer et al., 2015) can occur at any time, but when comparing samples, the impact of these artifacts will depend on the similarity among those samples. In this case, it is feasible to assume that the potential artifacts associated to the methodology overwhelmed the natural differences between the communities in these closer stations. For picoeukaryotes, in general samples grouped by station as expected, indicating that the sequencing biases, if any, were minor. However, there is one replicate from Station 4 that differs substantially from the other replicates. Problems during PCR amplification or degradation of the DNA could explain this difference.

Venn diagrams showed that less than half of the total OTUs were equally retrieved by both methodologies. However, the non-shared OTUs correspond to very rare contributors of these microbial communities. The concept of the rare biosphere has attracted a lot of attention in the last years (Pedrós-Alió, 2012; see reviews by Lynch and Neufeld, 2015). Microbial communities are dominated by a small number of species that account for most of the biomass and a large number of species that are represented by only a few individuals (rare members). The development of HTS has allowed accessing at least some of these rare microbial species. However, it is known that some of the rare OTUs retrieved in microbial diversity surveys correspond to sequencing errors (Kunin et al., 2010). We discarded the singletons (OTUs represented by a single sequence in the whole dataset) to avoid potential artifacts in diversity estimates. Nevertheless, still over half of the OTUs were only retrieved by one methodology. Part of it can be explained because rare OTUs may or may not appear in a dataset only by random chance but we cannot discard that part of this diversity is due to sequencing errors. For that reason, for the purpose of finding indicator species, we decided to focus only on those OTUs that represented >1% of the total relative abundance. Regardless of the differences in the rare OTUs, the two sequencing technologies revealed very similar profiles when grouping OTUs at the class and family levels (Figure 3). Relative taxa abundances were consistent across technologies and thus, the view of the community composition was fairly comparable. The results show that, due to the improvement in the length of Illumina sequence reads, Illumina tags offer similar classification efficiencies than 454 tags at a much lower cost (Glenn, 2011), being therefore a cost efficient approach for biodiversity monitoring.

Diversity and trophic state are two quality descriptors for evaluating ecosystem function in the MSFD. Despite the main goal of this work was to compare HTS methodologies for biodiversity monitoring, we further explored whether picoplankton biodiversity can be used as an alternative indicator of environmental status. The Mediterranean Sea is a valuable paradigm to assess anthropic pressure, because of the contrasting nature of its offshore and coastal areas. The offshore waters of the Mediterranean Sea are among the most oligotrophic areas of the world. In these waters, nutrient availability is low and inorganic phosphorus concentrations limit primary production. On the contrary, coastal areas are nutrient rich, as they receive river discharges, runoff from populated areas, and submarine groundwater, but they are also influenced by offshore oceanographic conditions. The coastal marine zone is therefore a transitional area characterized by strong physical, chemical, and biological gradients that extend from land to sea. Here, biological production is closely coupled to processes that deliver nutrients to surface waters. Anthropogenic forcing clearly influences the absolute availability of these nutrients and their stoichiometry, both of which impact phytoplankton productivity and species composition (Camp et al., 2015). The studied transect is expected to have a decreasing degree of anthropogenic pressure as the distance from the coast increases (from Station 1 to 5). Concentration of inorganic nutrients, as indication of eutrophication, showed indeed a decreasing concentration. We determined common alphadiversity indices as possible descriptors of the environmental status since pressures can lead to changes in microbial composition (Torsvik et al., 2002; Smith and Schindler, 2009) and those could reflect variations in biodiversity. For bacteria, Shannon and Simpson indices showed a similar trend with higher values at intermediate stations of the transect. The observed trend could be explained by the “intermediate disturbance hypothesis” (Connell, 1978), which suggests that intermediate intensity of disturbance maximizes diversity, and therefore systems with low and high disturbance, such Stations 5 and 1 in terms of nutrient load, can harbor similar levels of diversity. In any case, as previously observed in other systems (Garrido et al., 2014) these indices do not seem promising as indicators to asses environmental status. Contrarily, a clear decrease in richness was observed from coast to offshore. A sharp decrease of richness from coastal to offshore locations in the NW Mediterranean has been previously documented (Pommier et al., 2010). On the other hand, an increase in Shannon and Simpson indices was observed along the transect for picoeukaryotes, indicating a higher diversity in more oligotrophic stations (Cheung et al., 2010). Furthermore, the most abundant OTU in all stations but Station 5, Micromonas, is known to be more common in coastal areas than open ocean (Not et al., 2005), possibly related to higher nutrient load in coastal waters. The results found here suggest that it may be worth exploring the links between bacterial and picoeukaryotic diversity and environmental status on coastal waters over time and space covering a wide range of impacts.

Traditionally, several species of plants and animals have been and still are being used as indicator species for different types of pollution in monitoring programs (Borja et al., 2000, 2008; Ferrat et al., 2003; Montefalcone, 2009; Marbà et al., 2013). Likewise, plankton indicators have been proposed for diagnoses of ecosystem state (Beaugrand, 2005). Most studies have focused on species of zooplankton (i.e., Calanus finmarchicus) or some phytoplankton bloom-forming species. For example, Phaeocystis sp. produces spring blooms in the North Sea which magnitude might indicate an excess of available N or P in relation to dissolved silica and thus, is considered and indicator for eutrophication (Tett et al., 2007). However, several flaws in the usefulness of using large phytoplankton to reflect significant pressure-impact relationships have been identified (Cloern and Jassby, 2008, 2010; Camp et al., 2015). Bacterial and eukaryotic picoplankton constitute the smallest but most abundant organisms of plankton and are key players in ecosystem functioning. Since disturbances can affect community structure and ecosystem functioning, the smallest members of marine plankton may be crucial in understanding the magnitude of these disturbances particularly because of their fast response to environmental change. In fact, microorganisms have been already proposed as indicators of marine environmental quality, and not only the presence of pathogens such as E. coli, commonly used as indicator of fecal contamination, but in relation to biodiversity and ecosystem functioning (Caruso et al., 2015). Here we tested the Indicator Species Value from Dufrêne and Legendre (1997) in the different sampled stations. The IndVal identifies indicator species based on OTU fidelity and relative abundance. Different bacterial and picoeukaryotic OTUs showed high scores for Stations 1 and 5, as well as for intermediate stations and could represent potential “indicator species.” Alternatively to “indicator species,” we explored the potential of using the abundance of certain taxa and the ratio between different groups of microorganisms as an alternative indicator of environmental status. These indices may also offer ecological information (i.e., species relative composition). In fact, this approach has been explored in other ecosystems; for example in reclaimed waters, the ratio between the Bacteroidetes, Gammaproteobacteria, and Nitrospira/Betaproteobacteria (BGN:β) seems a possible alternative indicator of water quality (Garrido et al., 2014). We tested the correlation of different taxa and the degree of eutrophication (i.e., nutrient concentration) and found significant correlations between certain picoeukaryotic taxa e.g., Ciliophora, and nutrient load; this taxa has been found previously in high abundances in eutrophic waters (Romari and Vaulot, 2004). Yet, the strongest correlations were with the ratio of Alphaproteobacteria/Gammaproteobacteria, Alteromonas/SAR11, and Alteromonas + Oceanospirillales/SAR11. Whether these “indicator species” and indices can be used as robust alternative indicators of environmental status remains to be explored in different locations subjected to contrasting pressures and over time. The challenge is to discriminate between antropogenic-induced changes and the confounding effects of the natural variability of the marine environment.

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. The handling Editor declared a collaboration with the authors and states that the process nevertheless met the standards of a fair and objective review.