Assessment of the ecological structure of Posidonia oceanica (L.) Delile on the northern coast of Lazio, Italy (central Tyrrhenian, Mediterranean)

The ecological structure of Posidonia oceanica (L.) Delile meadows was evaluated on the northern coast of Lazio, Italy (central Tyrrhenian, Mediterranean sea). This is an infra-littoral zone with a wide range of anthropogenic activities and high geo-morphological variability, which reflects heterogeneity in shoot density, leaf morphology and biomass in fragmented patches. Genetic variability in populations corresponds to the formation of 3 sub-clusters, in the diverse impacted zones (north, centre and south), being correlated to the geographical distance between sites. AMOVA estimated a high genetic variation showing 43.05% individual differences within populations with a marked differentiation among the populations (56.9%) indicated by Fst value (0.57). These results revealed the role of the genetic structure of seagrasses for determining selectivity of fragmented habitat, in response to natural drivers. They showed that site-specific self-recruitment is related to biodiversity capacity and to the geo-morphological characteristic of the coast.


Introduction
Coastal areas are characterized by environmental disturbances due to both natural process and anthropogenic activities, with consequent impacts on marine ecosystems (Benoit and Comeau 2005). Posidonia oceanica (L.) Delile meadows are considered as indicators of these impacts (Pergent et al. 1995;Pergent-Martini et al. 2005;Martinez-Grego et al. 2008;Lopez y Royo et al. 2011;Bennett et al. 2011;Personnic et al. 2014), which may influence their natural structural variability, causing a general decline of biomass, density and coverage (Duarte 2002;Orth et al. 2006;Boudouresque et al. 2009;Marbà et al. 2014). In the Mediterranean Sea, P. oceanica forms widespread ecosystems with high levels of biodiversity and productivity (Boudouresque et al. 1984;2006; and, like other seagrasses, it is important from the ecological, geological and economic point of view (Costanza et al. 1997;Spalding et al. 2003;Vassallo et al. 2013). Moreover, this plant colonizes different types of substrates (e.g. sand, rocky, etc.) and generates meadows that occur in a wide range of morphotypes, ranging from continuous beds on the seafloor to fragmented patches of different shape and size (Molinier and Picard 1952;Den Hartog 1970;Semroud 1993;Boudouresque et al. 2012;Bonhomme et al. 2015). The monitoring of structural and functional descriptors of P. oceanica meadow health status (Buia et al. 2004) provides important information about the vitality and dynamic of meadows (Pergent-Martini et al. 2005), as well as the human influence on the environment (Montefalcone 2009). Recently, a new suite of descriptors has been proposed to assess the ecological status of P. oceanica meadows (Rotini et al. 2013), taking into account the physiological features of the plants by means of biochemical (Arnold et al. 2012) and genetic markers (Micheli et al. 2005(Micheli et al. , 2015Procaccini et al. 2007;Macreadie et al. 2014) as standard indicators.
The objective of this work was to assess the ecological structure and genetic patterns of P. oceanica meadows along the northern coast of Lazio to contribute to the conservation of seagrass by future management activity. For this aim, we used functional descriptors (phenotypic variables, including biomass and density) and RAPD (Random Amplified Polymorphic DNA) molecular markers. The PCR-based RAPD marker technique was recently used in seagrasses for estimating the influence of environmental disturbances on phenotypic variables and genetic diversity of populations (Kim et al. 2019), and for distinguishing genetic variation of populations at individual level (Dilipan et al. 2017). At multiple spatial scales, RAPD analyses was conducted to assess the genetic variability of seagrass, evidencing genetic differences closely correlated with coastal currents (Jones et al. 2008). In previous studies, RAPD allowed us to examine the genetic structure in meadows located in the Ligurian Gulf (North Tyrrhenian coast) that showed a decrease in genetic diversity along an anthropogenic disturbance gradient at small scale (Micheli et al. 2005), revealing the role of currents (Micheli et al. 2010b;Rotini et al. 2013). Furthermore, we were also able to evaluate changes in genetic structure of the meadows and their resilience over a decade (Micheli et al. 2012).
Here we used RAPD to detect the genetic structure of P. oceanica meadows in the fragmented habitat of the northern coast of Lazio, in order to compare it to that of other Tyrrhenian populations previously studied. Such information was detected in an area where natural drivers occur at exceptional conditions and where no similar works have been carried out before.

Study area
The study area between Marina di Tarquinia and Santa Severa (Lazio, Italy, Mediterranean Sea; Fig. 1) is located within two physiographic units, one extending from Monte Argentario to Capo Linaro, the other from Capo Linaro to Capo d' Anzio. Here the appearance of the coast is intimately related to the morphology of the seabed as the isobaths show a very uneven underwater sea-floor. In fact, according to Anselmi et al. (1978), different coastal geo-morphological types (coastal morpho-type) are present in the study area.
The northern part hosts the Mignone river floodplain, which is characterized by small sandy beaches and a rocky coastal terrace. From Civitavecchia to Santa Marinella, the coastline is dominated by the "Tolfa Mountains", which form a promontory characterizing the coastal morpho-type (terraces coasts). This promontory separates the southern physiographic unit from the northern one, and is crossed by small streams (e.g. Marangone stream) with local continental contributions. The southern part presents a small portion of coastal terraces and beaches with a prevalently sandy coast.
This coastal area is characterized by the presence of the littoral currents having a prevailing direction from south to north following a coastal local dynamic connected to high geo-morphological variability of the sea bottom, which can generate barriers to gene flow. In fact, in the study area the prevailing wind events come from to the southeast (data provided from the weather station of C-CEMS), inducing a sea current direction to north (Bonamano et al. 2015;. However, P. oceanica flowers appeared in the meadows, and seeds were found in the central zone, near the Marangone stream, during November 2013, as was observed in other sites of the area studied (Gnisci 2014;Cognetti de Martiis 2016).
The area is also characterized by a very large port (Civitavecchia harbour), two important power plants located in the northern part of Civitavecchia, and a dense urban environment constituted by the municipalities of Civitavecchia and Santa Marinella, which altogether form a single urban aggregate (Fig. 1).

Field and sampling work
Field work activities were carried out during the late spring 2013, along 40 km of coastline. By SCUBA diving, shoots were collected into 18 sampling areas (6 shoots per sites, 3 sites per each station) from May 3 rd to June 19 th 2013 ( Fig. 1, Table 1). For field work, we considered sufficient this temporal range as meteorological conditions  were constant in the study area, and because P. oceanica meadows undergo changes in a larger time scales. In each station, three sites were selected to estimate density and coverage of the meadow. By hierarchical sampling design, shoot density was calculated as shoot number/m² counting the number of shoots in nine random squares (40 × 40 cm) (Buia et al. 2004;Panayotidis et al. 1981).
On the seafloor, the percentage in coverage of P. oceanica was estimated counting the number of sub-squares occupied by the plants in three gridded squares of 1 × 1 m (Boudouresque et al. 2007).
At the same time, eighteen shoots were randomly collected in each station and stored in sea-water at 4 °C, for laboratory analyses.

Morphological measures and biomass
In the laboratory, leaves of the shoots were washed in distilled water and biometric variables, such as number, length and width of juvenile, intermediate and adult leaves per shoot were measured in each station, according to Giraud's classification (Giraud 1977). Leaf Area Index (LAI) average was calculated according to Buia et al. (2004). Furthermore, the leaves were cleaned (by a razor blade to remove epiphytes) and dried at 80 °C for 72 hours, to obtain the dried weight, i.e. the dry biomass (g dm/shoot).

RAPD genetic analysis
Genetic analyses of P. oceanica populations were performed on the same shoots collected in each of the 18 sampled stations. According to Micheli et al. (2005), the young leaves were washed in distilled water, frozen in liquid nitrogen at -180 °C and stored at -80 °C before the extraction of genomic DNA, as reported in Dellaporta et al. (1983). Subsequently, DNA was amplified by TAQ Polymerase (Applied Biosystems), using PCR (Polymerase Chain Reaction) conditions similar to those described by Echt et al. (1992) with some modification involving reaction buffer (100 mM Tris-HCL, 500 mM KCL and 25 mM MgCl2). Amplification reactions were carried out in a thermal cycler (Perkin-Elmer) with the following temperature program: 5 min denaturation cycle at 94 °C and, subsequently, 40 cycles of 30 s at 37 °C annealing temperature and 6 min at 72 °C synthesis temperature. For each DNA (1 μg/μL) amplification, the AmpliTaq DNA Polymerase (0.4 μg/μL) and 10 different easily repeated oligonucleotides (10 mM) were used for their capacity to discriminate bands and scoring them as present/absent. Amplification products were separated by agarose gel electrophoresis (1.4%), stained with ethidium bromide, visualized on UV light and photographed. The different DNA fingerprints obtained were analyzed for the presence (1) or absence (0) of the bands, in order to determine the percentage of polymorphism calculated as the number of polymorphic bands out of the total number of bands.

Statistical analyses
Morphological and structural variables of P. oceanica meadows were statistically analyzed by means of PAST (Hammer et al. 2001) and R Development Core Team (2008) software. Three way Analysis of Variance (ANOVA) was employed to test the effect of three main factors, i.e. location, depth and substrate typology, on meadow structure parameters and leaf morphology. Furthermore the non-parametric Kruskal-Wallis test was performed to check the level of significance (p < 0.05, p < 0.01 or p < 0.001) of the main factors (location, depth and substrate typology) effect on number of juvenile, intermediate and adult leaves in each shoot, which data no present a normal distribution. The Mann-Whitney post-hoc pairwise comparison was performed to valuate significant differences (p < 0.05, p < 0.01, p < 0.001) among different sampling area (North, Centre and South). Intra-and inter-population genetic distance has been determined using the NTSYS 2.0 program (Rohlf 1993) generating a similarity matrix (Nei and Li 1979). Cluster analysis was performed on the similarity matrix with the SAHN program using UPGMA (Un-weighted Pair Group Method with Arithmetic Mean) generating the dendrogram. A cophenetic correlation has been used to measure goodness of fit for the cluster analysis (r) by comparing it with the original pairwise distance matrix (Euclidean squared distance) using a Mantel test (MXCOMP program).
Then, Principal Coordinate Analysis (PCoA, NT-SYS software; Rohlf 1993) was performed to elucidate the distribution of the samples by deriving genetic distances.
Subsequently AMOVA (Analysis of Molecular Variance) was performed with AR-LEQUIN 3.5 program (Excoffier and Lischer 2010; http://cmpg.unibe.ch/software/arle-quin3) to elucidate the genetic variation and relationships between individuals within and among populations. Using AMOVA, we estimated the number of genetically homogeneous groups of individuals based on their genotypes, at multiple loci, using the Bayesian approach. The method assumes a model in which there are K populations, each characterized by a set of allele frequencies at each locus. For each individual, the software computes the proportion of the genome that can be assigned to the inferred K populations with a high probability. Because RAPD markers are dominant, we assumed that each band represented the phenotype at a single biallelic locus (Williams et al. 1990). The heterogeneity of band frequencies across all populations was tested using the Fixation index (Fst) in ARLEQUIN software. Fst is a measure of population differentiation (Meirmans 2006), and represents the most frequently used method to analyze population structure concerning gene flow.

Leaf morphology, biomass, density and coverage of the meadows
Results of leaf morphology (number of leaves/shoot, width and length), Leaf Area Index (LAI), biomass, shoot density and coverage of the plants collected in the 18 sampling stations, are reported in Table 2. The presence of highly fragmented meadows (evidenced in Suppl. material 1: Fig.  S1), in a large range of depths (from 4.9 m to 13.4 m. Table 1), was confirmed by the high variability of coverage found among stations, ranging from 6 to 98 % (St 16 and St 14 respectively), and by the coefficient of variation (72.4%) ( Table 2).
Three-way ANOVA (Suppl. material 3: Table S1), was tested to detect significant factors in a multifactor model, in which there is one dependent variable (shoot density of P. oceanica) and three independent variables (depth, location and substrate typology).
The three way ANOVA showed differences in sea-bottom coverage (P < 0.05) due to different substrate types (Table 3). In all areas, the mean shoot density estimated was 276 ± 133 SD (number of shoot/m² ± Standard Deviation), ranging from 141 ± 62 SD to 653 ± 145 SD (St 8 and St 5, respectively). Biomass mean value 0.70 ± 0.213 SD (g dm/shoot) showed the highest values in the stations at intermediate depths, ranging from 1.03 ± 0.49 SD to 0.40 ± 0.11 SD.
All the shoots showed similar values in leaf number (5.92 ± 1.03 SD mean value) and width (0.9 ± 0.05 SD mean value).
Due to the differences in depth and water transparency of the 18 sites, we found a highly variable leaf length with values ranging from 24.24 ± 14.22 SD to 49.44 ± 17.11 SD (mean values of 35.05 ± 7.3 SD), showing statistically significant differences in the three populations (ANOVA, P < 0.05) ( Table 2, Table 3, Suppl. material 3: Table S1).
In the Table 3, the multi-factors three-way ANOVA performed on shoot density of P. oceanica, showed no significant difference among the sites located in the north, central and south areas (ANOVA, p = 0.1868); on the contrary significant differences were found among stations located both on different depths (ANOVA, p < 0.001) and on different substrates (ANOVA, p < 0.05).
The non-parametric Kruskal-Wallis test was performed to verify the presence of significant differences in number of juvenile and intermediate leaves among the three sampling location (North, Centre and South). The test highlighted significant differences in the number of juvenile (p < 0.001) and intermediate (p < 0.001) leaves. Following Kruskal-Wallis results, the Mann-Whitney post-hoc pairwise comparison test was performed; it highlighted the lowest values of the juvenile leaf number, and the highest values of the intermediate leaf number in the central area compared to the other two locations. The adult leaves number showed no significant differences among the different locations (Kruskal-Wallis, p = 0.0525) (Suppl. material 2: Fig. S2).

Genetic variability of the population
By UPGMA dendrogram (Fig. 2), based on SAHN clustering of the data, the similarity coefficient (Nei's index) among the P. oceanica samples was 0.61 ± 0.01 based on 140 RAPD total bands (mono-and polymorphic) with fragments ranging in size from 0.3 kb to 3.8 kb. The percentage of polymorphism obtained in the population was 90.14% (Table 4). Mantel test, which compares Nei's distance and cophenetic matrices, returned statistically significant values, with a matrix correlation of r = 0.95 (normalized Mantel Statistic Z); t = 25 (Mantel T-test) and p = 1.00000 (Probability random Z < observed Z).
PCoA analysis showed that the individuals sampled in the 3 different impacted zones fell into 3 distinct groups (north, centre and south populations) (Fig. 3). The first two axes accounted for 57.8% of variation, with axis 1 and 2 explaining 39.7 % and 18.1% of the variation, respectively. AMOVA (Suppl. material 4: Table S2) showed 43.05% individual differences within populations with a marked differentiation among the populations (56.9%) indicated by Fst value (0.57).

Discussion
In the meadows, the genetic structure of P. oceanica populations was clustered into three main groups located into three different zones (Fig. 2). The first population has been identified in the northern area, corresponding to sampling sites between Civitavecchia harbour and Marina di Tarquinia, an area mainly dominate by the presence of the Mignone river floodplain. As observed in UPGMA analysis, meadows located in this northern area are divided in two sub-populations, which are crossed by Mignone river. The second population (in the centre) corresponds to the sampling stations between Civitavecchia harbour and Capo Linaro with the presence of low cliffs and a very heterogeneous sea floor, with a prevalently rocky bottom. The third population, located in the southern part of the study area, from Capo Linaro to Santa Severa, is characterized by a prevalent sandy sea floor. By AMOVA, a high level of genetic diversity was recognized within (43.05%) and among the populations (56.95%) and was confirmed by PCoA (Fig. 3) when the variation was applied to genetic distance among individuals (57.8%).
In this highly heterogeneous coastal area, we found 90.14% of total polymorphism ranging from 66.67 % to 100 % in the individuals analyzed (Table 4). This high percentage of polymorphism demonstrates the high variability displayed by each individual, which is also evidenced by the UPGMA cluster analysis (Figure 2) showing a low similarity value (0.52) within each population.
In the UPGMA cluster analysis, the 0.52 similarity value is lower than the one found for Monterosso meadow (0.66), where natural and anthropogenic pressure were low, and lower than the one found at the Mediterranean basin scale (0.81) (Micheli et al. 2005).
Moreover, in the same sub-cluster ( Figure 2) are grouped the individuals of population with the same morphological features of the leaves linked to their genetic capacity to fit the same environment.
In the Table 2, the structural and functional descriptors have confirmed a high variability in all sampling areas, due to the different depth, which may particularly influence the morphology of leaves. According to other authors (Di Maida et al. 2013;Balestri et al. 2004;Pergent et al. 1995), the meadows that we found are fragmented, discontinuous and are characterized by a wide range of shoot density and substrate coverage which reflect a spatial heterogeneity of P. oceanica meadows, probably due to the complex nature of the seafloor (Tables 1, 2). Our results highlight the influence of natural drivers, as the geo-morphological features of the fragmented habitat, in the connectivity of seagrass populations. The heterogeneity patterns of the meadows that we found along the northern coast of Lazio, describe the ability of plants to respond to the environment. This depends on the availability of genetic variation improving plant growth, development and resilience (Hughes and Stachowicz 2004;Reusch et al. 2005). At this regard, Jahnke et al. (2015) suggested not only a greater resistance and resilience of individuals of higher genetic and genotypic diversity under disturbance, but also that the high diversity in high impacted sites reflects the mismatches with respect to pre-environmental impact conditions. especially because flowering and sexual recruitment are seldom observed. In fact, by the genetic diversity metrics study, they observed that the absence of low and medium levels of genetic variation in impacted locations were probably due to local extinctions of individuals that already exceeded their resistance capacity. Therefore, they suggested to use this valuable tool for restoration and mitigation projects, because the exceptional longevity of individuals (i.e. P. oceanica population) could create a temporal mismatch that may falsely suggest good meadow health status, while gradual deterioration of allelic diversity might go unnoticed (Jahnke et al. 2015).  According to the above studies, our genetic data have highlighted the influence of different habitat conditions on plants found in the three different areas, which parallel the high variability of all structural and functional parameters ( Table 2). The incidence of multiple environmental cues led plants to a more plastic response, differing among individuals, and the genetic distance between populations correlates to the distance between sites showing different environmental conditions (substrate, marine dynamic, nutrients, etc.) and subject to both natural drivers and different stressors.

Conclusion
The first investigation on the ecological/genetic structure of P. oceanica growing along the coasts of the northern Lazio has called the attention on the higher genetic variability of these meadows than others previously studied (Micheli et al. 2005;2010b;Rotini et al. 2011;. The observed genetic diversity and differentiation (Fst = 0.57) found among the P. oceanica populations suggest an adaptation capacity of the species to share different environments (rocky and sandy). This reflects in the formation of different local meadows with a spatial heterogeneity distribution along the studied area. UPGMA cluster of genetic distance revealed a high variability of seagrass population in the relatively close meadows as the result of a low level of gene flow among them. By using a dispersal simulation analysis, Serra et al. (2015) suggested that such pattern of genetic structure probably reflects historical events of recolonization from relict glacial areas and of past vicariance, which seems to persist because of the low connectivity among populations via marine currents. Whilst Jahnke et al. (2015) observed that the absence of low and medium levels of genetic variation in impacted locations are probably due to local extinctions of individuals that already exceeded their resistance capacity.
In the Mediterranean, P. oceanica plays a crucial role for the ecosystem productivity through favoring biodiversity, revealing a whole range of intra-specific levels of diversity from complete clonality to high variability (Ribeiro et al. 2016). In this littoral zone of the northern coast of Lazio, natural drivers like the prevailing current direction towards north (Bonamano et al. 2015; could affect seed transport favoring the colonization of new habitats whilst the incidence of multiple environmental impacts and the intra-meadow genotypic heterogeneity drive the plant to a more plastic response among individuals (Micheli et al. 2015;. Posidonia oceanica is considered a vulnerable species, largely susceptible to perturbations, and is thus included within Habitats Directive (92/43/EEC), which supports the Natura 2000 network of European protected areas. From a conservation perspective, and considering future increase in both anthropogenic and climatic pressures (Mannino 2018), the assessment and management of environmental quality require an adequate knowledge of the ecosystems that can be obtained only by multidisciplinary investigations (Cozza et al. 2019;Domina et al. 2018). Therefore further studies focused on the intra-variability of the meadows, by analyzing the allelic richness (Jahnke et al. 2015), are needed. These studies should include the changes of private (neutral) allelic in order to control the population demography vs genetic relationship stability necessary to reinforce their resilience as we have observed in other protected areas in the Mediterranean (Micheli et al. 2015;. Finally, since seagrasses are globally declining, and their loss is often due to synergies among stressors (Ceccherelli et al. 2018), the ecological implication of the results obtained here, such as the genetic features and the environmental conditions, could be further aspects to be considered in the monitoring surveys.