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).