Trends in records and contribution of non-indigenous and cryptogenic species to marine communities in Danish waters: potential indictors for assessing impact

We assess the trends and influences of non-indigenous and cryptogenic species (hereafter simply referred to as “NIS”) on Danish marine community compositions using three decades of quantitative monitoring data. Since the initiation of the Danish marine monitoring programmes in the 1980s, the number of marine NIS recorded in Denmark increased from 30 to 77 in 2014. Thus, of the total 77 marine NIS known from Denmark, 56 (73%) were captured in the standardized monitoring program, while the remaining 21 species (27%) were not detected (in particular gelatinous zooplankton, shallow water fish, parasitic invertebrates, and littoral angiosperms) because of limited spatial-temporal sampling efforts as well as methodological limitations. Significant exponential increases in records of non-native phytoplankton, benthic invertebrates and macroalgae (only in one region) relative to the total species records in the database were observed. Multivariate analyses of presence-absence data, indicated that the contribution of NIS to total community similarity increased over time, highlighting that NIS is becoming an increasingly important component of Danish marine communities. While the presence of NIS generally explained less than 10% of long-term changes in the community similarity across large regions, NIS presence showed local dominance. A correlation analysis indicated that changes in overall species composition within functional groups (phytoplankton, zooplankton, macroalgae and benthic fauna) were related to changes in NIS contribution but more strongly influenced by salinity, highlighting a well-described general positive relationship between salinity and species richness in Danish waters. The process described within this study could form the basis for the analysis of impact of NIS in marine water as required by the Marine Strategy Framework Directive (MSFD).


Introduction
Human activities have over centuries broken down natural dispersal barriers, that in the past separated biota from different oceans, seas, rivers and terrestrial environments (Thomsen et al. 2015;Caplat et al. 2016).
Growth in global transport and deliberate introductions of new species have also facilitated the spread of non-indigenous species (NIS = alien, exotic, non-native, introduced) (IUCN 2000; Thomsen et al. 2008a). Many NIS have caused ecological, economic and public health impacts globally (Ruiz et al. 1997;Ojaveer and Kotta 2015;Thomsen et al. 2016) and/or changed the structure and dynamics of ecosystems (Vilà et al. 2011), although the direct causality is often uncertain (Didham et al. 2005). Economic impacts range from financial losses in fisheries to expenses to industries for cleaning flow-pipes and structures from biofouling (Black 2001;Williams et al. 2010) whereas public health impacts may arise from the introduction of pathogens or toxic organisms (Thomsen et al. 2015;Turner et al. 2018). Finally, ecological impacts include changes to habitats, community structures, food web function, and in extreme cases, loss of native species (Galil 2007).
Studies of the ecological impacts of marine NIS are relatively scarce. In the Baltic Sea, 132 NIS and cryptogenic species have been reported with clear differences between regions (Ojaveer et al. 2017). For this region, impact assessment only exists for 43 NIS, mainly covering the eastern parts of the Baltic Sea (Zaiko et al. 2011). This study like other reports on NIS introductions include cryptogenic species, of unknown origin which is often the case, as verification of a species origin is difficult and may require genetic analysis (Gollasch et al. 2009). As cryptogenic species may constitute a significant proportion of the newly arrived species, these are often therefore often considered in assessment of NIS (e.g. Thomsen et al. 2009;Ojaveer et al. 2017;Zaiko et al. 2011). While there is no doubt that invasive species have modified marine ecosystems, evidence for most of the reported impacts is weak, as it is based on expert judgement or dubious correlations, with limited documentation from controlled or natural experiments . Improved management of NIS accordingly requires better evidence on marine biological invasions, which goes beyond a simply tally of NIS, by providing information on their distribution and impact. The most important and simple proxy for a NIS' success and impact is its abundance (Catford et al. 2009;Parker et al. 1999; Thomsen et al. 2011). An important first step to evaluate NIS success and impact can, therefore, be to document the relative abundances of NIS compared to native species, as well as temporal changes in relative abundances, in different regions.
A recent account estimated that 787 non-indigenous taxa have been found in European marine systems (Tsiamis et al. 2019) several of which require special attention as they are highly invasive with strong impact on marine ecosystem services and biodiversity, often causing adverse effects on environmental quality (Wallentinus and Nyberg 2007;Ojaveer et al. 2017). Due to the threats NIS pose in European seas, a series of legislative measures, including the Marine Strategy Framework Directive (MSFD;EU 2008EU , 2010EU , 2017 has been developed, with the aim to achieve "good ecological/environmental status" (Cardoso and Free 2008;EU 2008). In the context of the MSFD, NIS are treated as a distinct Descriptor (D2) of GES "Non-indigenous species introduced by human activities are at levels that do not adversely alter the ecosystem". Three criteria have been defined by which GES is to be assessed (CD 2017/848), the primary (D2C1) examines the number of NIS introduced per assessment period (6 years) per region, the secondary (D2C2 and D2C3) use abundance and spatial distribution of invasive NIS, and their impact on indigenous species communities and groups in addition to broad habitat types (EU 2017). Within the EU, and Regional Sea Conventions (RSC) such as OSPAR and HELCOM, there is ongoing work to develop and refine indicators based on these defined criteria to assess GES. The aim of the approach is to be able to assess GES at different geographical scales from areas (such as single country) to the regional (such as that of an RSC). It is suggested that monitoring for NIS is incorporated into broader biodiversity monitoring (CD 2017/848). Given the broad range and variability of monitoring methods used within a single country's biodiversity monitoring programme, in addition to differences in approach between countries, estimating abundance for a broad range of NIS in a reliable and comparative manner is difficult. Therefore, developing approaches for assessing GES using the most reliable and commonly available form of data (in this case absence/presence; P/A) will facilitate robust and comparative assessments to be conducted.
The objective of this study is to use Danish national monitoring data to (1) make a first long-term (three decades) assessment of changes in the relative contribution and influence of NIS on indigenous species and (2) correlate changes in NIS numbers to environmental conditions to identify conditions favourable for NIS establishment. Our analysis covers all marine functional groups across large-scale Danish marine regions, and make assessment based only on P/A data in order to make our findings relevant for future assessments of NIS in the framework of descriptor D2 under the MSFD. The assessment presented herein, provides a simple yet effective process on which indicators to assess the impact of NIS on a regional scale could be based.

Study area
The Danish waters were divided into 5 regions: Region 1 = North Sea and the Skagerrak, Region 2 = Kattegat, Region 3 = Limfjorden, Region 4 = the Belt Sea, and Region 5 = Western Baltic Sea (Figure 1). These regions roughly correspond to the HELCOM (2-5)/OSPAR (1-3) boundaries except for the border between the Kattegat and the western Baltic Sea, where we separated the Kattegat and Sound basins by a more southerly line to separate hydrographical regions and key environmental variables (Table 1).

Data sources, NIS (including cryptogenic species) in Danish waters
Our analysis combines historical records from experts and literature with long-term monitoring data  from fjords, estuaries and coastal and open water sites scattered across Denmark. These sampling sites have been regularly monitored as part of the Danish National Aquatic Monitoring and Assessment Programme (NOVANA, but also referred to as DNAMAP -see Riemann et al. 2016) and the annual fish survey programme (ICES 2012(ICES , 2014. To reduce data-variability that follows from changes to sampling protocols and sampling effort (i.e., the number of stations and samples varied between years and regions) we restricted our quantitative analysis to data collected since the nationwide monitoring programmes were fully establishment in 1989. We also extracted environmental data to correlate changes in NIS community contribution to changes in water quality. We based our analysis of NIS and cryptogenic species on the species listed in the European Alien Species Information Network (EASIN) recently published by Tsiamis et al. (2019). Furthermore, we consulted the AquaNIS (2018)

Sampling protocols
Phytoplankton. Quantitative samples from the upper mixed layer were collected and immediately fixated with Lugol (final concentration 2-5%), from 23 stations visited every two to four weeks year round. From these samples, species were identified and quantified using an inverted microscope according to a standard protocol following largely CEN (2011CEN ( , 2015. Data were standardized to individuals per litre. A total of 23 stations were analysed during the 1989 to 2014 period, represented by 3, 8, 3, 5 and 4 stations, respectively, in regions 1 to 5. A total of 561 taxa were recorded in the entire data set, most identified to species level, but in some cases to genus (e.g., "Gymnodinium sp."). Details of the Danish phytoplankton monitoring program can be found in Fossing et al. (2019). Zooplankton. Samples were collected by slowly withdrawing a pumping plankton pump from 25 m depth towards the surface (or, in shallower waters, from the bottom to the surface). The pumped water were collected on a screen, sieved though a 60 μm mesh and finally preserved in Lugol's (Final conc. 2-5%) until analysis (Fossing et al. 2019). Species were identified and counted using a dissection microscope following a standard protocol (Anonymous 2014). Data were standardized to individuals per litre. A total of 35 stations were analysed during the 1989 to 2014 period with 13, 2, 7 and 13 stations respectively in regions 2 to 5. The stations were typically visited six times each year. Note that some stations were closed whereas others are opened since the initiation of the program. Currently 6 stations are open.
Benthic invertebrates. Benthic invertebrates were covered by merging two data sets: 1) data covering hard bottom habitats for 22 stations in the open parts of the inner Danish waters and 2) a larger data set covering Danish soft bottom habitats in fjords and coastal regions. Data from the open waters are obtained by taxonomically skilled divers and verified with samples in the laboratory (Dahl and Lundsteen 2018). The open water data are balanced with respect to sampling season and spatial coverage and cover a fairly homogeneous habitat. By comparison, data from the Danish fjords and coastal areas are more heterogeneous. However, the community analysis performed on the merged data set was based on presence/absence data and thereby a less biased comparison of community structure across the data sets. Samples on soft bottom habitats were collected with grabs from 258 stations most of which were visited annually following standard sampling and processing procedures (Hansen and Josefson 2014). Soft bottom fauna was recorded as species specific abundances of individuals > 1 mm. For the period 1989 to 2014, 23, 38, 27, 149 and 21 stations were represented by in region 1 to 5, respectively. Details of the Danish benthic invertebrate monitoring program can be found in Hansen and Josefson (2014) and Dahl and Lundsteen (2018).
Macroalgae. Species-specific data were recorded in situ by divers in different depth intervals (Høgslund et al. 2013). A total of 898 sampling stations were monitored from 1989 to 2014 with a total of 29, 151, 68, 563 and 99 stations, respectively, in regions 1 to 5. Most stations in region 1 were located in The Wadden Sea area. Details of the Danish macroalgae monitoring program can be found in Høgslund et al. (2013).
Fish. Data were provided by DTU Aqua based on six annual surveys. Two surveys per year were conducted in the North Sea (first and third quarter with R/V DANA), two surveys (per year) were conducted in the Kattegat (first and fourth quarter with R/V Havfisken) and two surveys (per year) were conducted in the eastern Baltic Sea during (first and fourth quarter with DANA). The surveyed areas in the Baltic, Kattegat and North Sea mainly covered depth strata from 20-120 metres where most commercially important fish species are found. The surveys used slightly different trawls, but all of the trawls are designed to target juvenile benthic fish. The Kattegat and Baltic surveys use TV3-trawls with a 20 mm cod-end, but scaled in two different sizes to fit the different vessels. The survey conducted in the North Sea used a GOV bottom trawl. However, the mesh size is similar to 20 mm in the cod-end. Fish abundance for each species is recorded as a number of fish per trawled hour. In this study, however, we only utilized information on presence and absence of species. Data span the period from 1991 to 2014. A total of 3661 trawl stations were analysed with 542, 1036, 615 and 1468 stations, respectively, in regions 1, 2, 4 and 5. These surveys do not cover region 3 = Limfjorden. Details of the fish monitoring program can be found in Storr-Paulsen et al. (2012).
Environmental data. Data for total nitrogen (TN), chlorophyll a (Chl), primary production (PP), Secchi depth, and water column values of temperature and salinity were obtained from the national Danish marine database (http://www.oda.dk/). Sampling and analysis of these parameters are described in Riemann et al. (2016). A total of 355 stations were included in the environmental data analysis with 60, 54, 35, 134 and 50 stations, respectively, in regions 1 to 5.

Data analysis
Time series analysis. To investigate long-term trends in the number of NIS and their relative contribution, we calculated yearly and regionally weighted estimates of species presence/absence. Biotic data were matched to corresponding data on nutrients, temperature, salinity, Chl, primary production, and Secchi depth. For species specific data, the analysis a 3-step procedure: 1) selection of species data from databases, 2) quality control (e.g., checking names, coordinates, station numbers, outliers), and 3) estimation of spatially aggregated data across regions for each year. Due to unbalanced sampling in time and space, we used a general linear model (GLM) to estimate the marginal means of each variable in each region on an annual scale. The GLM approach was used on both species composition data and water quality monitoring data to produce annual time series within each region. Simple linear regression analysis was finally applied to investigate long-term trends in the number of NIS and their relative number as compared to the total number of species.
Community analysis. To investigate the importance of NIS on temporal trends in community composition, we first quantified multivariate community composition across the entire data set, based on Bray-Curtis similarity index of presence/absence data (grouped by region and sampling year). A presence/absence transformation was chosen to reduce problems associated with changes in sampling effort and spatial distribution over time. The similarity among species communities was also quantified using the Bray-Curtis similarity index in the "RESEMBLANCE" routing. Also, the contribution of NIS to % of the total community similarity within regions and periods was analysed with "SIMPER". All community analysis was carried out in the software package PRIMER (Clarke and Gorley 2015).
Correlating NIS contribution to environmental data. We used Spearman Rank correlation analysis to investigate relationships between the relative importance of NIS to native species numbers (%NIS, based on presenceabsence data) and environmental variables (total nitrogen, chlorophyll a, primary production, Secchi depth, water temperature and salinity) expected to influence the distribution of marine species. To enable comparison among regions, the temporal variation for each variable was first normalised to a single mean value in each region over the entire sampling period (1989 to 2014). To reduce the probability of analysing Once the observation of a species has been confirmed, it is assumed that the species remains part of the Danish biota regardless of lack of observations in some of the following years. spurious correlations, for example caused by unmeasured environmental or co-varying drivers, the analysis was performed on detrended time data. Detrending was done by using linear regression on each variable as a function of years. The regression coefficients (slope and intercept) were then used to estimate each environmental variable as a function of years, and these estimates were then subtracted from the original values. As a result, the detrended variables did not correlate with years allowing us to exclude temporal effect from relationships between NIS community contribution and environmental conditions. Finally, the routine BEST (Clarke and Gorley 2015) was used to identify environmental factors with strongest correlation to phytoplankton, zooplankton, macroalgae and benthic invertebrate communities. For each possible combination of environmental factors, the BEST function calculates a dissimilarity matrix based on normalized Euclidean distances. Correlations between the biotic and environmental matrices were calculated using the Spearman rank coefficient.

Overall changes in NIS numbers and contribution in Danish waters
Since the establishment of the full national monitoring programme in 1989, the number of NIS recorded by this process almost doubled from 32 to 56 in 2014 ( Figure 2 and Table 2). The linear regression of accumulated NIS had a slope of 0.90 (NIS accumulated = 0.90 × Year − 1760, p-slope < 0.001), indicating approximately 9 new marine NIS have been introduced to Danish waters per decade since the late 1980s. Most of this increase was attributed to new phytoplankton species (slope = 0.43), followed by benthic invertebrates (slope = 0.22), macroalgae (slope = 0.11), fish (slope = 0.06) and zooplankton (slope = 0.06). In addition to the accumulated number of   NIS detected in the NOVANA nationwide monitoring programme, at least 21 more marine NIS have been recorded in Danish waters through a range of other means (Supplementary material Table S1) providing a total accumulated number of 77 NIS in 2014. Out of these NIS, phytoplankton were the dominant group and accounted for 40%, followed by benthic invertebrates (29%), macroalgae (16%), fish (5%), parasites (5%), zooplankton (4%) and flowering plants (1%). The 21 NIS recorded outside of the monitoring programmes were represented by one angiosperm, two macroalgae, five phytoplankton, one zooplankton, six benthic invertebrates, four invertebrate parasites and two fish species. Six of these additional NIS were recorded before 1989 while the remaining 16 were observed later, demonstrating that the national monitoring program does not capture all NIS. Although 77 NIS have been observed in Danish waters, the actual number of NIS recorded each year in the monitoring programme ranges between 14 (Baltic Sea) and 31 (Limfjorden), highlighting that some NIS are not recorded consistently every year after their first observation. Furthermore, the monitoring data showed an uneven distribution of the total accumulated number of observed NIS among regions with more NIS in the transition waters (Kattegat, Limfjorden, and the Belt Sea) between the North Sea and the Baltic Sea (Figure 3).

Temporal trends in the number and relative abundance of marine NIS in five Danish regions
Phytoplankton. A total of 162, 232, 218, 239 and 195 species were recorded in regions 1 to 5 respectively since 1989 through 2014. Of these, 26 phytoplankton species were identified as NIS after analysis of the entire monitoring dataset (Table 2) and an additional 5 NIS were documented through other activities (Table S1). 18 species were recorded before 1989 (Tables 2 and S1). Of the ten highest ranking NIS, six are dinoflagellates, three flagellates and one diatom. Importantly, six out of the ten most frequently observed species are also classified as species that can form "harmful algal blooms" (HABs). The increase in phytoplankton NIS at the national scale (Figure 3) is also observed in the regional analysis of the phytoplankton communities ( Figure 4A and 4B). The highest average number of NIS (eight) as well as the highest relative occurrence of % NIS (12.4%) pooled across all sampling years was observed in the high saline North Sea/Skagerrak region. The observed increase in %NIS was 0.14 to 0.21% per year (slopes of regression lines) across all regions ( Figure 4B). More specifically, the highest rates of increase in % NIS were observed in the Kattegat and Belt Sea regions with 0.21 and 0.20% per year, respectively.
Zooplankton. A total of 250 taxa (most recorded to species levels) were recorded in the entire data set, with 212, 103, 150 and 155 recorded species from regions 2 to 5. Only two zooplankton NIS were observed in the entire data set since 1989 (Table 2 and Figure 4A); the calanoid copepod Acartia tonsa and the marine daphnia Penilia avirostris. One additional zooplankton NIS (Mnemiopsis leidyi) was recorded outside the monitoring programme (Table S1). None of the temporal trends for zooplankton was significant ( Figure 5B). Benthic invertebrates. A total of 1304 species were recorded from 1989 to 2014, with 525, 853, 465, 939 and 452 species from regions 1 to 5 respectively. A total of 16 benthic invertebrates NIS were observed in the entire data set (Table 2) and an additional six species were recorded outside the monitoring programme (Table S1). The occurrence of benthic invertebrate NIS varied between years and was highest in Limfjorden and lowest in the Baltic Sea region ( Figure 6A). The relative contribution of NIS to the total species number (% NIS) increased significantly over time in all five regions ( Figure 6B). The largest current (2014) NIS contribution was recorded in Limfjorden (10%), followed by the western Baltic Sea (6%) with the other three regions having 4 to 5% NIS. The observed increase in % NIS was approximately 0.3% per year (slopes of regression lines) for all regions ( Figure 4B). Excluding the long-time introduced M. arenaria showed that % NIS contribution in recent years was highest in Limfjorden (5%), followed by the western Baltic Sea (4%), the Belt Sea (2%) and with less than 1% in the North Sea/Skagerrak and Kattegat regions. Analysis of data including Mya arenaria, however, elevated the NIS contribution by 1% in Limfjorden and the Belt Sea and 3% in the western Baltic region.
Macroalgae. A total of 327 species was recorded in the entire monitoring data set, with 57, 299, 141, 253 and 140 species from regions 1 to 5 respectively. Ten macroalgal NIS were observed in the entire data set ( Table 2) with an additional two species recorded outside the monitoring programme (Table S1). Of the ten NIS recorded in the monitoring programme, eight were observed in the Limfjorden, which is the only region where there has been an observed increase in NIS numbers since 1989 ( Figure 7A). Based on presence/absence data, none of the regions showed significant changes in the relative contribution of macroalgal NIS to total species number ( Figure 7B).
Fish. A total of 140 species was recorded in the entire data set with 112, 104, 79 and 60 species from region 1, 2, 4 and 5, respectively. Only one recording of a NIS was made in the entire data set (

Impact of NIS on community structure over the last three decades
Phytoplankton. An ordination plot of phytoplankton based on nonparametric multi-dimension scaling (nMDS) of data from the sampling period 1989-2014 showed that observations grouped together region by region had different communities in all regions ( Figure 8A). However, data from the North Sea and from the western Baltic Sea varied more than observed for the Kattegat, the Belt Sea and in the Limfjorden. The similarity between regions was related to the North Sea -Baltic Sea salinity gradient such that the North Sea phytoplankton community was most similar to the communities in the Limfjorden and the Kattegat which, in turn, were most similar to the Belt Sea communities.
The western Baltic Sea and the North Sea had the lowest community similarity. Clustering within regions was not strongly related to the relative community contribution of NIS which generally accounted for less than 20% of the phytoplankton records. The SIMPER analysis showed that the phytoplanktonic NIS accounted for less than 13% of the within-region similarity (Table 3). When the analysis was applied to each of the periods 1989-1997, 1998-2005 and 2006-2014, a consistent pattern emerged showing that phytoplanktonic NIS contributed increasingly to the withinregion community similarity through the three consecutive periods in all five regions (Table 3). Phytoplanktonic NIS contributed most to the similarity in the North Sea samples and least in the western Baltic Sea. NIS contribution to the diversity (gamma-diversity) followed the same pattern with a gradual increase during 1989-2014 (data not shown).  Table 3. The contribution of NIS to total community similarity within regions for each of the three periods for the different taxonomic groups investigated. Fish were not included because only two NIS were recorded in the monitoring data set (see Table 2). Zooplankton. The nMDS ordination for the 25-year time series, showed a regional grouping, except for the North Sea area, where no zooplanktonic NIS were recorded ( Figure 8B). However, the regional communities did not segregate as distinctly with zooplankton as with the phytoplankton, and macroalgal (see later) communities. The regional clustering was furthermore not related to NIS contribution (% NIS). The SIMPER analysis showed that zooplanktonic NIS generally explained less than 6% of the total within community similarity and only a weak tendency to increase their importance during the observation period (Table 3).
Benthic invertebrates. In the open water data set, only three species were recorded; Caulleriella killariensis, Crepidula fornicata and Mya arenaria.
The former two were only recorded twice. The contribution of benthic invertebrate NIS to species diversity in this dataset was marginal as they only increased alpha diversity (quantified as species per 143 cm 2 sample) by 0.1% (from 10.16 to 10.17, n = 2627). NIS contribution to the overall Shannon-diversity index was even smaller (from 2.593549 to 2.594961, n = 2627). While richness is not influenced by benthic invertebrate NIS, the nMDS analysis of similarities in species composition showed a clear clustering with respect to the region ( Figure 8C). Although the clustering was not related to the relative contribution of benthic invertebrate NIS, the SIMPER analysis showed a consistent pattern with NIS becoming a more important part of the community composition over time ( Table 3). The relative importance of NIS to within-region similarity increased from ca. 2% in 1989-1997 to 5% in 2006-2014 when comparing all five regions (Table 3).
Macroalgae. The community composition of the macroalgal communities showed similar clustering as seen for the phyto-and zooplankton communities ( Figure 8A, 8B and 8D). Macroalgal NIS were reported in all five regions, but with a very low contribution in the North Sea and Skagerrak and western Baltic Sea regions ( Figure 8D). The communities from each of the regions were different from each other, however, a low number of sampling stations in the North Sea and the western Baltic Sea caused these communities to appear more scattered than communities in the Kattegat, the Belt Sea and the Limfjorden. Nevertheless, the clustering of community data followed the overall salinity pattern such that the North Sea macroalgal community was most similar to the communities in the Kattegat and Limfjorden. The SIMPER analysis showed that macroalgal NIS were only important for the community composition in Limfjorden where an increase in their importance was observed during 1989-2015 (Table 3).
Fish. We found no changes in fish community composition associated with NIS, simply because only 1 NIS-species was found in only a few samples.

Characterisation of environmental conditions favouring NIS presence
We found weak and mostly non-significant correlations across regions between environmental variables and the relative contribution of NIS to total species richness (Figure 9). There was a tendency for the phytoplanktonic and macroalgal NIS groups to increase in prevalence (% NIS) under high chlorophyll a levels, low water clarity and high-temperature conditions (see Figure 9), suggesting that macroalgal NIS represent a higher proportion of the entire species-pool under warm and eutrophic conditions. Performing the analysis separately for each region also showed low, non-significant correlations (data not shown). Also, there was a tendency for benthic invertebrate NIS (% NIS) in the North Sea/Skagerrak and Kattegat to be high when associated with low primary production, low total nitrogen and low water temperature but also high water-clarity (Figure 9). The environmental  Figure 10. Spearman rank correlations between presence/-sampling matrices for each taxonomic group and corresponding matrices of environmental factors. Only variables which contributed significantly (p < 0.05) were included. Environmental and taxonomic variables are the same as in Figure 9. Also, we included the relative contribution of NIS (%NIS) in the analysis.
conditions characterising high NIS presence were, however, opposite in the Belt Sea and western Baltic (Figure 9). Given that very little zooplankton was available at the regional scale, correlations were all non-significant, and tendencies here were non-conclusive.
Finally, species community structure correlated strongly with gradients in salinity and % NIS (Figure 10). Data on water clarity, primary production and temperature showed a weak but significant correlation with the community data, whereas total nitrogen showed no significant (p < 0.05) correlation ( Figure 10).

Discussion
A total of 56 NIS were recorded in the Danish marine monitoring programmes over a 25 year period, and an additional 21 NIS have been  observed through other data collections and observations. The increase in NIS observed in the monitoring program alone since 1989 (26 NIS see Figure 2) is much higher than the average number of NIS recorded in Baltic Sea countries, but lower than the 132 NIS reported for the entire Baltic Sea (Ojaveer et al. 2017). Within Danish waters, most NIS were found in the inner Danish waters compared to the brackish western Baltic Sea region and North Sea region. This supports recent large-scale observations, that there are large differences within the Baltic Sea in terms of the number of found introduced species (Thomsen et al. 2008b;Ojaveer et al. 2017;Tsiamis et al. 2018Tsiamis et al. , 2019).

Data quality
Knowledge and data gaps, including uncertainties of which species are NIS (e.g., Tsiamis et al. 2019;Gómez 2019), and significant differences in sampling design, effort and focuses across the various monitoring programmes, are considered problematic for comparison of NIS trends within and between countries (Palialexis et al. 2014). According to Gómez (2019) this includes species such as Odontella rhombus, Peridinella danica, P. catenata, Pseudosolenia calcar-avis, Trieres regia and maybe also Chaetoceros peruvianus as the most notable examples listed accordingly to our ranking. Although the Danish marine monitoring programmes were not established to investigate trends in NIS, the database have previously been used to investigate the impact of benthic invertebrate NIS (Thomsen et al. 2008b) and macroalgal NIS (Staehr et al. 2000;Thomsen et al. 2007).
Our study is, however, the first which uses the national monitoring data to perform a quantitative analysis of long-term changes to marine NIS covering all biotic elements of Danish marine monitoring. To enhance the robustness of our analysis, the applied NIS lists were compared with lists produced by European expert groups (EASIN and AquaNIS). As presence/absence is the most robust level of data collected in both the Danish and other monitoring programmes, it is appropriate for indicators to be developed and assessments conducted with this data. Our number of phytoplankton NIS is comparatively larger than for the other taxonomic groups. Here we agree with Tsiamis et al. (2019) that "unicellular plankton NIS should be treated with caution (e.g. flagged with high uncertainty) until further research clarifies their enigmatic status".
Due to the heterogeneous data formats we performed the quantitative analysis for the 1989 to 2014 period on highly aggregated data (grouped into five large regions) to reduce bias comparing different monitoring programs. Although the applied GLM modelling is commonly used to reduce the importance of spatial and temporal heterogeneities in the sampling efforts (e.g. Fahrig 1992), it cannot be excluded that our results are affected by changes to the sampling design (including a number of sampled stations and frequency of sampling) from 1989 to 2014. The increase in NIS numbers occurred parallel to an increase in the total number of recorded species, both of which increased with the number of samples recorded over time most notably during the first years of sampling (see Figure S1). This pattern is always expected from species accumulation curves where the total number of recorded species versus the cumulated sampling effort follow a saturation curve (e.g. Witman et al. 2004). Thereby it is impossible to disentangle time of introduction from the rareness of a NIS. However, analysing trends in the relative contribution of NIS (% NIS) is robust to changes in sampling effort because the majority of the observed new species are not identified as NIS at time of sampling. Accordingly, while sampling effort is likely to have affected the total number of NIS observed, we do not expect that effort has affected our assessment of the impact of NIS measured as relative contribution of NIS to overall species number and to changes in community structure.

Long-term trends in NIS
Out of the 77 documented NIS in Danish waters, 30 species were observed before the establishment of the marine monitoring program in the 1980s, followed by a rapid increase hereafter. This supports recent reports of a gradual low rate of introduction during the first part of the 20 th century, followed by a much higher rate in recent decades mostly related to shipping, deliberate stocking and secondary spread of NIS already established in adjacent waters (Staehr and Thomsen 2012;Ojaveer et al. 2017). Also, the long history of introductions in Danish waters (e.g. Mya arenaria was likely introduced by the Viking expansion to the North American continent (Jensen and Knudsen 2005)) suggesting that many NIS are well established, as observed in other European Seas (Tsiamis et al. 2018). Our analysis of recent temporal trends in the relative contribution of NIS revealed significant and exponential increases in phytoplankton and benthic invertebrate species in all five regions (but no trend for macroalgae, zooplankton or fish). To our knowledge, such strong increases have only been documented for a few taxonomic groups in the Baltic Sea region (e.g. Staehr et al. 2000;Thomsen et al. 2009). Our analysis shows that ranking NIS according to their relative contribution based on P/A data (% NIS) gives a good indication of how well established NIS are within a region and allows assessment of possible trends. In comparison, ranking NIS according to relative abundance (e.g. measured as number of individuals, coverage, biomass) of the total community can be used to assess the potential impact of the NIS, e.g., through competition with other species, effects on other trophic levels, or on cycling of carbon and nutrients and oxygen demand (e.g. Staehr et al. 2000;Pedersen et al. 2005;Riisgård et al. 2015). Given that true abundance data exists, such data may provide a stronger measure to assess impacts of NIS on the indigenous communities.

Impact of NIS on community structure
Our multivariate analysis showed that NIS explained up to ten percent of the annual changes in the community composition, indicating that NIS populations have established to an extent where large scale community structures of entire Danish marine communities have been impacted. While this result agrees with recent analysis for the Baltic Sea (Ojaveer and Kotta 2015) the impact, however, varies greatly among the taxonomic groups and regions. Overall, the community analysis showed that in particular the phytoplanktonic NIS, have become much more dominant and become a more characteristic community component which also suggests an increasing impact on the plankton community structure and function.
Furthermore, the multivariate analysis showed that NIS contribution to total community similarity increased over time, highlighting that NIS has become a more characteristic component of Danish marine communities in general. A common result for all taxonomic groups was furthermore that regions typically clustered in space such that the NIS-component of the communities in the North Sea data resembled that of the connecting region; the Kattegat and Limfjorden. This suggest that the connectivity of water masses (regions) play a key role for introduction and spread of NIS among and within regions where the North Sea seems to be the "window" for introduction of NIS to the broader Baltic region. Other human factors such as shipping traffic together with natural environmental conditions are furthermore generally found to be responsible for the spread and distribution of NIS at a regional sea scale (Kotta et al. 2016).

Importance of some dominating NIS
Although the community analysis indicated that NIS presence had not caused a major shift in the overall community structures, there are several examples that NIS still can dominate communities locally on shorter time scales. A few recent examples in Danish waters of highly invasive NIS are Pseudochattonella sp. (phytoplankton) (Jakobsen et al. 2012) Sargassum muticum (macroalgae) (Staehr et al. 2000;Pedersen et al. 2005;Thomsen et al. 2007), Marenzelleria viridis (invertebrate) (Kotta et al. 2001), Mnemiopsis leidyi (ctenophore, only in the Belt Sea) (Riisgård et al. 2010(Riisgård et al. , 2015, Acartia tonsa, (zooplankton) (Ojaveer and Kotta 2015) and Neogobius melanostomus (Kotta et al. 2016) (fish, only in the Belt Sea) (Azour et al. 2015). The relatively low occurrence of these NIS at the national level in the entire monitoring program can be due to several conditions: 1) their dominance patterns are cancelled out when samples are pooled across entire regions and years; 2) impact will only be seen when analysing abundance data, or 3) for some of the species we simply don't have regular monitoring data available. In the following, we provide some examples on selected dominant NIS and their documented impact in Danish waters.
Pseudochattonella sp., a phytoplankton species which produces fish toxins, has since its first observation in 1998 caused several blooms during late winter to early spring in the Kattegat region, where it can contribute most of the phytoplanktonic C-biomass (Jakobsen et al. 2012). In comparison, Karenia mikimotoi (syn: Gyrodinium aureolum), another toxic phytoplankton species, have been reduced markedly over the last decade, co-occurring with a significant reduction in nutrient load of the Danish coastal waters (Riemann et al. 2016). While some species (e.g. Pseudochattonella sp.) continues to bloom, most phytoplankton biomass as chlorophyll are decreasing as nutrient levels decrease (Riemann et al. 2016), suggesting that functionally similar species may interact very differently with changes in the environment. Also, two phytoplankton NIS (Phaeocystis pouchetii and Lepidodinium chlorophorum) which may contribute significantly to the overall biomass, occurs rarely as they are not among the ten highest ranking based on % NIS. This indicates that while some NIS become a more stable part of the communities, others are found in extreme biomasses during ephemeral blooms. Overall, these results support recent findings that closely related NIS do not necessarily behave ecologically similar but can function differently depending on their recipient ecosystem (Evangelista et al. 2019).
The macroalgae, Sargassum muticum, has since its introduction to the western part of Limfjorden in 1984, successfully established in most of Limfjorden and by far become the most dominating species in this region, significantly impacting the macroalgal community structure (Staehr et al. 2000) and altered the carbon and nutrient cycling (Pedersen et al. 2005). This species is a good example of how powerfull a continuous monitoring program is to document changes in NIS abundance and impact on local and regional scale.
Mnemiopsis leidyi has since its introduction in 2007 become an important component of the zooplankton community where it can form blooms on smaller spatial and temporal scales (Riisgård et al. 2010(Riisgård et al. , 2015. Mnemiopsis leidyi has a high filtration capacity, compared to native jellyfish and can therefore exert top down control on several larvae of native planktonic species (Riisgård et al. 2015). Furthermore their abundance is sensitive to irregular environmental changes (Schnedler-Meyer et al. 2018) and the impact of Mnemiopsis leidyi do not appear to be as dramatic initially predicted, based on its native ecology in the Black Sea (Huwer et al. 2008;Kideys 2002). The copepod Acartia tonsa is another example of a "semi"-cryptogenic NIS which was recorded as early as 1927 in European waters (Rémy 1927). Although its first record in the DNAMAP program was in 1991, it has most likely already been well-established then in parts of the Danish waters (Jespersen 1949) and the Baltic Sea (Brylinski 1981) where it is suggested to have ecosystem impacts (Ojaveer and Kotta 2015). Acartia tonsa is however by now considered a common zooplankter in the Baltic (Brylinski 1981) where it in some areas have taken an important key role as grazer on phytoplankton and prey for benthic invertebrates and fish (Tiselius et al. 2008) thereby being an important contributor to the food web efficiency and nutrient cycling.
The non-indigenous fish species Neogobius melanostomus was only observed once in the National data set, but is nevertheless reported as a very dominant NIS in Danish waters (Azour et al. 2015). The species lives in shallow low saline waters of the western Baltic Sea where it outcompetes indigenous species and it spreads northward at a rate of ca. 30 km per year (Azour et al. 2015). Its preference to shallow waters explains why the species is not caught regularly in the more offshore commercial fish surveys.

Importance of environmental conditions
Invasiveness of NIS depends largely on the environmental characteristics (such as salinity, temperature range, depth, etc.) of an ecosystem (Zaiko et al. 2007) in combination with life history characteristics of the invader. Furthermore, NIS typically invades habitats and ecosystems that are strongly modified by human activities (Thomsen et al. 2016). For example, certain ecosystem processes, such as resistance to invasions, may be impeded in stressed low-diversity systems to favour NIS. In the Danish waters, eutrophication, physical disturbance of the seafloor, hazardous substances, and climate changes may all function to lower resistance against invasions. An area such as the Limfjord is known to be strongly affected both by eutrophication and physical disturbance of the seafloor (e.g. Christiansen et al. 2006). This may partly explain why regions such as Limfjorden are so strongly impacted by NIS. In addition, high number of NIS here may also reflect the intensity of introduction pathways and connectivity to the North Sea, suggesting that the gradient in NIS with salinity reflects natural dispersal. Connectivity to donor regions could therefore be an important factor to consider in relation to our findings. Understanding the locations most at risk and focusing monitoring at these points is an important element in assessing rates of introduction of NIS as well as facilitating rapid response to invasive new arrivals (Tidbury et al. 2016).
To elucidate the causes for the increasing contribution of NIS to Danish marine communities, analyses at more local scales are needed in order to link the occurrence of NIS to environmental variables with greater certainty (Carlton 2002). Such analyses will furthermore clarify potential speciesspecific effects on ecosystem and habitat levels and should be pursued.
The use of de-trended time series allowed us to disentangle changes and gradients in environmental conditions from the longer-term co-varying year-to-year changes. This approach decreases the risk of getting a "false positive" correlation (type I error) due to the time trend caused by other factors than those investigated in this study (Kurnaz 2004). The analysis aimed to evaluate if communities influenced by NIS were characterised by a specific environmental conditions, such as temperature, water clarity or nutrient levels. We found, not surprisingly, a general salinity-driven shift in the composition of most taxonomic groups, with decreasing species richness from the North Sea towards the low saline Baltic Sea (Middelboe et al. 1997;Josefson and Hansen 2004;Telesh and Khlebovich 2010). Similarly, our correlation of environmental conditions to changes in % NIS and community structures also demonstrated this strong influence of salinity for both benthic and pelagic communities in the Baltic Sea (Middelboe et al. 1997;Josefson and Hansen 2004;Telesh and Khlebovich 2010). Nevertheless, most correlations were weak, and the conditions affecting NIS prevalence cannot easily be interpreted in terms of the selected environmental parameters. Although NIS number increases with salinity, the NIS impact on indigenous species may be higher in low salinity regions, as these are considered more vulnerable from greater availability of niches (Paavola et al. 2005). Suitability of environmental conditions, invasiveness of a region and availability of dispersal vectors will therefore influence the correlation between NIS impact and environmental parameters (Paavola et al. 2005).

Conclusions and recommendations
It is difficult to document effects of NIS on community structure and ecosystem processes because effects depend on not only presence but also abundance of a NIS, which vary in space and time, and per-capita effects are highly context dependent (e.g. Parker et al. 1999;Thomsen et al. 2016). Nevertheless, using commonly available presence/absence data, our analyses showed that the contribution of NIS to species richness and community structure on a nationwide scale has increased over time for almost all functional groups. Comparing different time periods we document that NIS represented up to twelve percent of the entire community for some groups during recent times (with much larger impacts at the local scale in some regions, especially the Limfjord). Where possible, the assessment of NIS influence on the indigenous species, provided in this study, should be supplemented with investigations of abundance based impact assessment of specific NIS in local areas. In Danish waters this has previously been done for Sargassum muticum in Limfjorden (eg. Staehr et al. 2000;Pedersen et al. 2005) and Mnemiopsis leidyi in Danish coastal waters (Riisgård et al. 2010(Riisgård et al. , 2015. Such data enable a more thorough examination NIS impact on habitat function. Only two zooplankton NIS and two fish and no parasitic NIS were recorded in the monitoring programmes. Other data have shown the presence of these types of NIS in Danish waters (Table S1) highlighting that the current sampling and analysis (e.g. microscopically exanimation vs. molecular analysis) methods are insufficient to monitor these NIS groups. Finally, because detection of new species and NIS are sensitive to sampling effort, it is important to maintain a relatively stable monitoring effort. The number of sampling stations and sampling frequencies should adequately cover the different regions in Danish waters and the seasonality of key organisms (planktonic organisms in particular).
The assessment processes we present within this study provide a robust means by which the influence/impact of NIS can be assessed at over long-term trends and at the community level. As this process uses presence/absence it can be easily adopted from existing marine biodiversity monitoring data without the need to collect abundance data. This process could therefore be effectively adopted as a means of assessing impact under the MSFD as required for criteria D2C3.