Discussion
In this study we investigated if current sequencing technology and eRNA
shotgun sequencing has the power to differentiate changes in
biodiversity of multiple microeukaryotes in bulk sediment samples. We
focused on nematodes, forams, and ciliates which are useful
bioindicators and known to change in diversity and community composition
in relation to environmental change (Gong et al., 2015; Ingels et al.,
2009; Martins et al., 2015; Jan Pawlowski et al., 2014; Ristau et al.,
2015). The results showed a difference in community structure for each
of the communities along the OC gradient in the study area. For example,
the non-selective deposit feeding nematode genera Sabatieria andAxonolaimus (Schratzberger, Warr, & Rogers, 2007) had a higher
relative abundance at the High OC stations. Potentially this could be a
beneficial feeding strategy at the deeper stations where the sediment
consists mainly of decayed organic particles and bacteria as food (and
is reflected in the nematode feeding type analysis; Figure. 6).Sabatieria are typical nematodes found in organic rich sediments,
and have been identified with other non-selective deposit feeders such
as Daptonema (Armenteros et al., 2009; Montagna & Harper, 1996;
Schratzberger, Warr, & Rogers, 2006). Interestingly, the generaDaptonema and Desmolaimus (also a non-selective deposit
feeder (Schratzberger et al., 2007)) had a higher relative abundance at
the Low OC stations, and they could potentially be directly competing
with other nematodes with a similar feeding behavior (such asSabatieria ). The Low OC stations had more nematodes classified as
selective deposit feeders and predator/omnivores, suggesting different
kinds of food (e.g. algae) and more competition in these environments.
These nematode genera have previously been detected in other basins the
Baltic Sea using 18S rRNA gene metabarcoding (Elias Broman et al.,
2019), and here their presence was confirmed by shotgun sequencing.
The foram genera Elphidium and Rhizammina showed
contrasting patterns in the dataset, with Elphidium having higher
relative abundance at High OC stations, and Rhizammina at the Low
OC stations. The species Elphidium excavatum has
previously been found to be dominant in a lagoon located on the Atlantic
coast of Portugal with high content of OC (2–4 % OC) compared to
adjacent bays with lower OC (Martins et al., 2015). BothElphidium and Rhizammina are known to exist in the
south-western Baltic Sea (Frenzel, Tech, & Bartholdy, 2005; Schweizer,
Polovodova, Nikulina, & Schönfeld, 2011), and to our knowledge, this is
the first study using molecular data to investigate diversity of forams
in the Gulf of Finland. Even though SILVA is one of the recommended
options to classify 18S rRNA sequences of protists (Creer et al., 2016)
we were still surprised to see many differences in classifiedForaminifera genera between the SILVA and NCBI NT databases. For
example, SILVA reported a high relative abundance of genera (e.g.Calcarina , up to 35% in the offshore stations) never detected in
the Baltic Sea (Supplementary Data 5). Almost 100 foram species have
been reported from the south-western Baltic Sea, but there is still a
large scarcity of studies investigating forams in the central and
northern Baltic Sea (Frenzel et al., 2005). It is therefore possible
that the differences we observed between databases are due to a limited
number of 18S rRNA Baltic Sea-like foram sequences in SILVA. There are
specific foram databases such as the PFR2 (Morard et
al., 2015), but this database is focusing on oceanic forams featuring no
data from the Baltic Sea (e.g. no data on Elphidium andRhizammina ). Nevertheless, we report good results derived from
the NCBI NT database.
Finally, we did not detect a difference in alpha diversity for ciliates
between High OC and Low OC stations. Neither was there a difference
between genera with a high relative abundance (exceptCryptocaryon ), and so the difference observed in beta diversity
would therefore have been due to differences in low abundant genera. The
ciliate genus Cryptocaryon was more prominent at the Low OC
stations;a marine genus known to include parasitic species targeting
fish (Wright & Colorni, 2002). However, low-saline (5–7 ppt) variants
of Cryptocaryon have previously been described (Yambot, Song, &
Sung, 2003). Cryptocaryon -like ciliates have previously been
detected phylogenetically in the more saline (~14 ppt)
deeper waters of the Baltic Sea (Stock, Jürgens, Bunge, & Stoeck,
2009), and both SILVA and NCBI NT (and manually checking classified
sequences via BLAST) confirmed Cryptocaryon -like ciliates in our
samples. Potentially because the Low OC stations were located in shallow
areas closer to the shore, the Cryptocaryon -like ciliates
detected in this study might be adapted to lower salinities
(~7 ppt) and related to host organisms residing in these
more diverse and euphotic habitats. In addition, considering that the
summer heatwave of 2018 was one of the most intense ever recorded in the
study area (Humborg et al., 2019), the warmer waters might attractCryptocaryon -like ciliates which are typically more common at
temperate and tropical temperatures (Colorni & Burgess, 1997). As far
as we are aware, this is the first time such species have been reported
from the Gulf of Finland based on molecular data.
Limitations of this study include the relatively small sample size for
RNA extraction (2 g, i.e. ~2 ml for high porosity
sediment). A previous study investigating the effect of sample size on
diversity estimations of microeukaryotes and metazoans using
metabarcoding, found that larger volumes of sediments are necessary to a
accurately estimate small-scale spatial heterogeneity in biodiversity
(Nascimento, Lallias, Bik, & Creer, 2018). Here we used the available
commercial extraction kit that could process the largest input of
sediment volume. It will be useful for future studies to develop larger
kits for eRNA extraction and investigate similar effects of sample size
on biodiversity assessment based on eRNA. Nevertheless, even with this
sample size we report clear differences in biodiversity of multiplied
communities including metazoans along the environmental gradients.
Organisms contain multiple ribosomes and while this is not an issue when
analyzing species richness or beta diversity with a presence/absence
index, for community composition organisms with a large number of
ribosomes could skew the proportions. This issue also exists in DNA
metabarcoding studies, with many eukaryotic organisms carrying multiple
genome copies per cell (i.e. polyploidy) (Edgar, Zielke, & Gutierrez,
2014) and for prokaryotes that can also be polyploid (Soppa, 2017). The
cost of shotgun sequencing on the latest platform is still quite high
($8815 USD per Illumina NovaSeq S4 lane) compared to metabarcoding of
marker genes ($1373 USD per Illumina MiSeq flowcell, prices based on
the cost at SciLifeLab, Stockholm, Sweden in 2019). However, there has
been a large decrease in sequencing cost over the past 20 years
(Wetterstrand, 2020) and if this trend continues, alongside streamlined
bioinformatic protocols, large scale eRNA shotgun studies could be a
future possibility in biomonitoring programmes.
Shotgun sequencing “catches” all organisms in the sample including
both prokaryotes and eukaryotes (Zepeda Mendoza, Sicheritz-Pontén, &
Gilbert, 2015). Considering that rRNA has a short lifespan in the
environment (Blazewicz, Barnard, Daly, & Firestone, 2013) the approach
used in this study very likely mostly targeted only living (or very
recently dead) organisms, as well as eRNA derived from organisms
recently passing through the sampling area but not located in the
collected sediment (up to ~13h or longer in biofilm, see
S. A. Wood et al., 2020). This makes eRNA sequencing a useful method to
study benthic communities, considering that a substantial portion of
sediment consists of long-lasting dead organic matter (Burdige, 2006).
For example, forams are a known microfossil group and with the use of
DNA extraction it has been possible to study the ancient DNA of these
organisms (Lejzerowicz et al., 2013). Shotgun sequencing of eRNA has
been used in a wide variety of marine studies, including investigations
of prokaryotic communities (Elias Broman, Sachpazidou, Pinhassi, &
Dopson, 2017; Elias Broman, Sjöstedt, Pinhassi, & Dopson, 2017;
Klindworth et al., 2014; Urich et al., 2014), marine viruses, (Culley,
Lang, & Suttle, 2006), sediment eukaryotic metatranscriptomes (Elias
Broman, Varvara, Dopson, & Hylander, 2017), and old marine groundwaters
in the deep terrestrial biosphere (Lopez-Fernandez et al., 2018).
However, there is paucity of studies using eRNA to assess the
biodiversity of microeukaryotes in sediment. Although 18S rRNA
metabarcoding has gained popularity to investigate such communities (see
e.g. Birrer et al., 2018; Comeau, Lagunas, Scarcella, Varela, &
Lovejoy, 2019; Rodríguez-Martínez et al., 2020), we have here shown that
eRNA shotgun sequencing is also a viable approach to detect differences
in diversity and community compositions for multiple communities as a
response to different environmental conditions. Even though not directly
compared in this study, shotgun sequencing avoids the main limitations
of metabarcoding such as 1) PCR primers amplifying species with
different rates and 2) the amount of cycles and type of polymerase used
influencing diversity and community composition (Kelly et al., 2019;
Nichols et al., 2018). In addition, eRNA shotgun studies also provides
information on all organisms over a large range of trophic stages in the
sediment, and it is also possible to study the RNA transcripts of
expressed genes to estimate oxidation and reduction processes of
nutrients (mainly related to the metabolism of prokaryotes) (see e.g.
Elias Broman et al., 2017). Finally, common shotgun sequencing
bioinformatic pipelines are intricate and include many different
software which increases the complexity of the data analysis. Here,
conversely, we followed a protocol with few and straightforward steps,
including: 1) Quality trimming (removal of Illumina adapters and phiX
control sequences, quality trimming, and verifying final quality); 2)
extraction of SSU rRNA sequences from the dataset; 3) taxonomic
classification of the SSU rRNA sequences using Kraken2 against the NCBI
NT and SILVA databases; and 4) estimation of relative abundance at genus
level using Bracken2 (for more details on the kraken2+bracken2
combination see D. E. Wood et al., 2019). These new bioinformatic tools
makes it less daunting and possible to classify large datasets
containing hundreds of millions of sequences within minutes to hours
which would previously have taken several weeks using traditional
aligners. As such this approach is closer to current metabarcoding
bioinformatic pipelines with relatively straightforward steps (see e.g.
the DADA2 pipeline Callahan et al., 2016).