Benthic biodiversity on old platforms , young wind farms , and rocky reefs

Joop W. P. Coolen*, Babeth van der Weide, Joël Cuperus, Maxime Blomberg, Godfried W. N. M. Van Moorsel, Marco A. Faasse, Oscar G. Bos, Steven Degraer, and Han J. Lindeboom Wageningen Marine Research, PO Box 57, 1780 AB Den Helder, The Netherlands Chair Group Aquatic Ecology and Water Quality Management, Wageningen University, Droevendaalsesteeg 3a, 6708 PD Wageningen, The Netherlands Central Information Services Department, Ministry of Infrastructure and the Environment, Rijkswaterstaat, Zuiderwagenplein 2, 8224 AD Lelystad, The Netherlands Ecosub, PO Box 126, 3940 AC Doorn, The Netherlands eCOAST Marine Research, DOK41, Voltaweg 11c, 4382 NG Vlissingen, The Netherlands Aquatic and Terrestrial Ecology Marine Ecology and Management Section, Royal Belgian Institute of Natural Sciences Operational Directorate Natural Environment, Gulledelle 100, 1200 Brussels, Belgium


Introduction
The introduction of artificial hard substrates in areas dominated by sandy seabeds increases habitat available to epifouling organisms (Zintzen and Massin, 2010;Lengkeek et al., 2013;De Mesel et al., 2015).As a result of the increasing number of offshore energy devices in the North Sea, the amount of artificial hard substrate increases steadily (Krone et al., 2013).The European Wind Energy Association estimates that the North Sea will hold up to 62 GW of offshore wind energy installations in 2030 (EWEA, 2015).Assuming 5 MW turbines, this is the equivalent of 12 400 turbine foundations, all providing artificial hard substrates to epifouling organisms.
Oil and gas companies have been installing platforms in the North Sea since the 1960s (Shepherd, 2015) and expanded into the Dutch continental shelf (DCS) from the early 1970s onward (EBN, 2014).The oldest North Sea offshore wind farm was constructed in 2002 (Leonhard and Christensen, 2006).Compared with oil and gas, offshore wind is still a young industry.Most oil and gas platforms in the North Sea were constructed using steel jacket foundations, which structurally resemble offshore windturbine foundations, offering similar substrate to epifouling organisms.It is therefore expected that offshore wind farms and oil and gas platforms house similar benthic communities when other environmental variables are constant.However, other variables often differ between locations.For example, epifouling communities on offshore installations evolve over time with dominance changing among species (Whomersley and Picken, 2003).Furthermore, species richness may increase with installation age (van der Stap et al., 2016).
In recent years the number of offshore wind farms in the North Sea increased strongly.At the same time, a large-scale removal of offshore oil and gas installations is expected over the coming years (Ahiaga-dagbui et al., 2017).Many oil and gas platforms are reaching their financial or technical end-of-life and will be decommissioned (Royal Academy of Engineering, 2013).To prevent high removal costs, loss of fishing habitat and to maintain the artificial reef effect of structures, alternative uses for the foundations of these platforms are proposed (Day and Gusmitta, 2016).These include so called rigs-to-reefs (RTR) schemes, in which foundations are to be left at sea as artificial reefs, either left-in-place, toppled over, or relocated (Picken and McIntyre, 1989;Macreadie et al., 2012).Similar proposals (named renewables-to-reefs; RTR) may be expected once offshore wind farms reach their end-of-life (Fowler et al., 2015;Smyth et al., 2015).To evaluate the ecological potential of such alternatives, more insight is needed into the current state of the epifouling communities present on the structures.Offshore reefs of natural origin in the Netherlands are primarily composed of rocky substrate (Coolen et al., 2016).These natural reefs, including their epibenthic community, qualify for protection under the European Habitat Directive (European Commission, 1992).If artificial habitats hold similar communities, they could be considered as desirable as well.This insight could then be used in the evaluation of future RTR proposals or in impact assessments for new installations.To evaluate the added value, if any, of these structures to the North Sea ecosystem structure and functioning, the following questions should be answered: Which organisms have colonised these artificial structures?What characteristics of these structures influence their community composition?What are the habitat preferences of the species currently found on these structures?Are the communities found similar to those on natural reefs?What habitats are particularly suitable for non-native species?
Many variables may be important drivers of species' presence on reefs (Coolen et al., 2016;Herbert et al., 2017).Depth, for example, is widely known to influence the vertical distribution of species on artificial reefs: some species are present only on deeper reefs (Coolen et al., 2016) and species richness varies with depth and age (van der Stap et al., 2016).Structural differences between reefs also explain part of the variation in the associated community (Coolen et al., 2015).Straight steel surfaces differ from rocky reefs with a more complex surface area including holes and smallscale variation in surface orientation.Therefore, it is hypothesized that species richness of benthic communities is higher on rock than on smooth steel surfaces.
Relations among species also influence species composition.Species with disproportionally large effects on their environment by modifying the habitat and associated community are defined as keystone species (Paine, 1969).Mills et al. (1993) described several types of keystone species, including keystone predators and modifiers.Keystone predators forage on species that otherwise would become dominant, reducing their abundance, which may increase local species richness.Alternatively, keystone predators may consume prey in such numbers that species richness is negatively affected while the predator becomes dominant.Keystone modifiers affect the survival of many other species.Their presence therefore could impact species richness.Several species present on artificial reefs may fit these types defined by Mills et al. (1993).The sea urchin Psammechinus miliaris (Mu ¨ller, 1771) could fulfil a role as keystone predator.In low numbers, such predators may increase local species richness by creating patchiness in the fauna cover, allowing new species to colonise the substrate (Menge and Sutherland, 1976).This behaviour has been observed on artificial structures in the North Sea as well (J.W.P. Coolen, G. van Moorsel, pers. obs.).The sea anemone Metridium dianthus (Ellis, 1768) is known to smother and thereby kill other species (Kaplan, 1984;Nelson and Craig, 2011).Coolen et al. (2015) noticed that samples dominated by M. dianthus held less species than samples without this sea anemone.Therefore, it is hypothesised that M. dianthus has a negative influence on species richness through space monopolization and predation of larvae, functioning as a keystone predator.However, M. dianthus could also be seen as a keystone modifier since it reduces substrate available to other species.Mussels Mytilus edulis increase habitat heterogeneity resulting in increased species richness (Drent and Dekker, 2013).This bivalve could therefore be considered a keystone habitat modifier.A similar effect has been shown from the hydroid Tubularia indivisa Linnaeus, 1758, which creates complex growth forms facilitating a large number of associated species (Zintzen et al., 2007).
The effects of some of these variables have been evaluated to some extent in previous research (Zintzen et al., 2007;Walker and Schlacher, 2014;Coolen et al., 2015Coolen et al., , 2016;;van der Stap et al., 2016).However, these assessments only included a small number of variables, e.g.distance to coast, age, and depth on platforms (van der Stap et al., 2016) or orientation on concrete artificial reefs in southern Portugal (Moura et al., 2008).Therefore, the work reported here aims to assess the combined effect of a larger set of variables.Using samples taken on old offshore oil and gas platforms (age 15-40 years), a young wind farm (5 years) and a natural reef on the DCS, the effect of depth, age, disturbance by marine growth removal, season, substrate type, and the presence/ absence of potential keystone species on species richness and composition was investigated.

Study area
Samples were taken on natural and artificial reefs on the DCS in the North Sea.Most of the DCS seabed consists of fine mud and fine to coarse sand [European Biodiversity data centre (BDC), 2016], interrupted by isolated gravel beds including cobbles and boulders on the Cleaver Bank, Texel Rough, Dogger Bank, and Borkum Reef Grounds (BRG) (Duineveld et al. 1991;Van Moorsel 2003;EMODnet 2015;Coolen et al. 2016; Figure 1).

Platform data
Five oil and gas production platforms, operated by ENGIE Exploration & Production Nederland B.V. (ENGIE) were sampled between April 2014 and October 2015 (Figure 1, Table 1).The platforms were deployed in the period 1972-1999 and are supported by steel jacket foundations with cathodic protection with anodes.The seabed surrounding these installations is between 22 and 32 m deep and has been covered with cobbles and boulders (0.05-0.5 m), to prevent scouring of the sediment around the legs (rock dump, also known as scour protection).The area covered by rock dump varies between 500 and 8800 m 2 per platform, depending on local currents and pipelines present (G.Menijn, ENGIE, pers.comm.).

Sampling methods
Platform samples were taken during regular inspection, repair, and maintenance (IRM) activities of the operator by a diver using a surface-supplied airlift sampler (Coolen et al., 2015) constructed of 50 and 75 mm standard PVC pipes.A 0.05 m 2 metal frame was attached to the structure using magnets and all fauna within the frame was scraped off using a putty knife.Removed fauna was then sucked up by the airlift sampler and collected in a net with a 0.5 mm mesh size.In between sampling, nets were replaced and sealed with a PVC screw lid.Three replicate samples were taken at 5 m depth intervals.At two locations, cobbles present on the seabed around the installation were collected by hand from a 0.05 m 2 metal frame and placed in zip-lock bags.Samples from this rock dump were taken randomly, but distance among them was at least 5 m.
Samples were pre-processed on board: They were removed from the nets using running seawater, placed in 2.5 l containers and covered with fresh seawater.An oversaturated mentholseawater solution was added and the samples were stored at 4 C for a minimum of 5 h.Next borax-buffered formaldehyde was added with fresh seawater to reach a final concentration of 6% formaldehyde.Within 7 days after sampling, all samples were rinsed to minimize formaldehyde residue and placed on ethanol 70% with 3% glycerol for long-term storage.
All samples were sorted into taxonomic groups.Taxa were then identified to the lowest taxonomic level possible, preferably  species.The World Register of Marine Species (WoRMS Editorial Board, 2017) was used as a reference of taxonomic nomenclature.

Additional data
For comparison of platform data with natural reefs, unpublished raw data were acquired from a study of the benthic diversity of the rocky reefs at the BRG carried out in August 2013.The BRG were sampled using identical methods as used in the platform study although the airlift used was fed by air from scuba equipment and divers followed a 50 m transect line.Samples were taken from large boulders present on the sandy seabed and in gravel fields by scraping off the epibenthic community (Coolen et al., 2015).
For a comparison with younger installations, raw data from the study of wind turbine foundations and rock dump at the Princess Amalia Wind Farm (PAWF) were used (Vanagt et al. 2013;Vanagt and Faasse 2014).The PAWF turbine foundations were deployed between November 2006 and March 2007 and surveys were carried out in October 2011 and July 2013.Two samples of 0.056 m 2 were taken from the intertidal zone and at 2, 5, 10, and 17 m using a putty knife and collected in 0.25 mm mesh size nets.Four wind turbine foundations were sampled per survey.Small rocks were collected from the scour protection at random locations around the foundations.
Taxonomy of the BRG and PAWF data sets was synchronized with the platform data using the WoRMS AphiaID service (World Register of Marine Species, 2015) in R version 3.3.0(R Core Team, 2016) using the worms package (Holstein, 2017) and combined with the platform data in a single dataset.For each sample, location name, depth (m), substrate type (steel/rock), orientation (horizontal/vertical/mixed for smaller rocks), sampling month, and age (years; only artificial structures) were registered.Age was considered to be the number of years since construction, or when cleaning of the construction had taken place: the number of years since cleaning of the substrate.Offshore operators regularly remove part of the epifouling down to À10 m depth to reduce wave load on the structure, depending on thickness of encrusting organisms.Since we expected cleaning might have a different effect on the community than age of the structure, a factor variable (yes/no) indicating whether the sampled depth had been cleaned in the last 5 years was also included in the dataset.

Data analysis
Data for all species were combined into a single dataset, consisting of 191 samples.In total, 161 samples originated from steel and 30 from rocky substrate.Within this set, species observations in 35 samples (all from steel) were only available as presenceabsence data.These were excluded during modelling but included in the list of taxa observed.Algae and copepods were removed from the dataset since these were not registered in some studies.Observations of juvenile or damaged specimens identified to a level above species were either combined with a species in the same taxon or removed when more than one species was present in the taxon.Only juvenile and damaged taxa in samples with no species in the same taxon were left at the higher level, following Coolen et al. (2015).Given the differences in sampled area (PAWF: 0.056 m 2 , platforms, and BRG: 0.05 m 2 ), abundance was converted to individuals per m 2 .Colony forming species were not quantified in all samples and in the different projects various quantification methods had been applied.Therefore all Hydrozoa, Bryozoa, Porifera, Ascidiacea, Entoprocta, and Alcyonium digitatum were converted to presence/absence.For data analysis, R version 3.3.0(R Core Team, 2016) and RStudio version 0.99.486(RStudio Team, 2015) were used.Both univariate and multivariate analyses were performed.

Univariate analysis
For the univariate analysis, species richness (S; number of unique species) was calculated for each sample and used as response variable.The protocol by Zuur et al. (2010) was used for data exploration.The presence of outliers, multicollinearity, and relations between species richness, other biotic, and abiotic variables was assessed using boxplots, Cleveland dotplots (Cleveland, 1985), pairplots, Pearson correlation coefficients, variance inflation factors, and multipanel scatterplots from the lattice package (Sarkar, 2008).A power analysis was performed to estimate whether the sample size was sufficient to test the combined set of variables, using the pwr.f2.test function from the pwr package (Champely et al., 2017).A power of 0.8, p ¼ 0.05, a medium effect size of 0.15 (Cohen, 1988) and 16 coefficients in the full model, including multiple coefficients for the different locations and 2 for each non-linear variable were used.
Species richness may have a non-linear relation with environmental variables, such as depth (van der Stap et al., 2016).To allow for easy calculation of such relations, all models were constructed using generalized additive models (GAMs).As species richness only results in positive discrete numbers, and initial models applying a Poisson distribution resulted in overdispersed residuals, a negative binomial distribution with log link was used.Models were created with the gam function from the mgcv package (Wood, 2011).Random location and orientation effects were included in the GAM to remove dependency among samples, viz.multiple samples taken at a single platform and to correct for possible confounding effects of orientation.Orientation was not included as prediction variable as differences between orientations were mostly observed between steel and rock.Month in year (January-December ¼ 1-12) was included as a variable to include seasonal effects.Potential keystone species were selected based on personal field observations of dominant species by the authors, and published observations indicating keystoneness of species.In addition to Metridium dianthus, Mytilus edulis, Psammechinus miliaris, and Tubulariidae (Menge and Sutherland, 1976;Zintzen et al., 2007;Drent and Dekker, 2013;Coolen et al., 2015), the starfish Asterias rubens was included as a potential keystone species (Dolmer, 1998).Depth, month, and all species excluding Tubulariidae (which was registered as presence/absence) were included as a smoothed (non-linear) variable.The maximum number of knots of the smoothers was set to three to reduce overfitting of the smoother and to allow for easier visual interpretation.
Starting from a model including all variables, a set of 11 alternative models was generated.These alternatives included: Models were then compared using Akaike Information Criterion (Akaike 1973).The number of samples per variable was <40 for most models considered and therefore AICc was used (Burnham and Anderson, 2004).The model including all variables was validated to assess if underlying assumptions of homogeneity of variance and normality of the residuals were met.Model residuals were plotted against all variables in and outside the model as well as fitted values to assess model fit.
The model including all variables took the following form: Where S is the species richness for sample i. Term f() marks a smoothing function and substrate is rock or steel.The residuals e i were assumed to be normally distributed with a mean of 0 and variance of r.
An additional model was created to assess the effect of age and cleaning of the structures, by using only data for which age was known: samples taken from steel at the platforms and the wind farm.This model took the form shown above, with the addition of age (years) and cleaned (yes/no) as variables.
The total predicted number of species per structure type (platforms, PAWF, BRG) was estimated with the specpool function (vegan package, Oksanen et al. 2008) using the Chao estimate (Chao, 1987).
As De Mesel et al. (2015) reported high fractions of nonindigenous species in intertidal zones of offshore wind turbines, the difference in fraction of non-indigenous species between the intertidal (depth <2 m) and deeper zone as well as between substrate and between reef types was investigated, using a generalized linear model (GLM) with binomial distribution.Species were assumed to be non-indigenous based on the Dutch list nonindigenous species (Bos et al., 2017), with corrections based on publications showing species on the list were indigenous or cryptogenic (Korringa, 1954;Ates, 2006;Gittenberger et al., 2010;Fofonoff et al., 2017).

Multivariate analysis
For the multivariate analysis, a dissimilarity matrix (Bray-Curtis dissimilarity index; Bray and Curtis, 1957) was created by nonmetric multidimensional scaling (MDS) using the function metaMDS from the vegan package (Oksanen et al., 2008) with a minimum of 100 tries.The number of dimensions needed for the MDS was assessed using Scree plots (Cattell, 1966).The goodness of fit of the MDS was evaluated using a Shepard plot (Shepard et al., 1972).In plots of the dissimilarity matrix, clustering of samples according to location group and substrate type were visualized using the ordiellipse function (vegan package).
A PERMANOVA (10 000 permutations, Bray-Curtis dissimilarity index; Anderson 2001) using the untransformed abundance data was applied to assess the effect of the same variables as used for the univariate analysis, including location as additional factor variable, using the adonis function from the vegan package (Oksanen et al., 2008).Using a MDS plot, the multivariate spread among groups (structure type) was assessed, since PERMANOVA is sensitive to differences in multivariate spread (Anderson, 2005).

Results
The complete dataset is available as Supplementary material and included 193 species and 105 taxa at a higher level, from 7 locations (Tables 2 and 3).In total, 143 species (n ¼ 88 samples) were observed on platforms, 95 species (n ¼ 92) on the foundations located in the PAWF, and 49 species (n ¼ 11) on the BRG.Total species richness on individual platforms varied between 41 on D15-A (n ¼ 6) and 86 on K9-A (n ¼ 22).Species richness per sample varied strongly between reef types (Figure 3a).
The vertical profile of frequent species on the platforms was similar to other installations in the North Sea (Krone et al., 2013;De Mesel et al., 2015).From the intertidal zone down to 10-15 m, Mytilus edulis was most observed, in densities up to 23 000 individuals per m 2 , changing into a Tubulariidae-Jassa herdmani community between 10 and 20 m.Jassa herdmani were found in densities of up to 1.45 * 10 6 individuals per m 2 .Below 20 m, the community on the steel structures was most frequently inhabited by Metridium dianthus, with up to 720 individuals per m 2 .
The following species were present in at least 50% of the samples on the steel of the platforms: the amphipods Jassa herdmani (Walker, 1893) and Stenothoe monoculoides (Montagu, 1815), mussel M. edulis, sea anemone M. dianthus, bryozoan Electra pilosa (Linnaeus, 1767), and anomuran crustacean Pisidia longicornis (Linnaeus, 1767).Of these common steel species, M. dianthus and E. pilosa were also present in >50% of the samples from BRG rocks and rock dump around artificial structures, whereas J. herdmani and S. monoculoides were common on rock dump, but never observed in BRG rock samples.M. edulis was found in 27% of the rock dump samples, with only a single observation on the BRG rocks.P. longicornis was found in 41% of the rock dump samples but never on rocks of the BRG.
Table 2. Full species list with presence (X), absence (À) of species per reef type, non-indigenous status of each species (X ¼ non-indigenous), minimum (min) and maximum (max) depth of observation in meters.

Continued
In total, 12 non-indigenous species (Bos et al., 2017) were present in the samples, excluding Diadumene cincta, Jassa marmorata, Monocorophium acherusicum, Balanus balanus, and Sabellaria spinulosa, which are registered as non-indigenous in the Netherlands but were deemed cryptogenic or indigenous based on literature review (Korringa, 1954;Ates, 2006;Southward, 2008;Gittenberger et al., 2010;Fofonoff et al., 2017).In total, 61 samples held no nonindigenous species while 88 held a single one and 42 samples held !2.Within the wind farm, 64% of the samples held one or more non-indigenous species and 11 non-indigenous species were found in total (9 on steel and 4 on rock dump).At platforms, 74% of the samples contained non-indigenous species, with a total of 6 observed (5 on steel, 2 on rock dump), while 2 non-indigenous species were found on the BRG, in 55% of the samples.The percentage of non-indigenous species per sample differed significantly (p < 0.001) between the intertidal (13 6 2% standard error) and the deeper parts (5 6 0.4% SE) but not between reef types (PAWF, platforms, BRG, p ¼ 0.9) or substrates (steel, rock; p ¼ 0.3).The most frequent non-indigenous species was the tunicate Diplosoma listerianum, which was most frequent between 5 and 25 m (present in 85 samples).Several species known for the broader region but not reported before in Dutch waters were observed, including the bryozoan Cribrilina punctata (Beukhof et al., 2016) and annelids Harmothoe aspera (Spierings et al., 2017) and Syllis vittata (Dias et al., 2017).
The total predicted species richness (only species, no higher taxa) for all platforms (S ¼ 172) was higher than for the PAWF (S ¼ 101) and more than double that of the BRG (S ¼ 74).The predicted richness on individual platforms ranged between 53 on D15-A and 107 on K9-A (Table 3).The average species richness per location was 70 (619 standard deviation).Beta diversity between locations was 2.75 [Sørensen index 0.29, N ¼ 7 (Sørensen, 1948)].The difference between observed and predicted species richness indicates that with 35% difference, the BRG would need much more samples to find most species present.In contrast, the difference between species found and expected is low for the PAWF (6%), indicating that the sampling performed here was thorough and not many species remained to be found.Found vs. predicted number of species on the platforms was intermediate, with a potential increase in species number of 20% with higher sampling effort.

Univariate analysis
The power analysis showed that a minimum of 142 samples would be needed, indicating that the sample size of 156 was sufficient.Model comparison showed that model M8, including all variables but A. rubens had the lowest AICc as well as the highest deviance explained (63.7%).DAICc of none of the other models was <2 compared with M8, while M0 (the full model), M3 (without substrate), M2 (without sample month), and M11 (without month and substrate) had DAICc between 2 and 4. All further alternatives had DAICc >4 (Table 4).Based on these results model M8 was selected for further evaluation.
M8 held 7 variables, and the deviance explained by the model was 63.7%, with an adjusted R 2 of 0.57 and 20.45 degrees of freedom.M8 took the following form, with a base model using factor Tubulariidae as absent and substrate as rocks: See Table 5 for the summary of this model.Model validation showed homogeneity and normal distribution of the residuals, Table 3. Observed and predicted species richness for all locations and the platforms combined with total species richness (S total ), predicted species richness (Chao) with standard error (Chao SE ), the increase in richness from observed to predicted (Increase), and number of samples (n).Table 4. Overview of full model and all alternative models considered (excluded variables indicated with -), p-value for all variables (p < 0.001:***, p < 0.01:**, p < 0.05:*, p < 0.1:., p > 0.1: ns) and estimate for linear variables (estimates for smoother variables can be derived from Figure 4).with no strong undesirable patterns.The typical non-linear relations between smoother variables and species richness are visualized in Figure 4.In general, the model shows the highest species richness in intermediate depths, with high abundance of P. miliaris, M. edulis, and presence of Tubulariidae, as well as low abundance of M. dianthus.
When assessing the impact of age and cleaning activities on species richness on steel structures, age and cleaning only explained <0.1% and 0.3% of the variation.Also, both had a non-significant effect (p > 0.1).Since both variables were nonsignificant and further results were similar to the model presented above, this model is not presented here.

Multivariate analysis
Based on the Scree plot, two dimensions were used for the nonmetric MDS, since there was no strong difference between two and three dimensions and two dimensions are easier to interpret.The MDS had a non-metric fit (R 2 ) of 0.96 with stress 0.21.The visualization showed a large overlap in steel and rock samples, and between the different locations.Part of the PAWF samples clustered strongly but this was within the more various spread of the platform samples.Overlap was present between PAWF and BRG and there was no separation between clusters (Figure 5).The ordination of the samples showed a strong relation with depth (p < 0.001).Samples taken from deep steel showed a higher similarity to BRG samples than shallow steel samples.
All terms included in the PERMANOVA had a significant (p < 0.01) effect on community composition with the exception of A. rubens (p ¼ 0.63) and M. dianthus (p ¼ 0.051; Table 6).Tubulariidae and location explained most of the variation in the dataset (R 2 ¼ 0.11 and 0.10, respectively), while depth had an R 2 of 0.07.R 2 of the other variables varied between 0.004 and  0.04.The total R 2 of the formula was 0.40, resulting in a residual R 2 of 0.60.The PERMANOVA on steel-only data showed a small significant (p < 0.01) effect of age (R 2 ¼ 0.02) but not from cleaning (p ¼ 0.48, R 2 ¼ 0.005) on the community.

Depth
The depth distribution of species richness clearly follows a nonlinear pattern with maximum species richness at intermediate depths between 15 and 25 m (Figure 4).Similar patterns have been observed on other reefs, both temperate and tropical (Clarke and Lidgard, 2000;Cornell and Karlson, 2000;van der Stap et al., 2016).Depth had a high explanatory value (7% of the total variation) in the multivariate analysis.Although not specifically investigated here, the maximum richness at intermediate depth may be explained by the intermediate disturbance theory (IDT;Connell 1978).At intermediate disturbance, diversity is at its peak.Halfway the water column, disruption from wave action is lower than near the surface but higher than at the bottom.IDT is known from wave disturbance on shallow marine sites (England et al., 2008) and was also suggested by van der Stap et al. (2016) as the main explanation for higher species richness at intermediate depths on offshore platforms.

Age and location
The univariate models did not show a significant effect of age on species richness.Although the platforms were all older than the PAWF foundations and total observed and predicted richness on the combined platforms were higher, predicted richness in PAWF fell within the range of predicted richness of the individual platform locations.These results are in line with findings by Wendt et al. (1989), who found no significant differences in species richness among shipwrecks of different ages.Their subjects, however, differed <6.5 years in age.In contrast, our multivariate analyses showed a small but significant effect of age on species composition.Whomersley and Picken (2003) studied marine growth patterns on four offshore production platforms over a period of 11 years and observed large variations in percentage cover in the dominant species.This large variation may result in nonsignificant results when trying to detect trends.Epibenthic communities on wind farms in the southern North Sea had a yearly richness cycle, with lowest richness in spring and highest in summer (De Mesel et al., 2015).This is in contrast with our results, showing a slight decrease from April followed by an increase in richness from July to October, although data outside this range were missing from our analysis.Short-term variation in richness may be much more important than long-term variation.Beta diversity indicated that there is differentiation between the studied locations.Therefore, it is likely that the higher total richness on combined platforms was caused by the scattered geographical distribution of the platforms compared with the limited distribution in the PAWF.The MDS plot (Figure 5) confirms this, as a large part of the PAWF samples cluster while the platform samples show much larger variation on average.

Key species
Some species are common on both natural and artificial reefs.Metridium dianthus and Electra pilosa are common on steel, rock dump and BRG rocks.The cosmopolitan E. pilosa (Nikulina et al., 2007) is one of the few species that was present on both substrate types.E. pilosa is a common epibenthic species in the southern North Sea (Jennings et al., 1999).It is often observed on artificial structures such as shipwrecks (Schrieken et al., 2013;Faasse et al., 2016) and oil platforms (Forteath et al., 1984;Beukhof et al., 2016).Species that are very common on steel but rare on the natural rocks are Jassa herdmani, Mytilus edulis, and Stenothoe monoculoides.These are all species that mostly occur at shallow or intermediate depths and are uncommon at the seabed around the artificial structures.Given that they are rare in rock samples, it is likely that these species are only present offshore due to artificial structures providing shallow habitat.Their presence at the rock dump around platforms is probably caused by M. edulis clumps falling off the structure (Krone et al., 2013).
Tubulariid Hydrozoa were shown by Zintzen et al. (2007) to positively correlate with species richness on a shipwreck.Results reported here confirm this.Tubulariidae create complex, bushlike environments, giving shelter to other organisms.This is the most likely cause of the higher amount of species found in samples with these Hydrozoa.The sea anemone Metridium dianthus had been suggested to have a strong negative effect on species richness (Coolen et al., 2015;van der Stap et al., 2016), this was indicated by the current results, where high abundance of M. dianthus correlates with low species richness.The species has a small effect on the community composition, shown by the PERMANOVA results, explaining 1% of the variation in the data.
The sea urchin Psammechinus miliaris had a positive correlation with species richness.This species consumes a wide variety of food sources (Hughes, 2006;Kelly et al., 2013).P. miliaris is known to clear patches of steel-substrate epifouling (Van Moorsel et al., 1989).Scraping of organisms by P. miliaris may increase heterogeneity in the local epifouling community by creating new settling opportunities for additional species.
Mussel Mytilus edulis showed a strong positive correlation with species richness in the univariate models.The structure of mussels, increasing local variation in surface orientation, is similar to that of gravelly reefs.M. edulis also provides secondary hard substrates, which augments biodiversity (Gutierrez et al., 2003;Norling and Kautsky, 2008).Likewise, presence of M. edulis increases habitat heterogeneity in the Wadden Sea and thereby positively influences species richness and composition (Drent and Dekker, 2013).

Rock vs. steel
Substrate type had only a small and insignificant effect on species richness but a significant effect on species composition.The effect varied between positive and negative in alternative models (Table 4).Although rocky substrate is much more complex than straight steel surfaces and its complexity can be positively related to species richness (Johnson et al., 2003;Kostylev et al., 2005), this was not apparent here.

Artificial vs. natural reefs
The MDS shows a gradient from shallow steel to deeper laying natural reefs.No clear separation of substrate types is present.There is a clustering of rocks of artificial and natural origin with the most overlap with deeper steel from platforms.However, no strong differentiation between the natural and artificial substrates can be concluded from the MDS.This is in contrast to earlier studies which showed communities on artificial and natural reefs differed significantly (Page et al., 2007;Wilhelmsson and Malm, 2008).
Non-indigenous species percentage was higher in the intertidal zone, which is in line with observations by De Mesel et al. (2015).Furthermore, native species that are uncommon offshore without intertidal structures were found on and around locations with intertidal zones.In case these species are unwanted, it is advisable for future RTR projects to remove the intertidal zone from abandoned installations in order to reduce the presence of intertidal species at offshore locations.

Study limitations and recommendations for future studies
The study presented here has limitations that should be kept in mind when interpreting the results.The data were collected during studies that were carried out with different objectives.The locations studied were geographically separated.This should be improved in future research by including natural reefs in the proximity of the platforms, e.g. the reefs at the Texel Rough which are located between the L15-A and L10 platforms.Some installations in the North Sea have concrete foundations.During the platform study we had the opportunity to observe footage of a concrete and a steel installation in close proximity.The dominant epibenthic species appeared to differ strongly.
We recommend concrete foundations to be included in future studies.

Optimizing artificial reefs
When reefed or newly built structures are to be optimized to provide habitat to epibenthic species that would normally inhabit natural rocky reefs, we suggest to add as much rock dump or scour protection in various sizes around the structure as possible, to increase local habitat complexity.We advise to prevent intertidal zones where possible, e.g. by removing the intertidal zone from reefed structures, cutting them well below the water surface.

Conclusion
In the Dutch part of the North Sea, there is a gradient in community differentiation from deep rocks to shallow steel.Variation in species richness correlates with depth and sampling month.Species richness is also positively correlated with abundance of Mytilus edulis, Psammechinus miliaris, and the presence of Tubulariidae.In contrast, Metridium dianthus abundance negatively correlates with richness.Adding additional rocky substrate around steel artificial structures is beneficial for epibenthic species that inhabit natural reefs.
(a) The full model with the removal of each variable, resulting in eight variations (b) A model with only the abiotic variables depth, month, substrate (c) A model with only the biotic variables M. dianthus, M. edulis, P. miliaris, A. rubens, and Tubulariidae (d) A model with only depth and the biotic variables M. dianthus, M. edulis, P. miliaris, A. rubens, and Tubulariidae

Figure 4 .
Figure 4. Modelled relation between non-linear smoothed variables and species richness (only species, no higher taxa) per sample (solid line) with standard error (dashed lines).Variables other than presented in the respective graphs were set to the following constants for prediction purposes: Depth: 15 m, month: 5 (May), substrate: steel, organisms: 0.

Figure 5 .
Figure 5. MDS plots of all reef data (stress ¼ 0.21), showing variables per sample.Black lines: isodepth lines.(a) By reef type (dark and circles: BRG, Borkum Reef Grounds [n ¼ 11], yellow and triangles: Princess Amalia Wind Farm [n ¼ 92], blue and plus sign: platforms [n ¼ 88]).(b) By substrate type (blue and diamonds: steel from platform and wind turbine foundations [n ¼ 161], black and X: rocks from BRG and rock dump around platform and wind turbine foundations [n ¼ 30]).Pair of points on the far right originate from the intertidal zone in PAWF and hold only 2 and 3 species.

Table 1 .
Locations of samples with maximum and minimum sampling depth, maximum age, location coordinates (decimal degrees, WGS1984), number of samples per substrate type and months in which sampling took place.

Table 5 .
Detailed model results model M8.

Table 6 .
PERMANOVA all data with species and higher taxa, based on Bray Curtis distance measure.