Introduction
Biodiversity is decreasing globally due to human alteration and pollution of terrestrial and aquatic environments (Brondizio, Settele, Díaz, & Ngo, 2019). Essential ecosystem services affiliated with human health, such as availability of food, clean water, and recreational areas are dependent on biodiversity (Cardinale et al., 2012; Pan, Marcoval, Bazzini, Vallina, & Marco, 2013). In addition to the provision of ecosystem services, biodiversity losses have also been linked to a decrease in ecosystem stability (McCann, 2000). Anthropogenic pressure on coastal aquatic ecosystems by e.g. climate change, eutrophication and contaminant pollution threatens the diversity of many organisms in these systems (Pan et al., 2013). Anthropogenic pressure in coastal ecosystems should be taken seriously because coastal zones are transitional areas directly adjacent to human settlements between land and sea, and impacted areas are predicted to increase in both number and area with a continued climate change scenario (Levin et al., 2001; Rabalais, Turner, Díaz, & Justić, 2009). It is therefore essential to understand how the diversity of organisms living in coastal zones respond to changes in environmental gradients and anthropogenic pressure (Snelgrove, Thrush, Wall, & Norkko, 2014).
Biodiversity assessments of benthic macrofauna are commonly used in national monitoring programs, including coastal zones, to determine various ecological indexes (Pinto et al., 2009). However, microeukaryotes present in sediment such as meiofaunal nematodes (< 1 mm body size) are also known to react to e.g. eutrophication status (Ristau, Spann, & Traunspurger, 2015), and the composition and quantity of organic carbon (OC) (Ingels, Kiriakoulakis, Wolff, & Vanreusel, 2009; Pusceddu, Gambi, Zeppilli, Bianchelli, & Danovaro, 2009). Furthermore, because nematodes are known to have different feeding behaviours such as bacterivory, detritivory or algal feeding (Moens, Traunspurger, & Bergtold, 2006; Wieser, 1953) changes in nematode assemblages are therefore likely to affect food web dynamics and ecosystem function (e.g. Nascimento et al., 2019; Nascimento, Karlson, & Elmgren, 2008; Nascimento, Näslund, & Elmgren, 2012). Other arguments for including meiofauna such as nematodes in national monitoring systems include their high diversity, short generation time, and ubiquitous distribution (Kennedy & Jacoby, 1999). However, these organisms are often neglected in monitoring studies (Kennedy & Jacoby, 1999), likely due to financial reasons derived from time consuming activities such as sieving, sorting, and microscopic morphological analyses.
In addition, the protist phyla Foraminifera (henceforth forams) and Ciliophora (i.e. ciliates) are well-studied as bioindicators of environmental state of aquatic ecosystems. The diversity and community composition of forams are known to change with anthropogenic pollution, fish farming, and decreasing water quality (Damak, Frontalini, Elleuch, & Kallel, 2019; Frontalini & Coccioni, 2011; Jan Pawlowski, Esling, Lejzerowicz, Cedhagen, & Wilding, 2014; Raposo et al., 2018; Uthicke & Nobes, 2008), and similar to nematodes, OC enrichment of the sediment also influences the diversity of forams (Martins et al., 2015). Ciliates are used as bioindicators in wastewater treatment plants (Chen, Xu, Tam, Cheung, & Shin, 2008; Foissner, 2016). In natural aquatic environments, the diversity and community composition of ciliates are influenced by e.g. salinity, pH, and anthropogenic pollution (e.g. Gong et al., 2015; Jiang, Xu, Hu, Warren, & Song, 2013). One of the main merits of assessing the diversity of protists as bioindicators is their documented rapid change to environmental conditions (Payne, 2013). The assessment of both meiofaunal and protist biodiversity is therefore a good proxy in monitoring programmes to study changes in ecosystems. However, these communities are rarely studied together and challenges still include being able to investigate multiple communities from bulk sediment samples without time consuming activities involved in studying the benthic meiofaunal fraction such as sieving, sorting, and microscopy.
In the last ten years, environmental DNA (eDNA) metabarcoding studies using the 18S rRNA marker gene have been extensively conducted to study microeukaryotes including nematodes, forams, and ciliates (Bik et al., 2012; Elias Broman et al., 2019; Carugati, Corinaldesi, Dell’Anno, & Danovaro, 2015; Fonseca et al., 2010; Forster et al., 2019; Lallias et al., 2014; Lara & Acosta-Mercado, 2012; Nascimento et al., 2019; J. Pawlowski, Lejzerowicz, & Esling, 2014; Peham, Steiner, Schlick-Steiner, & Arthofer, 2017). Such tools have drastically reduced the time needed to taxonomically classify organisms compared to morphological taxonomic techniques, that also involves sieving and sorting of organisms (Carugati et al., 2015). However, limitations exist with DNA metabarcoding such as non-optimized PCR protocols and primer bias when targeting multiple taxa (Kelly, Shelton, & Gallego, 2019). In addition to using eDNA that will assess the biodiversity of both living organisms plus non-degraded DNA from dead organisms, environmental RNA (eRNA) is targeting only living organisms or RNA derived from organisms of recent origin in the environment (Cristescu, 2019; S. A. Wood et al., 2020). With new bioinformatic tools that can taxonomically classify hundreds of millions of sequences within minutes to hours (e.g. D. E. Wood, Lu, & Langmead, 2019) and estimate relative abundances at species or genus level (e.g. Lu, Breitwieser, Thielen, & Salzberg, 2017), it is valuable to investigate how eRNA combined with the latest sequencing technology can be used to assess differences in biodiversity of active multiple communities from highly diverse and dense environments such as soils and sediments, while avoiding the biases associated with PCR amplification.
Here we assessed the biodiversity and community composition of three microeukaryotic groups in sediment samples: nematodes, forams, and ciliates, along an OC gradient in a coastal archipelago in the Gulf of Finland, Baltic Sea. The aim was to investigate if eRNA shotgun sequencing, without any sieving or sorting of samples (i.e. bulk sediment), could be used to detect differences in biodiversity of multiple microeukaryotic communities. Additionally, we assessed if changes in nematode functional ecology (feeding type) as a response to the OC gradient could be detected. We expected that nematode deposit feeders would have higher relative abundance in stations with higher OC. This approach was coupled to the latest sequencing platform (Illumina NovaSeq S4, yielding ~87 million sequences per sample), and new bioinformatic tools to estimate taxonomic classifications and relative abundances from data of this size (Kraken 2 + Bracken 2 combination). The Gulf of Finland is characterised by strong environmental gradients associated with eutrophication (Andersen et al., 2015; Villnäs et al., 2019). This contributes to spatially heterogenous benthic macro-communities in terms of diversity and composition in this ecosystem (Bonsdorff, Laine, Hanninen, Vuorinen, & Norkko, 2003). The Gulf of Finland is therefore a well-suited system to investigate if a similar heterogeneity exists in active microeukaryotic communities.