Coexistence of Pacific oyster Crassostrea gigas ( Thunberg , 1793 ) and blue mussels Mytilus edulis ( Linnaeus , 1758 ) on a sheltered intertidal bivalve bed ?

The invasive Pacific oyster, Crassostrea gigas Thunberg, 1793 was introduced in Denmark for aquaculture in the 1970s. Presently, feral populations are found in many parts of the country, with the largest populations established on existing beds of blue mussel, Mytilus edulis Linnaeus, 1758. This study was conducted in the Limfjord estuary, at Agger Tange, where C. gigas was introduced in 1972. The study site is a large cluster of raised intertidal bivalve beds inhabited by C. gigas and M. edulis in a sheltered part of the estuary. The two bivalves have some of the same living requirements, and as C. gigas have been present in the ecosystem for more than 40 years, we hypothesize that the presence of C. gigas has altered the spatial and temporal distribution of M. edulis by inducing a niche separation. The spatiotemporal development of the bivalve bed was determined using orthophotos. C. gigas and M. edulis were collected from the bivalve bed, shell lengths were converted into biomass, which were interpolated to create biomass contours and combined with modelled topography of the bivalve bed to study niche separation. The bivalve bed slowly extended northwards over a period of 11 years, where it also became more fragmented. The northern part of the bed was composed of mussel mats on top of soft sediment. This area was dominated by M. edulis, while areas in the south were dominated by C. gigas. In the southern part, the bivalve bed was composed of thick and compact sediment suggesting it represent the oldest part of the bivalve bed. There were no differences in the conditions of C. gigas and M. edulis from old or newly established areas, and there were no difference in the vertical distributions of the bivalve species. Thus, spatial and temporal separation of the two species is not pronounced at present, and thus unable to explain why they


Introduction
When introduced into new ecosystems, nonindigenous species can become invasive if they are able to establish populations, affect the structure and functioning of the recipient community, and subsequently disperse (Reise et al. 2006). The Pacific oyster Crassostrea gigas (Thunberg, 1793) originated from Japan and Southern China, but has been introduced all over the globe, where it has established populations with varying degrees of success (Ruesink 2005). It display many of the traits that characterise successful invasive species, as it exhibit a significant physiological plasticity, early sexual maturation, high reproductive output, high growth rates, and a high dispersal capacity (Troost 2010).
Crassostrea gigas is often found in the same intertidal areas as the native blue mussel Mytilus edulis Linnaeus, 1758 (Nehls et al. 2006), and tends to settle on existing mussel beds and conspecifics . The larvae have gregarious settling behaviour, and preferentially settle on shells of conspecifics (Diederich 2005), which is cued by the presence of adults (Bonar et al. 1990;Tamburri et al. 2007).
Both C. gigas and M. edulis are ecosystem engineers (Jones et al. 1994;Reise 2002; Gutierrez et al. 2003) that can create reef structures protruding from the seabed. When altering the habitat, ecosystem-engineering species often create a positive feedback, which increases their own fitness (Jones et al. 1997;Hui et al. 2004;Wright and Jones 2006). For example, settlement of C. gigas increases bottom roughness and hence water turbulence in the benthic boundary layer (Reise 2002;Markert et al. 2010;Styles 2015), which enhances food availability (Lenihan 1999). Furthermore, oysters reefs and mussels beds tend to change predatorprey interactions (Bartholomew et al. 2000;Grabowski 2004), as predation is reduced by structural complexity of the habitat (Dolmer 1998;Hill and Weissburg 2013;Waser et al. 2015). Oyster reefs therefore function as a refuge for a wide range of species (Eschweiler and Christensen 2011;Kingsley-Smith et al. 2012;Norling et al. 2015).
Crassostrea gigas has a competitive advantage over M. edulis as it has a higher absolute growth rate, reaches a larger size, and is less susceptible to predation (Dare et al. 1983;Troost 2010). In terms of competition for food C. gigas and M. edulis have comparable weight-specific clearance rates (e.g. Troost et al. 2009), and can selectively filtrate for different algae species (Bougrier et al. 1997). However, being a larger animal, and therefore having a higher average individual clearance capacity, C. gigas might be able to reduce the food availability for M. edulis. Furthermore, being large also allow C. gigas to protrude further into the water column thereby potentially increasing its food availability (Diederich 2006).
Competitive superiority may cause C. gigas to spatially displace M. edulis to less favourable locations, as is observed at different locations in the Wadden Sea (e.g. Diederich 2005; Eschweiler and Christensen 2011). Our study site is an intertidal bivalve bed located in the western part of the Limfjord, Denmark, which is inhabited primarily by C. gigas and M. edulis. C. gigas was introduced in 1972 for aquaculture purposes (Jensen and Knudsen 2005). At present, feral populations exist at several sites in the Limfjord, and these populations consist of several size cohorts, which indicate recurring recruitment (Wrange et al. 2010;Groslier et al. 2014;Holm et al. 2015). In a previous paper, we show that C. gigas, despite being in the ecosystem for more than three decades, has only had moderate establishment success (Holm et al. 2015). When resources are limited, species that occupy the same niche cannot occur in the same ecosystem (Hutchinson 1957), as one species will exclude the other (Harding 1960). One of the mechanisms that allows for coexistence of species, with an overlapping fundamental niche, is separation along one or more niche axis (Amarasekare 2003). We hypothesize, as the two species seemingly coexist at Agger Tange, that there has been a spatial and/or temporal displacement of M. edulis. Because C. gigas is competitively superior, M. edulis therefore needs to use an alternative spatial/temporal niche in order for the two species to coexist. We tested this hypothesis by using a thorough analysis of the spatial and temporal distribution of the two species.

Study site
The study site is a sheltered 12,000 m 2 cluster of intertidal bivalve beds situated at Agger Tange in Nissum Bredning (Limfjord, Denmark, N 56°43.3′; E 8°15.4′) (Figure 1), which is primarily inhabited by M. edulis and C. gigas. The study site is surrounded by extensive sand and mud flats that are intersected by deep (~2 m) and narrow tidal channels. The average water depth of the bivalve bed (0.16 ± 0.15 m), combined with a tidal peakto-peak amplitude of 0.25 ± 0.03 m, leaves up to 40% the bivalve bed exposed to air at low tide (unpublished data).

Sample design and data collection
The field aspects of this study were conducted during June 2010. All bivalve samples were collected by hand. A 40 × 40 m virtual grid was laid over a map of the study site using GIS software (ESRI ArcMap, version 10.1, California), and four sample points were randomly generated within each grid cell. In addition, three smaller fractions of the bivalve bed were sampled to capture small-scale changes. Sample points were selected by randomly assigning a direction and distance every fifth meter in the longitudinal direction of the three fractions of the bivalve bed ( Figure 2). However, the two different sample designs gave similar results and data were pooled (n = 155). No samples were collected in the northernmost part of the bivalve bed because it was being used for another experiment. No sample point was sampled twice because the position of each sampling point was determined using GPS that had ±3m resolution.
At each sample point, a 0.25 m 2 ring was placed on the bottom, and all mussels and oysters within the ring were collected and counted. The size distribution was determined in the field by  measuring the shell length of all C. gigas and a minimum of ~100 M. edulis from each sampling point. Shell length was measured from umbo to the longest diameter (nearest mm) using an electronic calliper. Modes in the frequency of shell length measurements were assumed to correspond to cohorts or age classes and were determined using Bhattacharya's Method in the FiSAT II software (FAO, Rome). The method graphically separates normally-distributed groups from a mixed group by using the natural logarithm on ratios of frequencies (Goonetilleke and Sivasubramaniam 1987).
A total of 180 C. gigas (Shell length range: 29 mm to 162 mm) and 193 M. edulis (Shell length range: 14 mm to 68 mm) were randomly picked from the samples (48 and 55 sample points, respectively) to determine condition (CI) and the relationship between shell length and weight. The CI used for further analysis was based only on the sampling points were CI had been determined. All individuals were given a unique identification number and stored at -18 °C until being processed in the laboratory a few weeks later.
After thawing the soft tissue was separated from the shells, and the shell lenghts were measured. For C. gigas only, the shells were cleaned and rinsed thoroughly. The flesh of both bivalve species, and shells of C. gigas, were dried at 105 °C until reaching a constant weight (after approximately 24 and 48 hours, respectively). The condition index (CI) of C. gigas was determined as the ratio between shell-free dry weight (SFDW in g) and dry weight of the shell (g) (Lucas and Beninger 1985). The CI for M. edulis (Petersen et al. 2004) was determined as the ratio between shell-free dry weight (mg) and shell length (cm): SFDW/L 3 Figure 3. Orthophotos depicting the extent of the bivalve bed from 1999 to 2010. In all orthophotos, it is possible to distinguish the bivalve bed from the surrounding sand and mud flats, thus they were included in the analysis of spatial development of the bivalve bed (Danish Digital Orthophoto, consultancy firm COWI, Denmark). The resolution of the orthophotos ranged from 13 × 13 cm to 40 × 40 cm.
The relationship between shell lengths (cm) and weight (SFDW in g) was estimated using a power function fitted to measurements: SFDW C. gigas = 0.0035L 2.634 , R 2 = 0.89 SFDW M. edulis = 0.01211L 2.134 , R 2 = 0.86 Water depth was measured at each sampling point. To correct for tidal fluctuations, readings on a vertical positioned ruler were combined with pressure changes. Pressure was measured by two pressure gauge DST CTD loggers (Star Oddi, Iceland), which were located next to the ruler, just above the bivalves, at the centre of the bivalve bed. Pressure changes (bar) were converted to changes in water level (cm) by using a linear relation between pressure and water level (cm = 893.48Bar -381.84, R 2 = 0.93). This measure for tidal fluctuations had no reference point. It was assumed that the water level at Agger Tange fluctuated around the same mean as the water level at Thyborøn harbour ~4 km away (provided by the Danish Maritime Safety Administration), which is related to the Danish Vertical Reference 90 (DVR90). This was used as a reference point for the water level measurements at Agger Tange by fitting the curve for water level changes to the curve for water level changes from Thyborøn. The depth measurements in every sampling point were corrected according to the fitted curve for tidal fluctuations.
The average CI, used as an indicator of the integrated living conditions, were compared with vertical position.

Spatial development
The spatial development of bivalve beds was examined using orthophotos from 1999, 2004, 2006, 2008 Orthophoto, consultancy firm COWI, Denmark) ( Figure 3). The orthophotos differed in quality and resolution, thus all were assessed for the possibility to discriminate between the bivalve bed and the surrounding sand and mud flats. The outline of the bivalve bed was sketched using a fixed zoom (1:60). Each sketch was confirmed by two independent observers to obtain the most reliable result. Areas smaller than one m 2 were not included in the analysis. Three equally sized zones (South, Middle, and North) were assigned to the sketched orthophotos (Figure 1), so that changes in the North-South direction could be identified. Together, the three zones covered the total extent of the study area.
Fragmentation of the mussel bed was examined using the ratio between area and perimeter of each fragment (> 1m 2 ). The stability of the bivalve bed was evaluated by creating a raster layer (1×1m), based on the coverage of the bivalve bed, and attributing cells with a value of 1 if there were no differences in the cover from year to year. Thus, cells could score a maximum value of 5 (stable) if it was covered in all the years studied (1999, 2002, 2004, 2006, 2008 and 2010), and score a minimum value of 0 (unstable) if it was only covered one year. The effect of stability of the bivalve bed on population structure and condition was examined by extracting sample points that were collected in areas that scored 0 to 1 (unstable) or 4 to 5 (stable) (n = 8 and n = 19, respectively).
The biomass of the two bivalves was used to map their distribution. The area specific shellfree biomass (g SFDW m -2 ) was calculated by summing up individual biomass in each sample point, determined by using shell length measurements and the length-weight regression and, for M. edulis, by multiplying by the fraction of the sample that was measured. Contours of the distribution of bivalve biomass across the entire bivalve bed were estimated by spatial interpolation using inverse distance weighting. This method assumes that bivalve distribution is continuous across the bivalve bed and that local influence diminishes with distance while allowing local observations to influence local estimations. The resulting biomass distribution was produced in a grid of a one meter resolution. The interpolation and all spatial analyses were done using standard GIS software (ESRI ArcMap, version 10.1, California).
The horizontal pattern of relative distribution of the two species on the bivalve bed was estimated in each grid cell (1×1 meter) by ranking area specific biomass from 0 to 1, with 1 representing maximum biomass. A normalized index was then calculated by subtracting the ranked area specific biomass of M. edulis from C. gigas, divided by the combined area specific biomass of the two bivalves (from here on denoted as a distribution index). The distribution index ranged from -1 to 1. Values close to zero indicated an equal mean area-specific biomass of the two bivalves. Positive values (>1.96 SD) indicated dominance by C. gigas and negative values (< -1.96 SD) dominance of M. edulis. Sample points within areas dominated Negative values indicate areas that were exposed to air at average water level in the study period (June 2010). by either C. gigas or M. edulis were extracted for analysis of CI and shell lengths. Shell lengths distributions were be used for determining the ages of C. gigas, and thereby the population structure (for more details see Holm et al. 2015). The population structure of C. gigas was then used as a proxy for the age of different parts of the bivalve bed.

Statistical analysis
Statistical analyses were performed using GraphPad software (GraphPad Prism 6, California). The significance level for all tests was set at α = 0.05. The influence of bed dynamics on shell lengths was tested using a nonparametric Mann-Whitney U-test. The effect of species dominance in different areas on shell length was also tested using a Mann-Whitney U-test. Spearman rank correlation was used to test the effect of vertical position by comparing CI with water depth. Pearson Correlation was used for testing the average condition of M. edulis as a function of total biomass of C. gigas. All mean values were presented ± 1 SD.

Results
There was no difference in the vertical distribution patterns of biomass between C. gigas and M. edulis (Figure 4). Both species were present within the same depth range (-10 cm to 60 cm), with decreasing biomass in the deepest and the shallowest parts of the bivalve bed (Figure 4). The contours of biomass show that C. gigas had the highest biomass in the southern part of the bivalve bed, and was only present with a low biomass at the northern part ( Figure 5A). M. edulis, on the other hand, had high biomass in areas all over the bivalve bed ( Figure 5B).
The CI of C. gigas and M. edulis did not change with depth ( Figure 6 A and B, Spearman rank correlation, r s = 0.004, N = 48, p > 0.05 and r s = -0.051, N = 55, p > 0.05, respectively). Furthermore, there was no significant correlation between the average CI of M. edulis as a function of total estimated biomass of C. gigas (Figure 7, Pearson Correlation, r = -0.091, N = 55, p > 0.05).
There were only minor differences in the total cover of the bivalve bed between years, and it only  increased 5% between 1999 and 2010 (11,600 m 2 to 12,200 m 2 , respectively). However, there were internal differences between the three zones of the bivalve bed. The northern and southern zones increased by 90% and 10%, respectively, whereas the middle zone decreased by 17%. The northern part of the bivalve bed was more fragmented (composed of small clusters of bivalves) than the southern part (average area/perimeter ratio 0.59 ± 0.46 compared to 1.06 ± 0.95, respectively). The northern part of the bivalve bed was primarily composed of mats of M. edulis on top of soft sediment that was easily re-suspended. The southern part of the bivalve bed was characterized by being a solid structure, where C. gigas and M. edulis had settled on a solid reef structure that was composed of thick and compact sediment in which C. gigas and M. edulis shells were embedded.
Labile areas were primarily located in the northern part of the bivalve bed; while the southern part was rather stable (Figure 8). The bed dynamics had no influence on the CI of C. gigas (ANOVA p < 0.05) or M. edulis (p < 0.05) ( Table 1). C. gigas was dominant in relatively few areas, primarily in the southern part, while M. edulis dominated the northern part of the bivalve bed. Furthermore, the relative distribution of the two bivalves was the same in the middle and southern parts of the bivalve bed (Figure 9). There were differences in the size structure of C. gigas between areas dominated by C. gigas and areas dominated by M. edulis. In areas dominated by C. gigas, two size cohorts were present (modes at 50 ± 11 mm and 125 ± 17 mm) (N = 100), whereas, only one size cohort (mode at 40 ± 7 mm) (N = 21) was present in areas dominated by M. edulis. However, there was no significant effect on the shell length of M. edulis, when compared to areas dominated by either C. gigas or M. edulis (Mann-Whitney U-test, U = 8, N 1 = 6, N 2 = 5, p > 0.05) ( Table 1).

Discussion
For M. edulis there are potential trade-offs with regards to being situated in a bivalve bed inhabited by C. gigas. On the one hand, the two bivalves are potential competitors for resources, such as space and food (Troost 2010); however, the complex reef structure of oyster reefs provides mussels with protection from predation from, for example, shore crabs (Carcinus maenas) (Waser et al. 2015). In some bivalve beds in the Wadden Sea, M. edulis occupies the crevices between individuals of C. gigas (Eschweiler and Christensen 2011). This micro-distribution seems to be a tradeoff between food uptake and predator avoidance, as the CI and predation pressure are both reduced when mussels are situated between individuals of C. gigas (Eschweiler and Christensen 2011). In the Wadden Sea, strong recruitment of M. edulis occurs after harsh winters (Beukema 1992;Beukema et al. 2001), presumably because this creates a mismatch between mussel recruits and their primary invertebrate predators (Strasser and Günther 2001;Strasser 2002;Beukema and Dekker 2005). Thus, following severe winters, settling mussel spat experience reduced predation pressure, which allows for this stronger than average recruitment success. This mechanism does not seem to be present at Agger Tange, as there was not strong recruitment of M. edulis after the harsh winter of 2009/2010 (Holm et al. 2015). Thus, avoiding invertebrate predators seems to be of less importance in our study area.
The bivalve bed at Agger Tange is located in a shallow intertidal area where up to 40% of the bed is exposed at low tide (unpublished data), and they therefore seem to be food limited, as their filtration capacity exceed food availability (Vismann et al., unpublished data). Vertical position could therefore have a significant effect on food availability, with the uppermost parts of the raised bed being unfavourable, as C. gigas and M. edulis in these areas would be unable to take up sufficient amount of food. In terms of maximum clearance rates and growth rates C. gigas seem to be competitively superior (Troost 2010). While the two species may not directly compete for the same food source (Bougrier et al. 1997;Dubois et al. 2007), they retain the same size spectrum of particles, and therefore reduce the food availability for the other species (Troost 2010). In an environment with low food availability, competition between C. gigas and M. edulis is expected to shape their population sizes (Troost 2010). If the two species directly competed for food at Agger Tange, and C. gigas being the competitively strongest species, then it is expected that M. edulis would be displaced to a less favourable location on the bivalve bed, have lower body condition, or have smaller body sizes. However, there were no detectable vertical separation between C. gigas and M. edulis and the CI was not related to vertical position. Furthermore, the CI of M. edulis was unaffected by the biomass of C. gigas, and thus C. gigas seem unable to supress the food uptake of M. edulis.
The overall coverage of the bivalve bed changed little from 1999 to 2010. However, coverage has been influenced by significant internal dynamics, as some parts expanded and some retracted during the studied period without changing the total coverage of the bivalve bed. The northern part became more fragmented, and was composed of mussel mats on top of soft sediment, and the coverage changed significantly from year to year. In the south, the bed structure was raised and compact, due to the accumulation of shells and sediment, and changed little in coverage from year to year. There was, however, no evidence that interaction between the two species could explain this difference, as there were no differences in the CI of M. edulis as a function of the biomass of C. gigas. Furthermore, there were no differences in the average shell lengths or CI between individuals of both species situated in the north or south of the bivalve bed. Hence, the differences between these two parts of the bivalve bed seem to be best explained by a physical stressor occurring along the bivalve bed or that the age of these two areas on the bivalve bed differs. The study area is frequently covered with ice during winter, thus mechanical damage from ice moving back and forth could destroy parts of the bivalve bed. However, the northern part is also one of the deepest parts of the bivalve bed, and it therefore seems unlikely that ice alone could explain the difference. Furthermore, Strand et al. (2012) showed that after the severe winter of 2009/2010 there was only an average mortality of approximately 18% at Agger Tange, thus it seems unlikely that harsh winters with ice cover should govern the cover of the bivalve bed. There were no significant differences in shell lengths between the north and south, which would be expected if predation should regulate this area. It is, however, characteristic for old and stable bivalve beds that they are composed of compact consolidated sediment, clearly recognisable and raised above the adjacent sea floor (Reise 2002;Hertweck and Liebezeit 2002;de Vlas et al. 2005). In contrast, newly formed mussel beds are often temporary and consist of a monolayer of mussels (Dankers et al. 2001). Hence, the southern part of the bivalve bed could reflect a stable and older habitat, resistant to physical disturbances, and therefore has changed little in terms of coverage by bivalves during the studied period. In contrast, the northern part could be young as it is composed of labile mussel mats on top of soft sediment. C. gigas only dominated areas in the southern part of the bivalve bed, which was also determined as stable. In these areas, two size cohorts of C. gigas were present (50 ± 11 mm and 125 ± 17 mm), whereas only one (40 ± 7 mm) was present in areas dominated by M. edulis. Holm et al. (2015) estimated individuals with a shell length of 50 ± 13 mm and 57 ± 9 mm, based on yearly shell increment and growth rings, to be one year old, and, based on growth rings, individuals with a shell length of 85 ± 13 mm and 162 ± 13 mm to be three and five years old, respectively. The two smallest size cohorts (40 and 50 mm) were therefore estimated to be one year old, and the largest three to five years old. The population structure of C. gigas between areas dominated by one or the other species therefore differs.
Ecosystem engineers often produce a positive feedback when they modify the habitat (Jones et al. 1997;Hui et al. 2004;Wright and Jones 2006). In the south, the raised bed is expected to increase turbulence, which increases food availability for filter feeding bivalves (Lenihan 1999). However, there was no difference with respect to the average CI of C. gigas and M. edulis between labile and stable areas (north and south, respectively), indicating no difference in bivalve feeding conditions; therefore, food availability seems unable to explain why there are no old C. gigas in the labile areas, which is also the areas where M. edulis is dominant. The difference could be explained by the process where M. edulis is the first species to settle and subsequently provide substrate for C. gigas to settle upon. In the Oosterschelde, the Netherlands, C. gigas colonise bare mudflats (Troost 2010), which they stabilise, as they cement to conspecifics and thus creates solid structures, which then strengthen as shell debris and sediment consolidates (Reise et al. 2008). Thus, it has the ability to further stabilise the bivalve bed, after M. edulis has settled. It is argued that C. gigas might stabilise the sediment on a longer time scale than M. edulis (Troost 2010). This process, where M. edulis is responsible for the first step in an expansion of the bivalve bed, and the subsequent stabilisation by C. gigas, induces a temporal separation between the two species. There were no differences in the CI of C. gigas and M. edulis between old and newly established areas and no vertical separation between the two bivalves. Thus, temporal and spatial separation of the C. gigas and M. edulis seem unable to explain why the two species seemingly coexist at this sheltered intertidal bivalve bed.