Seagrass of Vasiliko Bay, Eastern Mediterranean: Lost Cause or Priority Conservation Habitat?

: Mediterranean coasts are a ﬀ ected by multiple mounting pressures. In Cyprus, marine ﬁsh farming has grown rapidly in the past decade and is concentrated in the west side of Vasiliko Bay. The east coast of this bay has ports, a power station, a desalination unit, a cement factory, a major new oil terminal, and gas storage facilities. The bay is earmarked to create the largest hydrocarbon processing, storing, and transport facility in the region. Here, we assess the status of Posidonia oceanica habitat in an understudied region at the upper thermal, and eastern limit, of this Mediterranean endemic seagrass. An extensive ancient seagrass meadow was revealed, covering about 200 ha across 10 km of coastline, over soft substrata at ca 10–30 m depth, and over hard substrata at ca 0–6 m depth. Seagrass shoot density and leaf surface area decreased, both with increasing depth and with proximity to industrial developments; part of the meadow had been destroyed by dredging to build a jetty. Close to ﬁsh farms the seagrass had higher epiphytic biomass as well as lower leaf number, mass, and surface area, all of which indicate adverse e ﬀ ects of eutrophication and increased turbidity. Despite these multiple stressors, most of the meadow was in good ecological status, with some of the highest shoot densities ever reported. Furthermore, iconic species like sea turtles, monk seals, and dolphins were seen during sampling. Posidonia oceanica meadows o ﬀ Cyprus are among the most valuable in the Mediterranean due to their tolerance of high seawater temperatures. Managers of future coastal developments in the region will need to adhere to European legislation and international conventions designed to secure the socioeconomic beneﬁts of seagrass beds.


Introduction
Seagrass meadows provide valuable ecosystem services but are being lost at 7% year −1 globally [1], with water quality and other stressors placing seagrass meadows among the most threatened ecosystems on Earth [2].Seagrass and algal beds cover less than 0.5% of the global ocean seabed but account for over 5% of the world s ecosystem services [3].For example, seagrasses help mitigate anthropogenic CO 2 emissions since they act as significant carbon sinks [4].
The endemic Mediterranean seagrass Posidonia oceanica (Linnaeus) Delile 1813, has the highest carbon burial rates, and forms the largest organic carbon stocks, of all seagrass species [5,6].It forms meadows that extend to depths of 45 m in the clearest waters; they cover a known area of 14,417 km 2 , while their presence along much of the southern Mediterranean and Levantine basin, remains uncharted [7,8].Posidonia oceanica stabilise sediments and help protect coastlines from erosion.They also sequester nutrients and contaminants, oxygenate the water with up to 14 Lm −2 d −1 , transfer nutrients to food webs, and provide an essential habitat for many protected and commercial species, contributing to fisheries and cultural heritage [9,10].The canopy and rhizomes of this habitat can be highly biodiverse with 660 species of epibionts reported in the literature [11], as well as 394 polychaete species, 120 crustacean spp., 171 mollusc spp., and 122 fish spp. reported from P. oceanica beds in the eastern Mediterranean [12][13][14][15].Additionally recorded in the meadows are many iconic and protected species, such as the large bivalve Pinna nobilis Linnaeus, 1758, the giant gastropod Tonna galea (Linnaeus, 1758), the seahorse Hippocampus hippocampus (Linnaeus, 1758), and the turtle Chelonia mydas (Linnaeus, 1758).
Seagrasses are widely used as bioindicators of ecosystem quality [16].In the Marine Strategy Framework Directive (MSFD, 2008/56/EC), P. oceanica is used as an indicator of good environmental status.It is on the International Union for the Conservation for Nature 'Red List', as it has suffered an estimated 34% decline in distribution over the past 50 years [7].Significant losses have been recorded around fish farms, where seagrass meadows shift from carbon sinks to carbon sources [17,18] although some meadows have remained stable over time [19,20].Posidonia oceanica can withstand seawater temperature variations, many contaminants, and alien species invasion, but it is highly sensitive to salinity fluctuations, turbidity and increased sedimentation, anchoring, and trawling [21].Habitat loss of this species is considered irreversible on human time-scales [22], because even if optimal conditions return, impacted areas take centuries to be recolonised [23,24].
There is uncertainty about the current distribution of P. oceanica in the eastern Mediterranean [7].It appears to be absent from the coasts of Syria, Lebanon, Israel, eastern Egypt, and eastern Turkey.It is thought that high temperature and salinity, plus turbidity due to coastal disturbance and river discharges, limits its distribution [25].Cyprus is a phenomenal exception to this pattern, despite having summer seawater waters at the upper thermal limit reported for this species [7,25].It ranges from 0 to >40 m depth around the perimeter of the whole island, and is the easternmost reported population of the species.The meadows cover an estimated area of 9040 ha in the 30% of the Cyprus coastline that has been surveyed [7].These P. oceanica meadows seem to be genetically isolated from those in the western and central Mediterranean [26][27][28].Other factors contributing to the enigmatic distribution of P. oceanica around Cyprus may be the occasional summer upwelling events that cool surface waters, the lack of permanent rivers, and the presence of reservoirs that minimise terrestrial run-off, and the relatively low anthropogenic and industrial footprint around the majority of the island's coastline.
Seagrass meadows are threatened by changes in water quality due to pressures such as urbanisation and coastal modification, fish farming, and dredging.Vasiliko Bay is the most industrialised coastline of Cyprus, although P. oceanica meadows extend from shallow water to approximately 30 m depth.Fish farms in the western side of the bay operated for almost two decades above seagrass beds, although since 2011 the farms have been located away from the meadows in deeper waters.On the other side of the bay, industrial pollution and dredging to position pipelines, construct or enlarge ports, marine jetties, and berths to accommodate large vessels, pose a growing threat to the meadows.
Here, we present results of the first scientific survey to assess the status of P. oceanica meadows in Vasiliko Bay.We correlate our data with human activities, and assess the ecological status of the bay.We also review whether P. oceanica descriptors are different compared to cooler and less oligotrophic sites around the Mediterranean.

Study Area
Vasiliko Bay is in southern Cyprus approximately 25 km east of the town of Limassol.The national mariculture industry is concentrated on the western side of this bay, with cages about 1-3 km from the shore (Figure 1); they produced 4000 tonnes of farmed fish in 2013, representing about two thirds of the total aquaculture production of Cyprus.Vasiliko Bay is in southern Cyprus approximately 25 km east of the town of Limassol.The national mariculture industry is concentrated on the western side of this bay, with cages about 1-3 km from the shore (Figure 1); they produced 4000 tonnes of farmed fish in 2013, representing about two thirds of the total aquaculture production of Cyprus.The central and eastern coasts of the bay have been entirely modified and heavily industrialised.In the central part of the bay the largest power plant in Cyprus burns hydrocarbons supplied from seabed pipelines.The east part of the bay has the largest cement manufacturing plant on the island, producing around 6 kT of clinker per day.There are also three ports in the bay (Figure 1).Vasiliko port on the eastern tip of the bay is a major importer/exporter of wheat, coal, perlite, cement, gravel, scrap iron, and others.Archirodon is a small port nearby, which includes a dry dock for ship maintenance and repairs, and a small shelter for fish farmers and fishing boats.Furthermore, there are active clay quarries, a recently constructed marine jetty, and the facilities of various multi-national corporations for the storage and management of petroleum and oil products, with a total storage capacity of storage of 858,000 m 3 .Government plans for further development are underway to expand industries in the region, such as a liquefied natural gas plant and associated marine jetties to transfer oil-based products.More dredging is anticipated, over or near seagrass meadows, to construct jetties and lay pipes, expand Vasiliko port, and construct a new port west of Vasiliko Bay to accommodate the fish farmers of the area.These major construction projects have the potential to jeopardise the integrity of local P. oceanica meadows, hence the need to increase awareness about the condition of this habitat and inform coastal zone management decisions.

Seagrass Mapping
We surveyed the distribution of P. oceanica meadows along 10 km of coastline.For the 0-15 m depth zone, IKONOS (1 m resolution, 2011) and GeoEye (0.5 m resolution, 2010) satellite images and higher resolution aerial images (2009-2010) from the Department of Lands and Surveys were compared, and the aerial images that had the greatest resolution were selected to digitise the shallow P. oceanica meadows using ArcGIS software.For mapping > 15 m depth, a Side Scan Sonar C-MAX The central and eastern coasts of the bay have been entirely modified and heavily industrialised.In the central part of the bay the largest power plant in Cyprus burns hydrocarbons supplied from seabed pipelines.The east part of the bay has the largest cement manufacturing plant on the island, producing around 6 kT of clinker per day.There are also three ports in the bay (Figure 1).Vasiliko port on the eastern tip of the bay is a major importer/exporter of wheat, coal, perlite, cement, gravel, scrap iron, and others.Archirodon is a small port nearby, which includes a dry dock for ship maintenance and repairs, and a small shelter for fish farmers and fishing boats.Furthermore, there are active clay quarries, a recently constructed marine jetty, and the facilities of various multi-national corporations for the storage and management of petroleum and oil products, with a total storage capacity of storage of 858,000 m 3 .Government plans for further development are underway to expand industries in the region, such as a liquefied natural gas plant and associated marine jetties to transfer oil-based products.More dredging is anticipated, over or near seagrass meadows, to construct jetties and lay pipes, expand Vasiliko port, and construct a new port west of Vasiliko Bay to accommodate the fish farmers of the area.These major construction projects have the potential to jeopardise the integrity of local P. oceanica meadows, hence the need to increase awareness about the condition of this habitat and inform coastal zone management decisions.

Seagrass Mapping
We surveyed the distribution of P. oceanica meadows along 10 km of coastline.For the 0-15 m depth zone, IKONOS (1 m resolution, 2011) and GeoEye (0.5 m resolution, 2010) satellite images and higher resolution aerial images (2009-2010) from the Department of Lands and Surveys were compared, and the aerial images that had the greatest resolution were selected to digitise the shallow P. oceanica meadows using ArcGIS software.For mapping > 15 m depth, a Side Scan Sonar C-MAX system was used.During 2012-2013, it was towed in transect lines (at 25 or 50 m intervals) parallel to the coastline, with complete overlapping of coverage in each transect.Some small areas in the vicinity of the fish farms with obstacles in the water column (e.g., ropes and chains) were not mapped.Digital data obtained from the towfish were processed in SonarWiz 5 software.Several features (sand, matte, living P. oceanica and rocky reefs) were ground-truthed by SCUBA dives.

Assessment of P. oceanica Meadows
The seagrass was monitored at 30 stations spread across the study area.Ten stations were set up in each depth zone: (i) 5 m, (ii) 15 m, and (iii) at the lowest limits of the meadows, which ranged between 22-33 m.The 15 m and deeper stations were about 1 km apart.In the shallows, the 5 m stations were on seagrass patches that were unevenly distributed.The 5 m and 15 m stations were sampled twice, the deepest sites were sampled once (Table S1 in Supplementary Material).

Shoot Densities
Shoot densities were counted in 0.16 m 2 quadrats, placed 1.5 m from the edge of the meadow investigated.In total, 500 quadrat measurements were recorded; 20 replicate measurements from each of the stations of 5 m and 15 m depth and ten from the lower limit stations to calculate shoot density m −2 .At the lowest limits, the proportion of horizontal (plagiotropic) shoots was determined to assess regression or advancement of the habitat [29,30].

Leaf Morphometrics and Epiphytes
Seagrass morphometric analysis was conducted using the techniques of Giraud [31], on 700 shoots collected using a non-destructive approach [32].Presence of epiphytes, herbivory marks, and necrosis were noted for 20, 40, and 10 orthotropic shoots from each sampling station at 5 m, 15 m, and lower limits, respectively.Leaf surface and the leaf area index (LAI) were calculated for all stations.The LAI expresses the total surface area (m 2 ) of P. oceanica canopy in 1 m 2 of benthic surface (LAI = mean foliar surface area shoot −1 × mean shoot density m −2 ).The LAI values were interpolated using the algorithm 'Nearest Neighbour' ArcGIS (v.10.2.2) software, and a heat map was generated.
The ratio of epiphyte to leaf dry mass was estimated for a total of 400 P. oceanica shoots, collected from the 15 m depth stations.The epiphytes were gently scraped off from the leaves with a scalpel and placed on dry pre-weighted Whatman GF/F.Leaves and epiphytes from each shoot were placed in a drying oven at 60 • C until constant weight.

Ecological Status
Seagrasses have a wide spatial distribution, essential ecological role, and high sensitivity to anthropogenic disturbance.The EU Water Framework Directive (WFD 2000/60/EC) declared the benthic macrophytes (macroalgae and angiosperms) as biological quality elements (BQEs) to assess the ecological status of coastal and transitional waters.The ecological status class (ESC) was assessed with WFD (2000/60/EC) seagrass indices, PREI (P.oceanica Rapid Easy Index) [33] and BiPo [29], applied in both seasons (autumn 2012 and spring 2013) at ten locations, and set at about 1 km intervals.The two indices, PREI and BiPo indices, are multimetric and are based on the deviation of the status of the BQE from its potential status under pristine conditions (reference conditions).The metrics that were used in the calculations included: (i) depth at the lower limits of the meadow, (ii) characterisation of the meadow type at the lower limit, (iii) shoot density at 15 m depth (mean of both seasons), (iv) foliar surface shoot −1 at 15 m depth (seasonal mean), and (v) dry weight of epiphyte/leaf (E/L) at 15 m depth, the latter used in the assessment by PREI only (seasonal mean).The reference conditions used in the calculation of the ESC scores were: (i) those suggested by Lopez et al. [29] for BiPo and applied in the central and western Mediterranean, (ii) the reference values set by the national authority of Cyprus for PREI, following the intercalibration MedGIG exercise, and (iii) the reference conditions adapted for the area of study and the time of sampling, set for the shoot density and foliar surface shoot −1 as the average of the three highest values (after the highest) measured during the survey [34].For E/L dry weight the reference was set at zero, considering that healthy seagrass leaves have minimal epiphytic biomass [33], and for the lower limit, the reference depth was set at 33 m, determined from the mapping with a CM2 C-MAX Side Scan Sonar and confirmation dives (for reference condition values used, see Table S2 in supplementary material).

Statistical Analyses
The spatial differences in the P. oceanica descriptors across sampling stations were identified by comparing stations of the same depth (5 m and 15 m) using a 1-way ANOVA, followed by a Tukey's honestly significant difference (HSD) test.The importance of season and the location of the stations on the development of epiphytes on P. oceanica meadows was examined via a 2-way ANOVA, using the E/L dry weight descriptor.To determine whether the two biotic indices (BiPo and PREI) followed the same pattern across stations, averaged EQR seasonal results were compared using a Pearson's Correlation analysis.For all the parametric tests, the assumptions for normality and equal variances were verified using a Shapiro-Wilk test and Bartlett test, respectively.When assumptions were violated, log 10 or square root transformations were conducted, depending on the skewness of the residuals.In the case of 1-way ANOVA analyses, when assumptions were not met, the analysis proceeded with the non-parametric Kruskal-Wallis rank sum test, followed with a post-hoc Dunn's test, with Bonferroni correction, to locate the pairwise differences.For all the statistical analyses the significance level α was set at 0.05, and their computation was carried out by R-studio (v 1.0.153).For the Dunn's test, the "Dunn.test"package was used [35], and graphical output was generated via the package ggplot2 [36].
Redundancy (ordination) analysis (RDA) was further carried out to identify trends and general patterns that could explain the variation of the P. oceanica descriptors (shoot density, number of leaves, foliar surface, E/L biomass, herbivory, and coefficient (A)), observed across the stations against proxies of anthropogenic pressures, depth, and time.The analysis consisted of three runs, where the first included the P. oceanica descriptors collected from all the sampling stations (shallow, intermediate, and lower limits).This analysis did not include time, since lower limits were sampled only once.The second run included only the shallow stations (5 m), and the last run only the stations from the intermediate depth (15 m).The proxies of anthropogenic pressures comprised: (1) number of anthropogenic pressures within 1 km radius from each station, (2) distance to the nearest anthropogenic pressure, (3) direction from the nearest anthropogenic pressure (excluded from the first part of the analysis), and (4) geographical position within the bay (bay location), defined by regions of anthropogenic pressure gradient.In the case of the shallow stations, the pressure gradient was defined as: outer bay to inner bay (low to high pressure), where the latter has the highest cumulative land-based anthropogenic pressures; intermediate depth, as western to eastern bay (high to low pressure), where the western side hosts almost all of the fish farm units; the whole study area, as western, central, and eastern bay.The RDA was conducted using the "vegan" package [37] in R-studio (version 3.5.0).Prior to computation, all P. oceanica descriptors were log (1 + x) transformed.When a right angle between a P. oceanica descriptor and predictor (arrow) was formed on the RDA plot, a zero correlation was assumed, and when it was lower than 20 • , a strong correlation was implied.If a descriptor was found at ~180 • from the predictor, it was then interpreted as a strong negative correlation.RDA was further complemented with two permutation ANOVA tests including: (1) a test to identify the significance of all the predictors together, (2) type III test (permutation = 500), to analyse the marginal effects when each predictor was eliminated from the model while the others remained.This allowed the identification of those that had a significant role in the variance of the model.

Mapping and Visual Observations
We mapped an extensive seagrass meadow covering about 192 ha, along the 10 km of coastline surveyed (Figure 1).About three quarters of this seagrass canopy was in the 10-30 m depth zone.The lower limit of the seagrass meadow ranged from 22 to 33 m.The shallowest lower limits were near Vasiliko port and near a marine jetty, which was recently constructed.In the shallows, the seagrass was found on limestone, which is the dominant substratum in the intertidal and upper subtidal, and extends up to 250 m from the coastline to 6 m depth.
In the shallow and some of the intermediate depth sites, the seagrass had a dense canopy with a few epibionts.A well-developed diverse coralligenous community, dominated by rhodophytes and bryozoans was pronounced at some sites (e.g., shallow station 2), but absent from the sites near industrial developments.The seagrass had dense shoots and long leaves with few epiphytes at the 5 m stations.At the western sites of intermediate depth and lower limits, near the fish farms, the meadows had significantly higher epiphytic vegetation and fine particulates covering the canopy (Figure 2).At the intermediate 15 m depth sites, the seagrass formed an almost continuous meadow, and erosive matte of >3 m vertical height was often observed.The dead matte of the seagrass meadows was more common at some intermediate, and all deeper, stations around the perimeter of fish farms and the Vasiliko port, and was usually covered by the chlorophytes Caulerpa prolifera (Forsskål) J.V.Lamouroux and C. cylindracea Sonder and turf algae.lower limit of the seagrass meadow ranged from 22 to 33 m.The shallowest lower limits were near Vasiliko port and near a marine jetty, which was recently constructed.In the shallows, the seagrass was found on limestone, which is the dominant substratum in the intertidal and upper subtidal, and extends up to 250 m from the coastline to 6 m depth.
In the shallow and some of the intermediate depth sites, the seagrass had a dense canopy with a few epibionts.A well-developed diverse coralligenous community, dominated by rhodophytes and bryozoans was pronounced at some sites (e.g., shallow station 2), but absent from the sites near industrial developments.The seagrass had dense shoots and long leaves with few epiphytes at the 5 m stations.At the western sites of intermediate depth and lower limits, near the fish farms, the meadows had significantly higher epiphytic vegetation and fine particulates covering the canopy (Figure 2).At the intermediate 15 m depth sites, the seagrass formed an almost continuous meadow, and erosive matte of >3 m vertical height was often observed.The dead matte of the seagrass meadows was more common at some intermediate, and all deeper, stations around the perimeter of fish farms and the Vasiliko port, and was usually covered by the chlorophytes Caulerpa prolifera (Forsskål) J.V.Lamouroux and C. cylindracea Sonder and turf algae.Dredging to construct a marine jetty in the eastern part of Vasiliko Bay occured between the two annual sampling periods, and severe impacts were observed on the underlying meadows e.g., sampling stations 8 and 9, at 15 m depth (Figure 3).Dredging to construct a marine jetty in the eastern part of Vasiliko Bay occured between the two annual sampling periods, and severe impacts were observed on the underlying meadows e.g., sampling stations 8 and 9, at 15 m depth (Figure 3).Dredging to construct a marine jetty in the eastern part of Vasiliko Bay occured between the two annual sampling periods, and severe impacts were observed on the underlying meadows e.g., sampling stations 8 and 9, at 15 m depth (Figure 3).

Posidonia oceanica Structural Descriptors
The mean values ± SE of all P. oceanica structural descriptors measured at the 30 sampling stations are presented in Table S3 in Supplementary Material.

Posidonia oceanica Structural Descriptors
The mean values ± SE of all P. oceanica structural descriptors measured at the 30 sampling stations are presented in Table S3 in Supplementary Material.
Large variation in shoot densities was measured among the 5 m depth stations (Kruskal-Wallis rank sum test, χ 2 = 61.84,df = 9, p < 0.05).The highest mean densities, exceeding 900 shoots m −2 , were reported at the sampling stations outside the bay (Figure 4).At the 15 m depth stations, the mean shoot density ranged between 395 and 478 shoots m −2 , and varied across the stations (1-way ANOVA, F = 2.88, df = 9, p < 0.05).Measuring the shoot density at station 9 in spring 2013 was not possible because the P. oceanica meadow studied in the autumn of 2012 had been dredged for the construction of the marine jetty.The mean shoot density at the lower limit stations ranged between 116 and 345 shoots m −2 (Table S3 in Supplementary Material).A low ratio of plagiotropic rhizomes (<30%) was measured at the sampling sites in the deeper stations.No comparisons among the lower limit stations were performed because the depth and light variable was not constant.
Both 5 m and 15 m depth stations showed inter-variation in foliar surface per shoot (1-way ANOVA, F5m = 3.74, F15m = 8.11, df = 9, p < 0.05).At the shallow sampling stations (5 m), the mean foliar surface per shoot ranged from 213-313 cm 2 , whereas at the sampling stations of intermediate depth (15 m) it ranged from 147-225 cm 2 (Figure 4).At the lower limits, mean foliar surface area values ranged from 102 to 162 cm 2 (Table S3 in Supplementary Material).About ten times higher LAI values (m 2 of seagrass canopy in 1 m 2 ) were recorded at shallow sampling stations compared to the lower limits (Figure 5).The highest values (>21-27 m 2 canopy m −2 ) were reported at the 5 m depth stations outside the Vasiliko Bay area.At the 15 m depth stations mean LAI ranged from 5.86-10.01m 2 canopy m −2 .At this depth the lowest values were measured at stations 4 and 5, which are nearest to fish farms.The mean LAI values at the lower limit stations ranged from 1.46 to 4.90 m 2 canopy m −2 (Table S3 in Supplementary Material).
About ten times higher LAI values (m 2 of seagrass canopy in 1 m 2 ) were recorded at shallow sampling stations compared to the lower limits (Figure 5).The highest values (>21-27 m 2 canopy m −2 ) were reported at the 5 m depth stations outside the Vasiliko Bay area.At the 15 m depth stations mean LAI ranged from 5.86-10.01m 2 canopy m −2 .At this depth the lowest values were measured at stations 4 and 5, which are nearest to fish farms.The mean LAI values at the lower limit stations ranged from 1.46 to 4.90 m 2 canopy m −2 (Table S3 in Supplementary Material).S4 in Supplementary Material).Confirming the visual observations, E/L varied across the Vasiliko Bay stations, both in spring (Kruskal-Wallis, χ 2 = 125.67,df = 9, p < 0.001) and summer (Kruskal-Wallis, χ 2 = 139.40,df = 9, p < 0.001) (Figure 6).There was a significant interaction between the main effects (station location and season) (2-way ANOVA, F = 11.87,df = 1, p < 0.05).Several stations had higher E/L in the western side nearest to fish farms and in spring (overall mean epiphytic biomass 0.124 g shoot −1 ) than in the autumn (overall mean epiphytic biomass 0.060 g shoot −1 ) (Figure 6).S4 in Supplementary Material).Confirming the visual observations, E/L varied across the Vasiliko Bay stations, both in spring (Kruskal-Wallis, χ 2 = 125.67,df = 9, p < 0.001) and summer (Kruskal-Wallis, χ 2 = 139.40,df = 9, p < 0.001) (Figure 6).There was a significant interaction between the main effects (station location and season) (2-way ANOVA, F = 11.87,df = 1, p < 0.05).Several stations had higher E/L in the western side nearest to fish farms and in spring (overall mean epiphytic biomass 0.124 g shoot −1 ) than in the autumn (overall mean epiphytic biomass 0.060 g shoot −1 ) (Figure 6).

Main Drivers Affecting P. oceanica Descriptors
Vasiliko Bay shows heterogeneous responses of P. oceanica descriptors attributed to both abiotic variation and anthropogenic pressure.The major driver within the whole area of study was depth, as it explained most of the variation in the dataset (F = 96.7,p < 0.01) and had a strong negative

Main Drivers Affecting P. oceanica Descriptors
Vasiliko Bay shows heterogeneous responses of P. oceanica descriptors attributed to both abiotic variation and anthropogenic pressure.The major driver within the whole area of study was depth, as it explained most of the variation in the dataset (F = 96.7,p < 0.01) and had a strong negative correlation with shoot density and foliar area (Figure 7).No clear patterns were detected with number of leaves, adult leaves that had herbivory marks or that had a broken apex (coefficient A), although they might be influenced by a combination of factors including distance to the nearest pressure, depth, and the pressure gradient (bay location) (Figure 7).At shallow stations (5 m), the most likely controlling factors were the direction of the station from the anthropogenic pressure (F = 8.54, p < 0.01), time (F = 16.87,p < 0.01), and bay location (F = 22.36, p < 0.001).The first factors affected foliar surface area, which increased from autumn 2012 to the summer of 2013, but was lower at stations located southwesterly from the industrialised shorefront (Figure 7).Moreover, the bay location seemed to influence the shoot density substantially.All stations located within the bay were characterised by lower shoot densities than the outer bay stations of the same depth that have little anthropogenic influence (Figures 4 and 7).
At sampling stations of intermediate depth (15 m) the variation of P. oceanica descriptors was driven by the cumulative effect of four out of five variables tested: time (F = 4.26, p < 0.05), bay location (F = 16.06,p < 0.01), direction (F = 42.71,p < 0.01), and number of pressures (F = 3.6, p < 0.05).The bay location and time (autumn 2012 and spring 2013) affected the E/L, whereby higher E/L values were reported, in the western side of the bay near the mariculture facilities, and during the spring of 2013 (Figures 6 and 7).At 15 m depth stations the number of leaves per shoot was affected mostly by the anthropogenic pressure gradient, and to a lesser extent by the number of pressures, showing lower values near fish farms (Figures 4 and 7).Shoot densities were not strongly affected by any variable (Figures 4 and 7).Foliar area strongly correlated with direction from the nearest pressure source, and had lower values in stations located inshore/northwards from an anthropogenic pressure.

Ecological Status Class
Using the metrics determined in this study (Table S5 in Supplementary Material) and the site- At shallow stations (5 m), the most likely controlling factors were the direction of the station from the anthropogenic pressure (F = 8.54, p < 0.01), time (F = 16.87,p < 0.01), and bay location (F = 22.36, p < 0.001).The first factors affected foliar surface area, which increased from autumn 2012 to the summer of 2013, but was lower at stations located southwesterly from the industrialised shorefront (Figure 7).Moreover, the bay location seemed to influence the shoot density substantially.All stations located within the bay were characterised by lower shoot densities than the outer bay stations of the same depth that have little anthropogenic influence (Figures 4 and 7).
At sampling stations of intermediate depth (15 m) the variation of P. oceanica descriptors was driven by the cumulative effect of four out of five variables tested: time (F = 4.26, p < 0.05), bay location (F = 16.06,p < 0.01), direction (F = 42.71,p < 0.01), and number of pressures (F = 3.6, p < 0.05).The bay location and time (autumn 2012 and spring 2013) affected the E/L, whereby higher E/L values were reported, in the western side of the bay near the mariculture facilities, and during the spring of 2013 (Figures 6 and 7).At 15 m depth stations the number of leaves per shoot was affected mostly by the anthropogenic pressure gradient, and to a lesser extent by the number of pressures, showing lower values near fish farms (Figures 4 and 7).Shoot densities were not strongly affected by any variable (Figures 4 and 7).Foliar area strongly correlated with direction from the nearest pressure source, and had lower values in stations located inshore/northwards from an anthropogenic pressure.

Ecological Status Class
Using the metrics determined in this study (Table S5 in Supplementary Material) and the site-adapted reference values (Table S2 in Supplementary Material), the ESC of the Vasiliko basin was estimated to be mostly 'Good' with both PREI and BiPo WFD indices.Exceptions were the sampling stations 5 (near a fish farm) and 9 (near the marine jetty) in the spring of 2013, which were estimated to have 'Medium' ESC with the BiPo index.The ESC assessed at the sampling stations 7 (in the central area of the bay) and 10 (east of Vasiliko bay, near Zygi) was assessed as marginally 'High' with the PREI (Figure 8).In both seasons, the ESC assessed with the PREI had slightly higher (0.064 ± 0.029) EQR values than the BiPo index, while the seasonal fluctuations in the EQR values were small (0.028 ± 0.024) for both indices (Figure 8).The EQR values between the two biotic indices were highly positively correlated (Pearson′s Correlation, ρ = 0.93, p < 0.05).The EQR scores using the site-adapted reference values were similar but consistently slightly lower (0.036 ± 0.026) than those assessed against the intercalibrated reference values for PREI, and slightly higher (0.029 ± 0.006) than those assessed against the reference values of BiPo, but with either method the evaluation of the ESC was the same.

Discussion
We mapped an impressive ancient Posidonia oceanica meadow in Vasiliko Bay (Cyprus) and, although it had been damaged in some areas, we hope this research is used to assist with responsible management to secure the ecosystem services that valuable seagrass habitats provide.
In the shallows P. oceanica grew in patches on rock, as is often the case in high-energy conditions [38,39].It had the highest densities and the longest leaves in shallow water, thanks to more light and increased productivity [40,41].This is surprising for Cyprus though, as waters above the summer thermocline reach 28.4 °C, which is the known upper thermal tolerance of this species [25].Cypriot P. oceanica beds might be an especially important genetic strain, as sea surface temperatures are rising rapidly in the region due to greenhouse gas emissions.Severe P. oceanica habitat loss (70%) is projected for 2050, with potential functional extinction by 2100, unless thermally resistant strains can be found [42].Cypriot strains could prove useful as transplants or as a source of seeds in other In both seasons, the ESC assessed with the PREI had slightly higher (0.064 ± 0.029) EQR values than the BiPo index, while the seasonal fluctuations in the EQR values were small (0.028 ± 0.024) for both indices (Figure 8).The EQR values between the two biotic indices were highly positively correlated (Pearson s Correlation, ρ = 0.93, p < 0.05).The EQR scores using the site-adapted reference values were similar but consistently slightly lower (0.036 ± 0.026) than those assessed against the intercalibrated reference values for PREI, and slightly higher (0.029 ± 0.006) than those assessed against the reference values of BiPo, but with either method the evaluation of the ESC was the same.

Discussion
We mapped an impressive ancient Posidonia oceanica meadow in Vasiliko Bay (Cyprus) and, although it had been damaged in some areas, we hope this research is used to assist with responsible management to secure the ecosystem services that valuable seagrass habitats provide.
In the shallows P. oceanica grew in patches on rock, as is often the case in high-energy conditions [38,39].It had the highest densities and the longest leaves in shallow water, thanks to more light and increased productivity [40,41].This is surprising for Cyprus though, as waters above the summer thermocline reach 28.4 • C, which is the known upper thermal tolerance of this species [25].Cypriot P. oceanica beds might be an especially important genetic strain, as sea surface temperatures are rising rapidly in the region due to greenhouse gas emissions.Severe P. oceanica habitat loss (70%) is projected for 2050, with potential functional extinction by 2100, unless thermally resistant strains can be found [42].Cypriot strains could prove useful as transplants or as a source of seeds in other Mediterranean locations where the seagrass might be lost due to warming.
At 10-30 m depth a belt of P. oceanica meadow covered about 200 ha of sandy substrata along ca 10 km of coastline where the seagrass matte often exceeded 3 m in height.As these mattes accrete at about 1 m/kyr [43,44], and P. oceanica clones can live for hundreds to thousands of years [45], we estimate that this meadow is at least 3000 years old and is an important blue carbon store.
The Department of Fisheries and Marine Research aims to designate the western waters of Vasiliko Bay as an aquaculture zone, so it is worth discussing our results from around fish farms in the area.Near to the farms the meadow was covered in particulate waste from the cages, although storms helped flush away this fine sediment in winter.Organic matter sedimentation rates above 1.5 g m −2 d −1 can kill P. oceanica [46].All our western 15 m depth stations had fewer and smaller leaves per shoot, but more epiphytes.Increased epiphytes on P. oceanica indicate eutrophication [47], such as the nutrient enrichment that occurs around fish farms.Our station 2 was surrounded by three mariculture units and had significantly higher epiphytic mass, and lower leaf number and mass, than all the other eastern sampling stations (6-10) at the same depth.In the absence of any other obvious organic enrichment, we consider the epiphytic loading on the seagrass at this station to be due to the fish farms.Dispersion of farm wastes can influence a large area, affecting the isotopic content of P. oceanica leaf tissues over distances of kilometres [48], and the cumulative effect of different stressors reduces the resilience of seagrasses (i.e., the capacity to resist and recover from stress) [49].
Our results indicate that fish farm wastes increased seagrass epiphytic biomass at distances at least two-fold higher than the 400 m reported in previous studies [50,51].In Vasiliko Bay, fish farms operated for almost two decades over P. oceanica meadows, but from 2011 onwards they have been gradually relocated into deeper waters (37-70 m) to help protect the seagrass.In deeper waters shoot density and leaf area decreased as expected due to lower light levels.At its lower depth limit large expanses of dead matte and few plagiotropic shoots were found, indicating that the meadow is dying back [52], especially close to Vasiliko port and near a cement plant where the meadow did not extend as far into deep water.Fortunately, seagrass persists near several fish farms and monitoring has shown that the meadow stopped dying back once fish farms were moved away into deeper water [53].
European seagrass indicators include biochemical, physiological, morphological, structural, demographic, and community measures [16,54].The two indices applied in this study (PREI and BiPo) have high correlations with human pressures [34].We found that PREI consistently gave slightly higher EQR values and that the ecological status class of Vasiliko Bay, using P. oceanica, was overall 'good', with little difference among sampling stations and seasons studied.These indices are designed to reflect water quality across basin scales and not to detect localised pressures/degradation.For example, at the lower limits of the seagrass, near fish farms, there was extensive dead matte and fine particulate matter covering the leaves of P. oceanica over large distances, which was not reflected in the assessments.Furthermore, many coastal pressures may not be felt at 15 m depth and the lower limits.Managers should consider that local seagrass degradation is not always reflected in the scores of WFD indices.Seagrasses are important components of ecological quality assessments, but the seagrass indicators to be monitored depend on the set objectives.
Maintaining biodiversity and the functional balance of the fauna within a seagrass meadow food web prevents detrimental trophic cascades [49].A highly diverse coralligenous community was found at shallow stations away from coastal developments, but was absent from the industrialised areas.The ecosystem services of seagrass beds can vary greatly, as two seagrass beds with similar shoot densities and leaf dimensions may be strikingly different, in terms of the associated communities they sustain and the ecological services they provide.The good quality of a water body and the apparent health of an emblematic species such as P. oceanica, are not always indicative of the good structure and functioning of the whole ecosystem [55].Thus, an ecosystem-based approach that considers the overall functioning of the P. oceanica ecosystem, from primary producers to top predators is more informative to managers than single-species WFD descriptors [56].
Shoot density is a reliable indicator for monitoring the vitality of P. oceanica meadows [57,58].Despite anthropogenic pressures in Vasiliko Bay, and the fact that summer water temperatures are near the upper tolerance limits of P. oceanica, shoot densities recorded at both 5 m depth and 15 m depth meadows are amongst the highest reported in the literature from across the Mediterranean Sea at the same depths (Table 1).This may be due to the extremely low nutrient levels and clarity of the waters, allowing more light to reach the benthos, and buffering the effects of nutrient and suspended solid wastes.Monitoring to assess the ecological status of seagrass meadows can identify the causes and potential synergistic effects of stressors and can improve the integrated coastal management and conservation of keystone species.Considering the mounting pressures and coastal modification in the central and eastern portions of the bay, and the aquaculture expansion westwards, careful planning is needed to ensure that the marine ecosystems of Vasiliko Bay retain their good water quality status and maintain the ecosystem services offered not so long ago.During this study, fuel jetty construction in the eastern part of the bay caused irreversible loss of P. oceanica habitat in a dredged area and increased turbidity, sedimentation and epiphytic loading on the meadow downstream.This demonstrates an urgent need to upgrade and improve the quality of marine Environmental Impact Assessments plus the evaluation and follow-up monitoring of large projects.Economic gain and the contribution of a development to a country s Gross Domestic Productivity should not monopolise decision making [87].Environmental Impact Assessments largely underpin decisions for coastal developments and are plagued by inconsistent methods and lack of independent evaluation, leading to inadequate scientific rigour [88].The assessments are often rushed, leading to poor decisions based on limited data, and so refining the EIA processes is singled out as one of the most significant strategies for helping to reverse seagrass meadow decline, and for bolstering their resilience [89].
The seagrass meadows of Vasiliko Bay are among the densest in the world, and support iconic wildlife such as sea turtles, monk seals, and dolphins.We hope this study will increase awareness amongst managers and stakeholders, now that further construction plans are underway in Vasiliko Bay.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2077-1312/8/9/0717/s1,Table S1: Posidonia oceanica sampling stations, geographic location, time of sampling and replicate numbers, Table S2: Reference condition values used for applying the biotic indices BiPo and PREI to assess the ecological water quality status of the coastal area studied.From the seasonal adapted reference conditions, only the leaf surface area was varied between the two seasons, Table S3: Posidonia oceanica descriptors determined at the 30 sampling stations.Mean values ± SE are presented, Table S4 Funding: This research was funded by the Research Promotion Foundation of Cyprus (Republic of Cyprus and European Regional Development Fund) (grant agreement SMEs/Product/0609/74).Data were collected during the research project 'A holistic approach for the evaluation of ecological status of coastal areas: the case of Vasiliko Bay', in 2011-2014.

Figure 2 .
Figure 2. (A): Station 4 at 15 m near a fish farm, P. oceanica covered with epiphytes.(B): Station 8 at 15 m, P. oceanica leaves had only a few calcareous epibionts (photograph taken on same day as (A)).

Figure 3 .
Figure 3. (A): Station 9 at 15 m, what remained of a P. oceanica meadow following dredging.(B): Large pieces of living matte ripped off and lying sideways.

Figure 2 .
Figure 2. (A): Station 4 at 15 m near a fish farm, P. oceanica covered with epiphytes.(B): Station 8 at 15 m, P. oceanica leaves had only a few calcareous epibionts (photograph taken on same day as (A)).

Figure 2 .
Figure 2. (A): Station 4 at 15 m near a fish farm, P. oceanica covered with epiphytes.(B): Station 8 at 15 m, P. oceanica leaves had only a few calcareous epibionts (photograph taken on same day as (A)).

Figure 3 .
Figure 3. (A): Station 9 at 15 m, what remained of a P. oceanica meadow following dredging.(B): Large pieces of living matte ripped off and lying sideways.

Figure 3 .
Figure 3. (A): Station 9 at 15 m, what remained of a P. oceanica meadow following dredging.(B): Large pieces of living matte ripped off and lying sideways.

Figure 4 .
Figure 4. Posidonia oceanica descriptors at 5 m and 15 m depth sampling stations.The mean number of leaves and foliar surface were calculated from 20 and 40 shoots for the 5 m and 15 m depth stations, respectively.The mean shoot density was calculated from 20 replicate measurements in 0.16 m 2 quadrats.The error bars represent the SE.Stations within each panel that do not share a letter are significantly different at p < 0.05.Pairwise comparisons in panels (B), (E), and (F) were generated by Tukey's (HSD) post-hoc test, and for panels (A), (C), and (D) by Dunn's post-hoc test with Bonferroni correction.

Figure 4 .
Figure 4. Posidonia oceanica descriptors at 5 m and 15 m depth sampling stations.The mean number of leaves and foliar surface were calculated from 20 and 40 shoots for the 5 m and 15 m depth stations, respectively.The mean shoot density was calculated from 20 replicate measurements in 0.16 m 2 quadrats.The error bars represent the SE.Stations within each panel that do not share a letter are significantly different at p < 0.05.Pairwise comparisons in panels (B), (E), and (F) were generated by Tukey's (HSD) post-hoc test, and for panels (A), (C), and (D) by Dunn's post-hoc test with Bonferroni correction.

Figure 5 .
Figure 5. Posidonia oceanica leaf area index (LAI) (m 2 of seagrass canopy in 1 m 2 of seagrass bottom) heat map, interpolated using the algorithm of the 'Nearest Neighbour', and the mean LAI values of all sampling stations.The dry mass of epiphytes and leaves of P. oceanica was estimated only for the sampling stations of intermediate 15 m depth (TableS4in Supplementary Material).Confirming the visual observations, E/L varied across the Vasiliko Bay stations, both in spring (Kruskal-Wallis, χ 2 = 125.67,df = 9, p < 0.001) and summer (Kruskal-Wallis, χ 2 = 139.40,df = 9, p < 0.001) (Figure6).There was a significant interaction between the main effects (station location and season) (2-way ANOVA, F = 11.87,df = 1, p < 0.05).Several stations had higher E/L in the western side nearest to fish farms and in spring (overall mean epiphytic biomass 0.124 g shoot −1 ) than in the autumn (overall mean epiphytic biomass 0.060 g shoot −1 ) (Figure6).

Figure 5 .
Figure 5. Posidonia oceanica leaf area index (LAI) (m 2 of seagrass canopy in 1 m 2 of seagrass bottom) heat map, interpolated using the algorithm of the 'Nearest Neighbour', and the mean LAI values of all sampling stations.The dry mass of epiphytes and leaves of P. oceanica was estimated only for the sampling stations of intermediate 15 m depth (TableS4in Supplementary Material).Confirming the visual observations, E/L varied across the Vasiliko Bay stations, both in spring (Kruskal-Wallis, χ 2 = 125.67,df = 9, p < 0.001) and summer (Kruskal-Wallis, χ 2 = 139.40,df = 9, p < 0.001) (Figure6).There was a significant interaction between the main effects (station location and season) (2-way ANOVA, F = 11.87,df = 1, p < 0.05).Several stations had higher E/L in the western side nearest to fish farms and in spring (overall mean epiphytic biomass 0.124 g shoot −1 ) than in the autumn (overall mean epiphytic biomass 0.060 g shoot −1 ) (Figure6).

J 20 Figure 6 .
Figure 6.Seasonal mean epiphyte to leaf dry weight ratio (E/L), calculated for 20 shoots collected from each station at 15 m depth.The error bars represent the SE.Stations from each season that do not share a letter are significantly different at p < 0.05.Pairwise comparisons are generated by Dunn's post-hoc test with Bonferroni correction.

Figure 6 .
Figure 6.Seasonal mean epiphyte to leaf dry weight ratio (E/L), calculated for 20 shoots collected from each station at 15 m depth.The error bars represent the SE.Stations from each season that do not share a letter are significantly different at p < 0.05.Pairwise comparisons are generated by Dunn's post-hoc test with Bonferroni correction.

J
. Mar. Sci.Eng.2020, 8, x FOR PEER REVIEW 12 of 20 to have 'Medium' ESC with the BiPo index.The ESC assessed at the sampling stations 7 (in the central area of the bay) and 10 (east of Vasiliko bay, near Zygi) was assessed as marginally 'High' with the PREI (Figure 8).

Figure 8 .
Figure 8. Mean (of two seasons) ecological quality ratio estimated with the PREI and BiPo WFD indices, using the adapted reference conditions, at 10 stations set at about 1 km apart.The error bars represent the SE.

Figure 8 .
Figure 8. Mean (of two seasons) ecological quality ratio estimated with the PREI and BiPo WFD indices, using the adapted reference conditions, at 10 stations set at about 1 km apart.The error bars represent the SE.

:
Dry weight of P. oceanica leaves and epiphytes, and E/L.Mean values ± SE are presented.Stations varied between them (Kruskal-Wallis, χ2 = 228.86,df = 9, p < 0.05) and those that do not share a letter are significant different at p < 0.05.Pairwise comparisons are generated by Dunn's post-hoc test with Bonferroni correction, Table S5: Data collected from the sampling stations and used for the application of the WFD biotic indices BiPo and PREI.Author Contributions: D.K.-Conceptualisation, funding acquisition, data collection and lead author; P.K.-Data collection and manuscript review; I.S.-Statistical analyses and review; M.J.A.-Manuscript revision and editing; S.C.-Preparation of ArcGIS maps; A.L.-Data collection; J.M.H.-S.-Supervision, manuscript editing and final review.All authors have read and agreed to the published version of the manuscript.

Table 1 .
Posidonia oceanica shoot density at 5 ± 1 m and 15 ± 1 m depths in Vasiliko Bay (mean range of ten different sampling stations in each bathymetric zone) and different Mediterranean sites, ranked by highest mean reported.