Corresponding author: Sirje Sildever ( sirje.sildever@taltech.ee ) Corresponding author: Satoshi Nagai ( snagai@affrc.go.jp ) Academic editor: Jianghua Yang
© 2021 Sirje Sildever, Peeter Laas, Natalja Kolesova, Inga Lips, Urmas Lips, Satoshi Nagai.
This is an open access article distributed under the terms of the Creative Commons Attribution License (CC BY 4.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Citation:
Sildever S, Laas P, Kolesova N, Lips I, Lips U, Nagai S (2021) Plankton biodiversity and species co-occurrence based on environmental DNA – a multiple marker study. Metabarcoding and Metagenomics 5: e72371. https://doi.org/10.3897/mbmg.5.72371
|
Metabarcoding in combination with high-throughput sequencing (HTS) allows simultaneous detection of multiple taxa by targeting single or several taxonomically informative gene regions from environmental DNA samples. In this study, a multiple-marker HTS approach was applied to investigate the plankton diversity and seasonal succession in the Baltic Sea from winter to autumn. Four different markers targeting the 16S, 18S, and 28S ribosomal RNA genes were employed, including a marker for more efficient dinoflagellate detection. Typical seasonal changes were observed in phyto- and bacterioplankton communities. In phytoplankton, the appearance patterns of selected common, dominant, or harmful species followed the patterns also confirmed based on 20 years of phytoplankton monitoring data. In the case of zooplankton, both macro- and microzooplankton species were detected. However, no seasonal patterns were detected in their appearance. In total, 15 and 2 new zoo- and phytoplankton species were detected from the Baltic Sea. HTS approach was especially useful for detecting microzooplankton species as well as for investigating the co-occurrence and potential interactions of different taxa. The results of this study further exemplify the efficiency of metabarcoding for biodiversity monitoring and the advantage of employing multiple markers through the detection of species not identifiable based on a single marker survey and/or by traditional morphology-based methods.
bacteria, Baltic Sea, metabarcoding, phytoplankton, rRNA genes, seasons, zooplankton
Zoo-, phyto- and bacterioplankton are important components of marine food webs as grazers, primary producers, and decomposers (
Long-term monitoring of plankton communities is needed to identify and model the influence of interannual variability and global change, as suggested for the phytoplankton (
To identify multiple species simultaneously and analyze several samples at once, metabarcoding and high-throughput sequencing (HTS) has been applied to investigate plankton diversity (
In the Baltic Sea, the HTS-approach has been previously applied for investigating plankton communities along the salinity gradient (
Water samples from 5 m depth were collected from 9 stations in the Gulf of Finland (GoF), the Baltic Sea, from winter to late spring (2013 Dec – 2014 May; 4 stations), and from summer to autumn (2013 May-Oct; 5 stations; Fig.
DNA extracted from 0.2 and 5 µm filters was mixed in equal volume and used as a template for paired-end library preparation. Eukaryotic species were targeted by universal primers for the 18S (V7-V9 region; SSR-F1289-sn, F: TGGAGYGATHTGTCTGGTTDATTCCG; SSR-R1772-sn, R: TCACCTACGGAWACCTTGTTACG,
The purified PCR products were used as a template in the second-round PCR. The reaction mixture composition was the same as in the first PCR, however, the total volume of the mixture was 50 µl, 2 µl of DNA template, and the following primers were used: 5’AATGATACGGCGACCACCGAGATCTACAC-8 bp index-ACACTCTTTCCCTACACGACGC (forward) and 5’CAAGCA GAAGACGGCATACGAGAT- 8 bp index GTGACTGGAGTTCAGACGTGTG (reverse). The PCR conditions were as follows: initial denaturation at 94˚C for 3 min, 8–12 cycles at 94˚C for 15 s, 56˚C for 30 s, and 68˚C for 40 s. PCR amplification was verified by agarose gel electrophoresis, and the PCR products were purified using an Agencourt AMPure XP. The amplified PCR products were quantified, pooled in equal concentrations, and stored at -30 ˚C until sequencing. Sequencing on Illumina Miseq 300 PE platform was ordered from the FASMAC, Co. Ltd., Japan.
Treatment of obtained sequences, selection of operational taxonomic units (OTUs), and taxonomic characterization of OTUs was done according to the workflow described by
The selected OTUs were taxonomically identified as follows: for the BLAST search, a subset of NCBI non-redundant nucleotides (NT) consisting of sequences that satisfied the below-mentioned conditions was prepared. The sequences used to identify OTUs were downloaded from the GenBank based on the following search criteria: “ribosomal”, “rrna” or “rdna”, and the search excluded keywords associated with “metagenome”, “uncultured” and “environmental”. The sequences of retrieved GenBank IDs were downloaded from the nucleotide database at the NCBI FTP server on March 23, 2019, and were used to construct a template sequence database. Subsequently, the taxonomic characterization for each OTU was performed using a BLAST search (Cheung et al., 2010) in NCBI BLAST+2.2.26+ (
For eukaryotes, only OTUs with blast top hit similarity of ≥ 0.99% and a minimum query cover of 70% were included in the further analysis. OTUs associated with multiple records from the same genus were merged if the multiple records consisted of sp. of the same genus or when a single species was clustered together with other records identified as sp. from the same genus. For prokaryotes, only OTUs with > 50 total sequence reads were included in the further analysis. This rule was not applied for cyanobacteria that were identified to species level. The phylum and class/order level taxonomic groupings were also confirmed based on the World Register of Marine Species (www.marinespecies.org) or the UniProt taxonomy database (www.uniprot.org). For OTUs associated with Dinophyceae and Bacillariophyceae, the class level grouping is shown instead of phylum level as both represent important and dominant phytoplankton groups. The autotrophic ciliate, Mesodinium rubrum, is included both in the phyto- and zooplankton results to show its dynamics concerning both groups as it is usually analyzed together with phytoplankton.
The presence of toxin-producing HAB species among the detected OTUs was investigated by comparing the list of OTUs associated with phytoplankton species against the IOC-UNESCO Taxonomic Reference List of Harmful Micro Algae (
The numbers of shared OTUs among different sampling years were counted and depicted as a Venn diagram by using the “vegan” package (
During winter-spring, salinity ranged from 5.11 g kg-1 to 6.56 g kg-1 with higher salinities at stations A1 and A4 located in the western GoF (Fig.
The number of eukaryotic OTUs detected was 1755, 551, and 229 based on 18S, 28S_E (universal eukaryote), and 28S_D (universal dinoflagellates) markers, respectively. At ≥ 0.99% BLAST top hit similarity level and a minimum of 70% query cover length, the highest number of OTUs (n=314) was identified by 18S followed by 28S_E (95) and 28S_D (76) primers (Table
Overview of eukaryotic OTUs associated with different phyla (> 10 OTUs).
Phyla | 18S | 28S | 28S_D |
---|---|---|---|
Arthropoda | 16 | 1 | |
Ascomycota | 1 | 26 | |
Bacillariophyta | 36 | 10 | 18 |
Basidiomycota | 6 | 12 | 3 |
Cercozoa | 10 | 4 | 1 |
Chlorophyta | 24 | 7 | 2 |
Ciliophora | 31 | 4 | 2 |
Cryptophyta | 11 | 2 | 1 |
Dinophyceae | 86 | 11 | 37 |
Haptophyta | 15 | 6 | 5 |
Ochrophyta | 19 | 1 | |
eukaryote SCGC | 11 | ||
Others* | 48 | 13 | 5 |
Total nr. of OTUs | 314 | 95 | 76 |
In total, 24 and 29 unique OTUs associated with micro-and mesozooplankton were detected. The majority of those were identified only based on the 18S rRNA gene (18 and 25 OTUs, respectively), whereas the 28S_D marker facilitated the detection of 2 unique mesozooplankton OTUs (Suppl. material
From phytoplankton species, 15 OTUs associated with common or dominant species (
Based on 16S universal primers, the presence of 893 prokaryotic OTUs (> 50 sequences) from 22 phyla and 63 order/class was detected (Table
Class/Order | Nr. of OTUs |
Alphaproteobacteria | 202 |
Bacteria candidate phyla | 11 |
Bacteroidetes/Chlorobi group | 173 |
Betaproteobacteria | 57 |
delta/epsilon subdivisions | 32 |
Gammaproteobacteria | 124 |
Micrococcales | 10 |
Oligoflexia | 12 |
Synechococcales | 35 |
unclassified Actinobacteria | 38 |
unclassified Bacteria | 73 |
unclassified Verrucomicrobia | 15 |
Others* | 111 |
Total nr. of OTUs | 893 |
3.3.1 Eukaryotes
From all the class/phyla detected, Bacillariophyta, Ciliophora, Chlorophyta, Dinophyceae, and Picozoa were present during all the sampling months based on different eukaryotic markers. Based on 18S and 28S_E, the relative sequence abundances of Dinophyceae generally increased from February with the highest values from March to May (Figs
Relative sequence abundances of different phyla detected by 18S in different sampling months. TOP 5 most abundant phyla/groups are shown for each station and sampling occasion (X-axis). Y-axis displays the relative sequence abundances for each phyla/group. Sequences belonging to the phyla/group that were not among the TOP 5 phyla/group for the particular sampling occasion are displayed as “Others”.
From 53 unique micro-and mesozooplankton OTUs detected by different markers, five were present in all sampling occasions and stations. Those OTUs were associated with ciliates (Askenasia sp., M. major, M. rubrum, S. lemnae) and a Picozoa, P. judraskeda (Fig.
From the 31 common/dominant/toxin-producing phytoplankton species detected, about 6 displayed a distinct seasonal pattern by being present mainly in spring or summer, e.g. Chaetoceros holsaticus, Heterocapsa triquetra, G. corollarium, N. spumigena, Pauliella taeniata, and Protoceratium reticulatum (Fig.
Relative sequence abundances of common/dominant/toxin-producing species detected based on different markers. The X-axis represents stations and sampling months, Y-axis displays species detected. Species marked with bold are known to produce toxins (Moestrup, et al. 2020). Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months. *no toxicity confirmed from the Baltic Sea (
Relative sequence abundances of zooplankton species detected based on different markers. The X-axis represents stations and sampling months, Y-axis displays species detected. Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months.
3.3.2 Prokaryotes
OTUs associated with the FCB group and Proteobacteria had high relative sequence abundances in the majority of the stations and sampling occasions (Fig.
Relative sequence abundances of different phyla detected by 16S in different sampling months. TOP 5 most abundant phyla/groups are shown for each station and sampling occasion. (X-axis). Y-axis displays the relative sequence abundances for each phyla/group. Sequences belonging to the phyla/group that did not belong to the TOP 5 phyla/group for the particular sampling occasion are displayed as “Others”.
Relative abundances of bacteria-associated OTUs (TOP 39). The X-axis represents stations and sampling months, Y-axis displays selected taxa. Order of taxa is derived from pair-wise similarity matrix based on r-values (Pearson). Station and month labels marked with “_p” indicate months with pooled samples. Different colors indicate associations with different phyla. The size of the bubbles indicates relative sequence abundances. Vertical lines differentiate between the sampling months.
Based on 18S, 28S_D and 16S markers, samples from winter-spring and autumn (December-April and October) clustered together on the NMDS plots (Fig.
Common/dominant or toxin-producing species displayed statistically significant (P < 0.05) correlations with 24 zooplankton species (Suppl. material
The aim of this study was to detect the diversity present and reveal changes in community composition during seasonal succession from winter to spring and summer to autumn. Based on four different markers, community composition, seasonal dynamics, and co-occurrence patterns were identified for bacterio-, phyto- and zooplankton. Two phytoplankton and 15 zooplankton species novel for the Baltic Sea were also detected.
4.1.1 Eukaryotes
In this study, the highest number of eukaryotic OTUs detected belonged to Dinophyceae (phytoplankton). An exception to this was the 28S_E marker, which detected the highest number of OTUs associated with Ascomycota (fungi). A high number of OTUs associated with Dinophyceae has also been reported in previous studies from the Baltic Sea as well as from other localities around the world (
4.1.1.1 Zooplankton
Based on the HTS-approach, 8 species of micro- and 7 species of mesozooplankton previously not recorded from the Baltic Sea were detected. Based on the zooplankton monitoring data from 2013/2014 (Estonian Env. Agency., 2021a), only 7 mesozooplankton species were detected by the traditional morphology-based approach compared to the 29 species detected by the HTS-approach. Microzooplankton is not a part of routine zooplankton monitoring, except for M. rubrum that is included in the phytoplankton analysis. The importance of including microzooplankton in the monitoring program has been recently emphasized due to their notable contribution (> 50 %) to the nano- and microplankton biomass (
In the case of mesozooplankton, improved detection of zooplankton species based on HTS-approach compared to the morphology-based detection has been reported before (
The HTS-approach also detected three of the most abundant mesozooplankton species in 2013/2014 based on the monitoring reports (Eurytemora affinis, Keratella quadrata, Limnocalanus macrurus; Estonian Env. Agency 2021a). However, some dominant species (Eubosmina maritima, Synchaeta baltica, Pseudocalanus acuspes, and Pleopsis polyphenoides) were not detected by any of the three eukaryotic markers used. This might be due to the lack (S. baltica) or low number (max. 6) of 18S rRNA gene sequences available in the database (
4.1.1.2 Common/dominant/toxin-producing species
From the OTUs associated with toxin-producing species (
Other new records of toxin-producing phytoplankton species were Dolichospermum planctonicum (previously Anabaena planctonica) and Pseudo-nitzschia pungens that have been previously reported from the southern Baltic, but not from the GoF (
From the spring bloom dinoflagellate species, G. corollarium from the DinoComplex and P. catenata, could be detected by the HTS-approach. However, two other members of the DinoComplex, A. malmogiense, and B. baltica were not detected or present at a low similarity level (<0.99%). In the case of the 18S rRNA gene, both species were classified together with other species, e.g. A. malmogiense was identified together with A. aciculiferum and B. baltica together with B. cincta and B. brevisulcata at ≥ 0.99 % similarity level. In the case of A. malmogiense, this can be explained by the lack of differences in the 18S rRNA gene sequences between the two species (
4.1.2 Bacterioplankton
Overall, class-level taxa prominent within the bacterioplankton community composition (BCC) and their seasonal succession, although providing a very general overview, complemented near-perfectly the previous findings from the GoF that have used the same size fractioned filtration system (
The highest number of OTUs were associated with Alphaproteobacteria and FCB group (Bacterioidetes/Chlorobi group), followed by Actinobacteria. In general, Cyanobacteria became more prevalent during summer and Flavobacteria during the spring bloom, although some members of the FCB group occur exclusively during summer bloom. Interestingly, the unicellular cyanobacteria, Synechocystis and Synechococcus, increased in relative sequence abundance also in April. Alphaproteobacteria were dominant at the surface layer throughout different seasons, compared to previous studies, SAR11 bacteria (Pelagibacterales) were less prevalent in their relative abundances, but they were still present throughout the year (
4.2.1 Zooplankton
From the microzooplankton species registered in this and a previous study from the northern Baltic Proper (
From the Mesodinium species detected, M. rubrum and M. major displayed almost identical appearance patterns. These two species have been identified based on their morphology and division into distinct clades based on rRNA genes and internal transcribed spacer region (
4.2.2 Common/dominant/toxin-producing species appearance patterns
Investigation of common/dominant/toxin-producing species appearances patterns based on their relative sequence abundances displayed similar patterns with 20 years of phytoplankton monitoring data (for the species identifiable under light microscopy; Estonian Env. Agency 2021c; Suppl. material
In this study, the majority of the common/dominant/toxin-producing species detected appeared in various months throughout the year, and 2/3 of those species had the highest relative sequence abundances within the months with the highest biomass reported (Estonian Env. Agency, 2021c). For example, all detected cyanobacterial species had higher relative sequence abundances in summer, which is consistent with the previous knowledge from published studies (
For species that were detected in this study, but not present in the long-term phytoplankton dataset, the appearance patterns were investigated based on literature or based on the genus level information from the phytoplankton dataset. For example, the toxin-producing dinoflagellate A. ostenfeldii was present during almost all the sampling months with the highest relative sequence abundances in July. However, in the Åland archipelago, the cells are present from May to September with the highest abundances in July, August, and September (
Contrary to A. anophagefferens, A. dexeteroporum was detected throughout the year with the highest relative sequence abundances in January, February, and June. Based on the previous studies, the species can also grow in a wide temperature range (0.6 ˚C to 28˚C), depending on the geographic origin of the strain (
Prymnesium parvum f. patelliferum is a form of P. parvum (
Based on the correlation analysis between the common, dominant, or toxin-producing phytoplankton and zooplankton species, the statistically significant (P < 0.05) associations probably reflect co-appearance and not interactions as the majority of the species from both groups were detected throughout the sampling period. This is further illustrated by the significant negative correlations between dinoflagellate, A. anophagefferens, and a ctenophore, Mertensia ovum, as the dinoflagellate was only present in the samples from July and August, whereas the ctenophore species was present in the majority of the samples, except during those months. A similar, but opposite pattern was displayed between a chrysophyte Plagioselmis prolonga and several zooplankton species with the cryptophyte missing from the July-August samples. The significant positive correlations may also reflect the zooplankton presence following the presence of its prey, e.g. copepod Limnocalanus macrurus and diatoms Melosira arctica and Skeletonema marinoi (
In the case of correlations between common, dominant, or toxin-producing phytoplankton species and bacteria, the majority of the significant correlations detected were positive, except for A. anophagefferens. However, as the negative correlations were detected between A. anophagefferens and bacteria belonging to various phyla and classes, the significant negative correlations might reflect differences in appearance patterns. On average, the highest number of bacterial OTUs from all stations were present from December to April (average number of OTUs: 506–603), whereas during the detection of A. anophagefferens in July-August, the average number of bacterial OTUs ranged around 380. The interactions between phytoplankton and bacteria can range from mutualistic to competitive or antagonistic (Cole, 1982;
The results of this study provide the first information on plankton diversity and appearance patterns obtained by HTS-approach over various seasons from the Baltic Sea. The HTS-approach facilitated the detection of 15 zooplankton species previously not known from the Baltic Sea, with more than half of those being microzooplankton species. In the case of phytoplankton, the presence of two toxin-producing species was also recorded for the first time from the Baltic Sea. The presence of a diverse bacterial community was also registered. The phyto- and bacterioplankton dynamics followed the seasonal patterns known from long-term monitoring and previous studies displaying the reliable detection of plankton appearance patterns with the HTS-approach. In the case of zooplankton species, the majority of the species were present throughout the sampling season, especially ciliates associated with the microzooplankton. The correlation patterns between the phytoplankton and bacteria/zooplankton reflected patterns in species appearances. Overall, the results of this study support the usage of HTS-approach for plankton biodiversity monitoring as well as for following patterns in their seasonal dynamics.
The captain and the crew of RV “Salme” are thanked for their assistance during the sampling campaigns. Prof. N. Lundholm, Drs. H. Hällfors, S. Lehtinen, A. Jaanus, and A. Enke are thanked for the discussion on the potential presence of Pseudo-nitzschia pungens in the northern Baltic Sea. This study was supported by grants from Japan Society for the Promotion of Science Short-term Postdoctoral Fellowship (PE18028) [SN, SS], a Grant-in-Aid (Marine Metagenomics for Monitoring the Coastal Microbiota) by the Ministry of Agriculture, Forestry, and Fisheries of Japan [SN], a Grant-in-Aid for Scientific Research (Kiban-B) by the Japan Society for the Promotion of Science (Grant number: 18KK0182)[SN], JST/JICA, Science and Technology Research Partnership for Sustainable Development (JPMJSA1705) [SN], institutional research funding (IUT 19-6) of the Estonian Ministry of Education and Research [IL, UL], Estonian Research Council grant PRG602 [UL], European Regional Development Fund and the program Mobilitas Pluss (MOBTP160) [SS].