Salinity drives meiofaunal community structure dynamics across the Baltic ecosystem

Abstract Coastal benthic biodiversity is under increased pressure from climate change, eutrophication, hypoxia, and changes in salinity due to increase in river runoff. The Baltic Sea is a large brackish system characterized by steep environmental gradients that experiences all of the mentioned stressors. As such it provides an ideal model system for studying the impact of on‐going and future climate change on biodiversity and function of benthic ecosystems. Meiofauna (animals < 1 mm) are abundant in sediment and are still largely unexplored even though they are known to regulate organic matter degradation and nutrient cycling. In this study, benthic meiofaunal community structure was analysed along a salinity gradient in the Baltic Sea proper using high‐throughput sequencing. Our results demonstrate that areas with higher salinity have a higher biodiversity, and salinity is probably the main driver influencing meiofauna diversity and community composition. Furthermore, in the more diverse and saline environments a larger amount of nematode genera classified as predators prevailed, and meiofauna‐macrofauna associations were more prominent. These findings show that in the Baltic Sea, a decrease in salinity resulting from accelerated climate change will probably lead to decreased benthic biodiversity, and cause profound changes in benthic communities, with potential consequences for ecosystem stability, functions and services.


| INTRODUC TI ON
Biodiversity underpins essential ecosystem services for human benefits such as food availability, provision of clean water, recreational areas and activities affiliated with human health, and play key roles in ecosystem processes such as nutrient cycling and secondary production (Pan, Marcoval, Bazzini, Vallina, & Marco, 2013). Climate change, eutrophication with associated algal blooms, hypoxic bottom zones, and changes in salinity are contemporary major threats for coastal biodiversity (Pan et al., 2013). Such impacts need to be understood in order to predict how marine ecosystems will respond to future changes.
The Baltic Sea is a brackish water system that contains strong abiotic environmental gradients in salinity, depth and temperature that structure its biodiversity and benthic community structure (Ojaveer et al., 2010). The Baltic Sea is also affected by multiple anthropogenic pressures like eutrophication (Conley, 2012) and climate change (Vuorinen et al., 2015). In its deeper basins, below the halocline, hypoxic and anoxic benthic zones are widespread (Conley, 2012). Low-saline areas (<6 ppt) have expanded in the Baltic Sea since the 1970s and are predicted to further increase with climate change due to increased freshwater runoff and increased water column stratification (Vuorinen et al., 2015). The Baltic Sea therefore presents an ideal ecosystem to study the impact of future climate change scenarios on biodiversity (Ojaveer et al., 2010) and concomitant effects on benthic structure and consequent benthic-pelagic coupling (Griffiths et al., 2017). Most knowledge on how benthic organisms in the Baltic Sea react to these pressures are based on benthic macrofauna, while meiofauna (animals < 1 mm) have been studied much less. Meiofauna is a much more abundant and diverse metazoan group in sediments than macrofauna and plays an important role in a number of ecosystems process (Bonaglia, Nascimento, Bartoli, Klawonn, & Brüchert, 2014;Nascimento, Näslund, & Elmgren, 2012;Näslund, Nascimento, & Gunnarsson, 2010). However, there are still large knowledge gaps regarding how meiofaunal diversity and structure is affected by environmental changes (Bik et al., 2012). Recent DNA and RNA techniques now offer new possibilities to better address such questions on larger geographical scales than previously possible with traditional techniques.
Meiofauna have a short life span and are known to stimulate bacterial growth (reviewed in Coull & Chandler, 2001). Meiofaunal diversity and community composition are structured by several interacting factors; both abiotic and biotic (Giere, 2009). Oxygen is important for meiofaunal survival and metabolism (Braeckman, Vanaverbeke, Vincx, van Oevelen, & Soetaert, 2013), with some exceptions for facultative anaerobes with anaerobic mitochondria (Tielens, Rotte, van Hellemond, & Martin, 2002). Additionally, meiofaunal species richness and abundance have been found to increase with increasing salinity (Coull, 1988). In benthic environments these organisms rework sediment particles through e.g., bioturbation (Cullen, 1973), and have been found to affect porosity and increase the transport of solutes in the sediment (Aller & Aller, 1992).
Meiofauna utilize many sources of organic substrates in the lower trophic food web, e.g., bacteria, and detritus such as settling algal matter from the pelagic water (reviewed in Schratzberger & Ingels, 2018). Furthermore, they have also been found to stimulate degradation of sediment organic matter (OM) and bacterial denitrification (Bonaglia et al., 2014), and may therefore be key players in sediment habitats influencing carbon and nitrogen cycles.
One of the most diverse animal groups on Earth are the roundworms, i.e., nematodes (Zhang, 2013), and they are also one the most abundant meiofauna in sediments (Coull, 1999). Nematodes have been found to enhance the oxygen production in diatom biofilms (Mathieu, Leflaive, Ten-Hage, De Wit, & Buffan-Dubau, 2007), and to enhance the mineralization of OM (Nascimento et al., 2012).
Because of their different feeding behaviours in sediments, nematodes have been widely used in functional analyses (e.g., Semprucci, Cesaroni, Guidi, & Balsamo, 2018;Vanaverbeke, Merckx, Degraer, & Vincx, 2011). An increased knowledge of nematode community composition in the Baltic Sea could therefore further elucidate the role of trophic interactions in sediments under anthropogenic stress and climate change scenarios.
Benthic macrofauna have been observed to control meiofauna populations (or limit in some cases) through for example predation (Olafsson, 2003) and competition of limited resources (Ingels, Dashfield, Somerfield, Widdicombe, & Austen, 2014;Nascimento, Karlson, Näslund, & Elmgren, 2011;Olafsson, 2003). There has been extensive work, mainly laboratory or in situ experimental approaches, conducted on meiofauna-macrofauna interactions using morphological approaches (Olafsson, 2003). Such studies have yielded a variety of mixed results, but also a general consensus that macrofauna bioturbation structures the meiofauna community (Olafsson, 2003). These ecological interactions have been shown to have an importance on biogeochemical cycles; however, studies that focus on meiofauna-macrofauna interactions in situ and over regional and ecologically relevant scales are scarce. Macrofauna diversity is generally higher in more saline regions (Gogina et al., 2016), and meiofauna-macrofauna interactions might therefore be more prominent in saline regions with higher diversity and species richness. Gaining such insights will help to elucidate potential trophic interactions in the sediment and how these may be affected by contemporary ecological and environmental pressures.
Studies using metabarcoding, i.e., high-throughput sequencing of taxonomically-informative marker genes, to investigate meiofaunal biodiversity is a growing field (Bik et al., 2012;Carugati, Corinaldesi, Dell'Anno, & Danovaro, 2015;Fonseca et al., 2010;Lallias et al., 2014;Peham, Steiner, Schlick-Steiner, & Arthofer, 2017), and opportunities to facilitate such insights and the investigation of 18S rRNA gene meiofauna community in the Baltic Sea are now emerging (Nascimento, Lallias, Bik, & Creer, 2018). Compared to traditional morphological taxonomic techniques, modern sequencing tools facilitate the study of regional patterns of meiofauna diversity in less time while requiring no specific expertise in morphological taxonomy (Carugati et al., 2015). However, caveats do exist, such as not being able to determine absolute abundance and limitations of reference databases to assign taxonomy (Carugati et al., 2015). The benthic meiofauna community of the Baltic Sea is still largely unexplored although many benthic habitats in the Baltic Sea are under stress from anthropogenic pressure.
In this study we aimed to assess Baltic Sea meiofaunal diversity and community structure at the ecosystem level. An additional goal was to improve our understanding of possible future trajectories of benthic coastal diversity by using the Baltic Sea as a model system.
We specifically tested the following hypotheses: (a) salinity is an important driver of meiofauna community structure in the Baltic Sea, and (b) biotic interactions with macrofauna play a more important role in structuring meiofauna communities in more saline areas coincident with higher macrofaunal species richness. To test these hypotheses we sampled sediment along a salinity gradient in the central Baltic Sea (Baltic Proper). In order to identify changes in community composition and diversity of benthic taxa, a combination of traditional taxonomic assessment for macrofauna and metabarcoding DNA analyses for meiofauna were used. Meiofauna community composition was then analysed together with macrofauna community composition and sediment abiotic parameters (sediment water and OM content, bottom water temperature, salinity, and dissolved oxygen). Finally, because of the large relative abundance and diversity of nematodes, data for the phylum Nematoda were analysed separately to investigate their functional ecology (maturity index and feeding type) along the salinity gradient.

| Field sampling, collection of macrofauna, and abiotic variables measurements
Soft bottom sediment of similar clay-muddy habitats and water samples were collected in May-June 2015, at 44 stations in the Baltic Sea from the Stockholm region to the southern Arkona basin proper, during the yearly Swedish national and regional benthic monitoring program ( Figure 1). Benthic macrofauna communities were sampled with a van Veen sediment grab (0.1 m 2 ) from each station (typically one replicate per station, except for nine stations that had three replicates due to a yearly monitoring programme: 4, 5, 8, 11, 13, 14, 33, 37, and 44). All macrofauna abundance and biomass data were normalized for m 2 sediment. Benthic meiofauna and sediment variables were measured by collecting sediment cores from the 44 stations using a Kajak gravity corer (surface area: 50 cm 2 , one core per station). To investigate large spatial scale variation, we sampled more stations within each region, rather than performing repeat sampling within stations. The latter strategy has been demonstrated to be effective at capturing both small and large spatial scale diversity of European meiofaunal communities Lallias et al., 2014). Consequently, sediment collected from stations within the same region were treated as ecological replicates for further analyses. For the meiofauna and sediment organic matter the top 0-2 cm layer of each sediment core was sliced and homogenized into a clean and rinsed 215 ml polypropylene container (207.0215PP; Noax Laboratory). Sampling and slicing equipment was rinsed with deionized water between each sample. The sliced portion was then divided into: (a) 15 ml transferred to a 90 ml polypropylene container (207.0090PP; Noax Laboratory) for measurement of water and OM content, and (b) the remaining portion kept for meiofauna extraction. Samples were frozen at -20°C while on the boat, put on ice during transportation to the laboratory (~2 hr), and finally stored at -20°C until DNA extraction. Sediment collected for macrofauna was sieved through a 1 mm mesh and the animals retained in the sieve were transferred to 100-1,500 ml polypropylene containers (Noax Laboratory) and conserved in 4% buffered formaldehyde for three months (EN 16655:2014(EN 16655: , 1992. Macrofauna abundance and wet weight biomass were counted visually and weighed according to the European standard (EN 16665:2014(EN 16665: : 1992. Sediment water content (%) and OM content (%) were analysed according to Dybern, Ackefors, and Elmgren (1976). In more detail, determination of water content was conducted by drying sediment at 80°C to a constant weight (at least for 12 hr, typically overnight). The OM content was measured by reweighing the dry sediment after loss on ignition (500°C for 2 hr). Bottom water was sampled at each station, approximately 20 cm above the sediment surface, with a modified Niskin bottle. On deck temperature and salinity were measured in the collected bottom water using a digital multimeter (WTW Cond 340i), and dissolved oxygen (O 2 ) was measured in duplicate samples using the Winkler titration method (EN 25813:1992).

| Collection of meiofauna, DNA extraction, and sequencing
The sediment collected for meiofauna analysis was thawed at the laboratory and meiofauna were extracted from the sediment using the procedure described by Nascimento, Karlson, and Elmgren (2008). Sediment samples were sieved through a sterilized 40 µm sieve (autoclaved, rinsed with 90% ethanol and MilliQ water between

| Statistics
To detect differences in community composition between sites nonmetric multidimensional scaling (NMDS) ordination was performed by loading Metazoa sequence variants data into the r package phyloseq 1.24.2 (McMurdie & Holmes, 2013) using r 3.5.1 (R Core Team, 2013). In more detail, NMDS plots of Bray-Curtis dissimilarity, based on the sequence variants relative abundances and presence/absence (Sørensen), were constructed using the "ordination" and "plot.
ordination" functions in phyloseq. To test for statistical differences in community composition, this was followed by statistical testing with pairwise PERMANOVA tests (9,999 permutations) using the adonis function in the vegan package (Oksanen et al., 2018). In addition, the "betadisper" function in the vegan package was used to find differences in multivariate homogeneity of beta diversity variance between regions (Anderson, Ellingsen, & McArdle, 2006). This was followed by PERMANOVA tests of the homogeneity variance between regions, and plotted using ggplot2 package as the average distance to the centroid. Alpha diversity indexes (ACE, Chao1, and Shannon's H) were based on all Metazoa sequence variations counts and were calculated in the software explicet. Before alpha diversity analysis, counts were subsampled to 2,200 counts for each station (lowest sample size; Station 14), except for one station (Station 33 Gotland) that was excluded due to having fewer counts than the amount of metazoan sequence variants in the data set (station 33:291 counts). Afterwards the data set was bootstrapped 100 times, alpha diversity was calculated, and the mean of each alpha diversity index reported. In addition, ACE alpha diversity was also calculated by using nonsubsampled counts using the fossil 0.3.7 package (Vavrek, 2011) in r.
Based on classified nematode genera that could be annotated according to functional traits, (a) the maturity index described by Bongers (1990) was calculated to identify habitat colonizers or persisters (based on a 1-5 scale per genera; values closer to 1 indicate colonizers), and (b) feeding type was determined according to Wieser (1953) for each genera based on available literature outlining their buccal cavity morphology. Statistics on alpha diversity, taxonomic groups, and nematode feeding types were conducted in the software IBM spss Statistics 25. The normality distribution of the data was tested with Shapiro Wilk tests, and nonparametric Mann-Whitney U and Kruskal-Wallis tests were used on data not following a normal distribution.
The function "bioenv" in the r package vegan was used to test which, or combination of, abiotic variables (based on euclidean distances) had the highest rank correlation explaining the Bray-Curtis dissimilarity distribution of sequence variants among the sampling stations (with the following parameters: method = "spearman", index = "bray", partial = NULL, metric = c("euclidean")). This was followed by Mantel tests (Mantel, 1967) of Bray-Curtis dissimilarity distances and abiotic variables (salinity and spatial distance) in r using the ade4 package and 9,999 permutations (Dray & Dufour, 2007).
To find potential biotic interactions between meiofauna and macrofauna, co-occurrences among meiofauna, and possible community niches based on abiotic variables we conducted correlation network analysis (Röttjers & Faust, 2018). Correlation network analysis was conducted by importing Metazoa genera sequence counts as primary data, and the measured values for abiotic variables and macrofauna abundances per sediment per m 2 as metadata using conet 1.1.1 (Faust & Raes, 2016) and visualized in cytoscape 3.6.1 (Shannon et al., 2003). The setup in conet consisted of normalizing sequence counts as proportions per sample; setting spearman correlations with ρ thresholds ≤-0.7 or ≥0.7, and Fisher's z p-value threshold <.05 with Bonferroni adjustment for multiple-test correction. We are aware that our data set included a complicated setup, less commonly used in network software (Röttjers & Faust, 2018), with 18S rRNA gene sequencing data combined with both abiotic and macrofauna data. However, we applied a number of recommendations outlined in Röttjers and Faust (2018) to minimize potential limitations of such an approach, namely: (a) data from meiofauna were physically isolated from sediments; (b) the DADA2 methodology that incorporates denoising algorithms was employed; (c) metazoan sequence variants were grouped into 125 groups (120 genera and five unclassified groups); and (d) differences in meiofaunal community composition between north and south sample regions were based on the NMDS Bray-Curtis. In combination with the bioenv analysis that identified salinity as a major factor of diversity and community structure, we divided the data into two clusters (north and south Baltic proper) to remove influences of heterogeneous local environmental factors.
Such precautions strengthen the correlation network analysis, and emphasises ecological relevance (as reviewed in Röttjers & Faust, 2018).

| RE SULTS
The DADA2 analysis of the raw sequence data resulted in 3,309 18S rRNA gene sequence variants of which 770 belonged to the Metazoa kingdom distributed over 120 genera. On average 23% of the sequences per sample were unassigned with blast, and could not be classified to a phyla in the silva database, and were therefore not included in further analyses. See Table S1 for a list of all DADA2 sequence variants, the taxonomic classifications and sequence counts, and Table S2 for a full list of metazoan genera.

| Meiofauna beta and alpha diversity
The  Table S3 for a full list of p-values for the geographic regions). In addition, the two regions in the SBP were significantly different from each other (i.e., Bornholm being lower compared to Arkona; betadisper, PERMANOVA, p < .01; Figure 2b). There was a relatively large abundance of pelagic Arthropoda in the 18S rRNA gene data set, and therefore, NMDS analysis was also performed without these sequence variants (mainly pelagic Copepod genera Eurytemora and Temora; see Table S3 for a full list of excluded genera). This analysis also showed two distinct clusters between the NBP and SBP ( Figure S2a;  Figure S2b and Table   S3). As such, the differences in meiofaunal homogeneity variance between regions were larger after the pelagic Arthropoda had been excluded from the data set.
A higher alpha diversity, based on all Metazoa sequence variants, was observed in the SBP stations compared to the NBP (p < .01 for all indexes (ACE, Chao1, and Shannon's H); One-way ANOVA; Figure 3).
When alpha-diversity was tested on the Nematoda sequence variants alone, there was also a significant difference (p < .01 for ACE and Chao1, F = 4.1 for both; Shannon's H not significant; Figure 3).
Similar results for the nematodes were also observed when ACE was tested on non-subsampled data (p < .01), although not when all metazoa sequence variants were tested (p = .08). These results show that a higher diversity of Metazoa sequence variants were obtained in SBP sediments. A full list of alpha-diversity indexes for each station for all meiofauna and Nematoda sequence variants is available in Table S4.

| Meiofauna community composition
Similar to the NMDS and alpha diversity analysis, there was a difference in relative abundance in phyla between the NBP and SBP, with Arthropoda having a higher relative abundance in the NBP compared to the SBP (p < .01, Mann-Whitney U test). In contrast, the phylum Nematoda had a lower relative abundance in the NBP (p < .01, Mann-Whitney U test; Figure 4a). Looking closer at the genera belonging to Arthropoda, the genus Eurytemora was dominant in the NBP compared to the SBP where Temora had the highest relative abundance (p < .01 for both, Mann-Whitney U tests; Figure 4b).
Nematodes showed a much higher diversity compared to the other major phyla, with 60 Nematoda genera compared to 28 and 19 genera belonging to Arthropoda and Platyhelminthes, respectively ( Figure 4c, a full list of all genera is available in Table S2). The phylogenetic placement of Nematoda sequences on a reference tree showed that the most dominant Nematoda sequences (Table   S1) aligned closely to NCBI reference sequences from Holovachov (Nematoda results are available in Figure S4 and  Figure 4c). In contrast, the genera Axonolaimus and Enoplolaimus had a higher relative abundance in the SBP (p < .05 and p < .01, respectively; Kruskal-Wallis test; Figure 4c). In addition, the relative abundance of unclassified sequence variants belonging to the Nematoda phylum was higher  (Table S1 and Figure S3).
Maturity index calculations, used to estimate nematode genera as habitat colonizers or persisters, showed that all observed nematode genera in the current study are classified closer to colonizers rather than persisters (maturity index < 2.7; to five indicate persisters with slow reproduction (Bongers, 1990).
Nematode genera were also classified into feeding type (according to Wieser, 1953), and showed that the most southern region Arkona had more predators/omnivores compared to all other regions (oneway ANOVA Tukey HSD post hoc test, p < .01; Figure 5d). Looking at the feeding types of nematode genera with a high relative abundance in the NBP the Cyatholaimus and unclassified sequence variants (potentially Chromadorita) were classified as epistrate feeders (feeding type 2A) (Table S5; unclassified sequence variants not included). In the SBP the genera Enoplolaimus was classified as predatory possessing large teeth (2B), while Microlaimus was classified as 2A (Table   S5). Other genera with a high relative abundance in the Nematoda data set such as Aphanolaimus, Daptonema, and Axonolaimus were classified as type 1A or 1B, being either selective or nonselective deposit feeders, respectively. A full list of maturity indexes and feeding type classifications is available in Table S5.
Looking at the Platyhelminthes the genus Odontorhynchus showed a significant difference with a higher relative abundance in the SBP, although with high variation (p < .05, Mann-Whitney U test; Figure 4d). In the two SBP regions the genus Placorhynchus was dominant in the Bornholm region while Odontorhynchus was more prevalent in the Arkona region (p < .05, Mann-Whitney U test).

| Macrofauna in the sediment
The Macrofauna data showed a higher species richness in the SBP than in NBP (on average eight species per station compared to four in the NBP; Figure 6). There were also more species belonging to the Annelida phylum in the SBP, e.g., Bylgides sarsi, Nepthys caeca, Pygospio elegans, and Scoloplos armiger ( Figure 6). The Bornholm re-  Table S6.

| Abiotic variables
Bottom water salinity increased as expected in the Baltic Sea (Ojaveer et al., 2010), from 5.3 ppt salinity in the NBP to 18.9 ppt in the SBP (Table 1). Bottom water temperature was generally low for most stations (average of ~6°C) except a few stations in the Östergötland region that had temperatures >10°C (average of ~11°C, stations 16-20;

| Correlations of meiofauna with abiotic variables and macrofauna data
Abiotic data from all stations were tested with Bray-Curtis dissimilarity of sequence variants, and the best explainable abiotic variables were longitude, latitude and salinity (ρ = 0.73). Mantel tests also confirmed that these abiotic variables were significantly correlated with the beta diversity measures (R 2 = 0.67 and p < .01, for both salinity and spatial location tested). The combination of abiotic variables latitude, sediment water content, and oxygen had the best rank correlation explaining the beta diversity among the stations in the SBP (ρ = 0.57; "bioenv" test in r package vegan). This was in contrast to the NBP where longitude, water depth, and oxygen were the best explainable variables (although with a low rank correlation, ρ = 0.32; in accordance to the lack of correlations with abiotic factors in the correlation network; Figure 7a).
Correlation network analysis was conducted on the NBP and SBP separately because the NDMS Bray-Curtis indicated differences in meiofaunal community structure between the two areas.
In addition, the bioenv analysis showed salinity to be a strong driver influencing meiofaunal community structure and diversity in the Baltic Proper. This precaution of removing sample heterogeneity in a larger ecosystem-wide data set is in accordance to Röttjers and Faust (2018)

| Abiotic explanatory variables of meiofaunal diversity
Salinity was the major explanatory variable of benthic meiofauna community composition in the Baltic Proper. In addition to our findings, salinity has been observed to influence macrofauna in the Baltic Sea (Gogina et al., 2016) and meiofauna community structure studied elsewhere (Coull, 1988;Lallias et al., 2014). Interestingly, similar findings were also recently discovered for sediment bacteria community composition along a salinity transect in the Baltic Sea (Klier, Dellwig, Leipe, Jürgens, & Herlemann, 2018). Bottom water oxygen also correlated with the difference in meiofaunal community composition, especially in the SBP. The role of oxygen is not surprising considering that oxygen is essential for the majority of meiofaunal organisms (Braeckman et al., 2013), and oxygen availability is known to cause shifts in the community composition of for example nematodes (Nguyen et al., 2018). The local regions as defined in this study ( Figure 1 and Table 1), also harboured significantly different communities of meiofauna (Figures 2 and 4). This difference could be attributable to specific salinity preferences, but also due to the sediment substrate and available food resources (Lee, Tietjen, Mastropaolo, & Rubin, 1977), and adult dispersal through water currents (Hagerman & Rieger, 1981). Marine meiofaunal communities have previously been shown to be heterogeneous both at large  and small spatial scales (Findlay, 1981). Our results indicate salinity to be a major barrier to dispersion of meiofauna species in the Baltic soft sediment, by limiting the dispersion of marine species to the north and of freshwater species to the south. Limitation to dispersion is an important factor driving community assembly in ecological systems (Vellend, 2010). Therefore the salinity gradient in the Baltic Sea influences sediment habitats with different kinds of food and predators that will in turn influence the meiofauna community composition and diversity.

| Geographical differences in community composition
Meiofaunal diversity was dominated by a large variety of Nematoda genera. This was not surprising considering that nematodes are highly diverse (Zhang, 2013), and typically the most abundant  (Coull, 1999). The SBP had a different Nematoda community composition, probably due to the higher salinity conditions that have previously been found to influence diversity and community structure in the Baltic Sea (Ojaveer et al., 2010). In the NBP there was a large proportion of unclassified Nematoda sequences (Figure 4c) and could possibly be due to the lack of freshwater-brackish species being classified in the reference databases (Holovachov et al., 2017). These unclassified sequences were indicated in the phylogenetic placement analysis to be affiliated with the genus Chromadorita ( Figure S3). This genus has previously been found in the Baltic Sea (Jensen, 1979) and contains species living on macrophytes (Jensen, 1979), free-living and feeding on diatoms (Jensen, 1984), and living inside cyanobacterial biofilms (Gaudes, Sabater, Vilalta, & Muñoz, 2006). The most southern region Arkona had not only a higher diversity but also a higher proportion of nematode predators/omnivores (Figure 5d), which could explain why there were more ecological correlations in the SBP (Figure 7b).
Even though the PCR primers used might have selected for certain eukaryotic species, and primer bias are likely to pervade all metabarcoding studies, we used the same biodiversity discovery method (i.e., metabarcoding primers) throughout. Despite the imperfect nature of metabarcoding (and other ecological sampling approaches), the difference among regions and areas in the Baltic Proper were statistically significant and showed stark dissimilarities in community composition.
In addition to nematodes, there was also a large relative abundance of Arthropods in the 18S rRNA gene data set, especially in the NBP. The majority of the Arthropoda belonged to the pelagic copepod genera Eurytemora and Temora in the NBP and SBP, respectively (Figure 4a,b). The hatching rate and development time of e.g., Eurytemora affinis is negatively affected by low salinity (Karlsson, Puiac, & Winder, 2018) which can explain the difference between the north and south regions. Possible additional explanations for copepods being in the sediment could be due to sinking marine snow containing carcasses, resting stages such as buried eggs or dormancy (Dahms, 1995). The high relative abundance of Arthropoda could therefore be derived from DNA being extracted from large amounts of copepod eggs or resting stages buried in the sediment surface. Considering that similar results have also been observed by Nascimento et al. (2018) from sediments collected in the Stockholm region, the large proportion of copepods is probably a trait for low saline waters (<10 ppt) in the Baltic Sea. Compared to the SBP where salinity is higher, the availability of copepod eggs in the low-saline NBP can be a larger source of energy for benthic macrofauna populations (Karlson & Viitasalo-Frösen, 2009). In addition, because the hatching rate is slower in low salinity (Karlsson et al., 2018) the accumulation of a seed bank followed by subsequent hatching could enhance the benthic-pelagic coupling. Our results highlight important geographic differences in meiofaunal communities that are only possible to uncover with modern molecular tools (Fonseca et al., 2010).

| Biotic interactions
Macrofauna species richness and meiofauna diversity were both higher in the SBP (Figures 6 and 3, respectively). Nascimento et al. (2011) found that a higher species richness of macrofauna increased interference competition among meiofauna and/or limited food availability in a laboratory study. Potentially, this could partly explain why there were more ecological connections between macro-and meiofauna in the SBP as indicated by the correlations network data ( Figure 7b). On the other hand, macrofaunal bioturbation can create more habitable niches and higher variety of food types allowing for a higher meiofauna diversity (Meysman, Middelburg, & Heip, 2006). The significant correlations included mainly Annelida as well as crustacean macrofauna, which are well-known bioturbators (Krantzberg, 1985). In addition, bottom water oxygen was one of the central nodes in the correlation network with connections to meio-and macrofauna ( Figure 7b). It is therefore possible that oxygen rich burrows made by annelids (Aller, 1988) or other modes of bioturbation by macrofaunal organisms (Krantzberg, 1985) stimulate bacterial growth and make specific niches and habitats favourable for meiofauna (reviewed in Olafsson, 2003). However, negative macro-meiofauna interactions have also been previously reported (reviewed in Olafsson, 2003). High macrofauna diversity can increase sediment oxygen consumption (Bolam, Fernandes, & Huxham, 2002), and interference competition with meiofauna by limiting its access to freshly deposited detritus (Nascimento et al., 2011). Such mechanism could explain some of the negative or nonsignificant correlations between macro-and meiofauna taxa found in our study. For example, we observed several Mollusca macrofauna in correlation clusters with meiofauna genera in the SBP (Figure 7b). However, this kind of interaction was not as prominent in the NBP. This is in accordance with previous experimental studies with sediments from the Sörmland region amended with bivalve M. balthica that showed no significant difference on the majority of meiofauna, including nematodes (Olafsson, Elmgren, & Papakosta, 1993). Considering that correlation network analysis can be a major strength to visualize and detect specific habitat niches (Röttjers & Faust, 2018), the meiofauna-macrofauna associations observed here could be indirect effects of shared niche preference. In addition, predation is an important mechanism structuring diversity in more stable and F I G U R E 7 Correlation networks of spearman correlations based on data from north (a) and south Baltic Proper (b). The correlations included meiofauna 18S rRNA gene data (each node represents one Metazoa genus), abiotic variables, and macrofauna abundance data. The mean was used for the two oxygen technical replicates. The colour of the lines denote ρ ≥ 0.7 (red) or ≤-0.7 (blue). All correlations are statically significant (p < .05). All abiotic nodes have been labelled as well as a few genera/macrofauna nodes according to the results presented in the text. Nodes with black borders denote unclassified sequences belonging to a certain phylum [Colour figure can be viewed at wileyonlinelibrary.com] tropically complex communities (Menge & Sutherland, 1976). The NBP had lower diversity and has a history of being more affected by eutrophication compared to the southern region Arkona (Andersen et al., 2015). The higher relative abundance of nematode predators in the SBP (Figure 5d) could indicate a relatively more stable environment where predation can maintain a higher diversity helped by more macrofauna-mediated niches, biodiversity and interactions.
Although, network correlations based on metabarcoding data need to be treated with caution (see Röttjers & Faust, 2018), our results clearly indicate that there are fewer, direct or indirect associations between meiofauna and macrofauna in low-saline areas in the Baltic Sea.

| Effects of climate change and future scenarios
The area of low saline regions in the Baltic Sea (surface water salinity < 6 ppt) has increased since the 1970s and are predicted to further increase with climate change due to elevated levels of runoff (Vuorinen et al., 2015). As indicated here a decrease in salinity might be accompanied by a decrease in meiofaunal biodiversity and biotic interactions in the Baltic Sea. Salinity strongly influences the community composition and diversity in other coastal systems (Lallias et al., 2014;Van Diggelen & Montagna, 2016) where similar effects can happen if salinity is reduced as a consequence of climate change.
Additionally, it is clear from our results that a continued expansion of hypoxic bottom zones will significantly alter benthic community structure. This may influence important ecosystem functions regulated by meiofauna, like OM degradation and nutrient cycling. Here, we show that multiple anthropogenic pressures like eutrophication (Finni, Kononen, Olsonen, & Wallström, 2001), expansion of hypoxic bottom zones (Meier et al., 2011), and of low-salinity areas (Vuorinen et al., 2015), will probably have profound impacts on benthic communities of anthropogenically stressed coastal systems. Ongoing environmental change will lead to lower benthic biodiversity and fewer biotic interactions. Such structural changes to benthic community composition will probably influence ecosystem functions and services, and decrease ecosystem stability (McCann, 2000).

ACK N OWLED G EM ENTS
The authors acknowledge support from the National Genomics

AUTH O R CO NTR I B UTI O N S
E.B. analysed data and drafted the manuscript, C.R. sampled in the field and together with C.S. conducted laboratory work. J.G. helped with field sampling and gave feedback on the manuscript, S.C. helped designed the study and gave feedback on the manuscript. F.N. designed the study, conducted laboratory work, analysed data and contributed to the manuscript writing. All authors gave final approval for publication.

DATA AVA I L A B I L I T Y S TAT E M E N T
The raw sequence data have been uploaded and are available on the NCBI database with the following BioProject number: PRJNA497177.