Tides, Topography, and Seagrass Cover Controls on the Spatial Distribution of Pinna nobilis on a Coastal Lagoon Tidal Flat

In the last two decades Pinna nobilis, the largest bivalve mollusk endemic to the Mediterranean Sea, has recolonized the tidal flats of some coastal lagoons along the Italian Adriatic coast. In this study, we investigate the influence of tides on the spatial distribution and density of a P. nobilis population developing on a tidal flat of the Venice lagoon (NE Italy) by exploiting remote sensing technologies. Our results show that there is a threshold topographic elevation (about 0.5 m below mean sea level for the studied tidal flat) above which the number and duration of emersions become limiting factors of P. nobilis abundance. Above this elevation, the population density decreases sharply. Densely populated areas tend to occur in tidal flat depressions, where the duration and frequency of emersions are low. We find, however, that the population density has large spatial variability in response to other factors, such as the seagrass percentage cover. The density of the shells increases with increasing seagrass percentage cover, and the dense P. nobilis population (0.8–2.68 N/m2) grows within the Cymodocea nodosa meadow that has a high seagrass percentage cover (>85%). However, within the seagrass meadow with the highest percentage cover, the bivalve preferentially colonizes the portion closest to the main channel, a possible source of nutrients. The shift to C. nodosa—dominated meadows in the Venice lagoon occurred during the last two decades has probably facilitated the observed colonization by P. nobilis. Our findings provide management information for protection and/or restoration of P. nobilis in coastal lagoons.

Besides nutrient loading, there are several geophysical and biogeochemical variables that influence the survival of P. nobilis. Among them, favorable hydrodynamic conditions are a key factor in the colonization of lagoons, since in these environments P. nobilis rarely experience high water velocities that may prove harmful to them . In fact, the survival of the mollusk depends on its ability to anchor through the byssus filaments to sediments such as stones, sand, rhizomes of marine phanerogams, bio-detrital fragments, and roots. Strong hydrodynamic forces put the stability of the mollusk at risk (García-March et al., 2007). Therefore, shallow lagoon waters, where tidal flats and channels are sheltered by the presence of barrier islands and salt marshes, are an ideal environment for P. nobilis development.
Hydrodynamic forces appear to also influence the orientation of the shells. In open sea environments, this factor influences the size and spatial distribution of individuals, and their growth rate (García-March et al., 2007). A common shell orientation within a population may be linked to the survival of only those individuals that are exposed to lower drag forces (García-March et al., 2007). However, another possible explanation may be linked to food availability. In fact, microtopography effects on currents and resuspension may influence the orientation of shells, favoring the food intake of individuals that open their shells toward the current (Addis et al., 2009;Coppa et al., 2019).
An important factor influencing hydrodynamics is the presence of seagrass meadows in the areas where P. nobilis develops. The mollusk is often found within seagrass meadows and favors those composed by Posidonia oceanica or Cymodocea nodosa (Zavodnik et al., 1991). Seagrass leaves (a) perform the important task of attenuating the hydrodynamic stress (Hendriks et al., 2011;Manca et al., 2012), (b) provide protection from predators especially for young organisms and facilitate the settlement of the larvae (Hendriks et al., 2008). Seagrass meadows, however, trap suspended sediments and nutrients (Folkard, 2005;McGlathery et al., 2007) and prevent the resuspension of sediments from the bottom (Carniello et al., 2014;Moriarty & Boon, 1989). Thus, a dense seagrass population may interfere with the feeding needs of P. nobilis. There is, therefore, a trade-off between the advantage of being protected by the seagrass leaves and the disadvantage that the dense plant population may interfere with the filtering activity. In fact, P. nobilis has been found to mainly colonize the edges of seagrass meadows, where individuals can maximize the filtering action while maintaining the advantage of being protected by the plants (Coppa et al., 2013(Coppa et al., , 2019Koch et al., 2006;Manca, 2010).
In addition to hydrodynamics, other environmental forcings can change the equilibrium and influence the development of seagrasses and consequently of the P. nobilis populations that develop within these meadows. One example is water temperature, which plays an important role in controlling site specific seasonal seagrass growth (Lee et al., 2007) and in lagoons and estuaries may influence the optimal depth at which seagrass populations grow (Carr et al., 2010). Generally, the optimal growth temperature for temperate species ranges between 11.5 and 26 °C (Lee et al., 2007). However, different species may prefer more limited temperature ranges and have different optimal temperatures. For example, the growth of C. nodosa is stimulated by higher temperature, with an optimal temperature of 24.5 °C. In fact, this species is included among the subtropical species (Lee et al., 2007). On the contrary, Zostera marina, a typical temperate species, has its optimal growth temperature at around 15 °C (Lee et al., 2007).
Although, as we have seen, there has been ample interest in the effects of hydrodynamics on P. nobilis, and in particular the dragging forces, one aspect linked to hydrodynamics has been completely neglected: tidal oscillations. Such an aspect is of key importance in coastal tidal systems such as lagoons and estuaries. In these systems, different portions of the basins are subject to different flooding periods and frequencies, depending on the bottom topography and the characteristics of the tidal wave. During propagation within shallow coastal basins, the tidal 3 of 17 wave is modified by the interaction with the bottom: while channels provide preferential paths for tide propagation, tidal flats and salt marshes attenuate the tide (i.e., reduce its amplitude), also creating a phase time-lag (e.g., Friedrichs & Aubrey, 1988). The typical nonlinearity of such tidal modifications produces high spatial variability in the time and frequency with which tidal flats and marshes are flooded or emerged (Silvestri et al., 2018). Despite the importance of tides for tidal systems, we could not find any papers specifically investigating the interaction between bottom geomorphology, tidal oscillation, and the colonization and survival of the P. nobilis population, which leaves a number of open questions. The purpose of this study is to investigate if there are specific patterns in the distribution of P. nobilis over a tidal flat in relation to its geomorphology and the seagrass vegetation percentage cover, defined as the percentage of the ground surface covered by vegetation. In estuaries and lagoon systems, tidal flats are often emerged exposing P. nobilis individuals to possible desiccation . Emersion times may be a strong limiting factor to the development of P. nobilis populations; hence, we analyze the frequency of emersions and the total emersion time during the year and in different seasons. The correlation between the density of individuals (i.e., the number of individuals per unit area) and the seagrass percentage cover is also investigated, as well as other environmental forcings. Our results contribute to detect some of the main geomorphological and environmental factors that likely drove the population growth in the last 15 years and provide management information for protection and/or restoration methods. P. nobilis, is a protected species under the European Council Directive 92/43/EEC since 1992 and has been recently listed as Critically Endangered on the IUCN Red List of Threatened Species (Kersting et al., 2019) due to pathogens that are causing mass mortality events across the entire Mediterranean region (Carella et al., 2019(Carella et al., , 2020. Recent surveys report that during the first months of 2021 the Haplosporidium pinnae pathogen has attacked several individuals of P. nobilis within the Venice lagoon (http://www.ismar.cnr.it/eventi-e-notizie/notizie/pinna-nobilis-si-aggrava-lemergenza) causing their death. Therefore, our study, based on a survey performed in June 2020, can be considered as a pre-epidemic reference to detect how bottom topography, emersion frequencies, and duration control the spatial distribution of P. nobilis on the studied tidal flat. Our findings can contribute to study, preserve, and eventually restore the P. nobilis population within the Venice lagoon and in other lagoon environments across the Mediterranean region. In other sites, in fact, restoration efforts have already started, involving larval dispersal from sites unaffected by pathogens and recolonization through recruitment of resistant juveniles (Kersting et al., 2020). Other examples of effective restoration activities are reported in literature and involve transplantation of individuals of different sizes (García-March & Vicente, 2006;Trigos and Vicente, 2016). To maximize the effectiveness of transplantation efforts, at least for lagoon environments, the knowledge of the geomorphic controls that leaded to pre-epidemic P. nobilis spatial distributions is of key importance and may contribute to increase the likelihood of survival of transplanted individuals.

Study Site
We selected a tidal flat colonized by P. nobilis located within the Venice lagoon (NE Italy). With an area of about 550 km 2 , the Venice lagoon is the largest lagoon of the Mediterranean Sea. The lagoon is separated from the north Adriatic Sea by two barrier islands and is connected to the sea by three inlets (Figure 1) that serve four subbasins. The Venice lagoon is characterized by an average depth of about 1.5 m and a microtidal tidal regime with maximum tidal excursion of about ±70 cm around mean sea level at the inlets.
The study site is located within the central subbasin near the Malamocco inlet (Figure 1), where a population of P. nobilis composed of just a few individuals was first detected in 2006/2007 and has increased over time (Russo, 2012). Even though there are no studies that determine the precise size of the population, Russo (2012) documented densities of up to 16 individuals per square meter. The mudflat selected for this study is delimited in all directions by small to large channels. Two very small channels to the North and to the East separate the mudflat from Pellestrina island. The channel to the South is a minor channel that has a very limited boat traffic, while the large channel to the West is the major channel with considerable boat traffic as it connects the Pellestrina town to the inlet and the city of Venice.
The mudflat has an extent of about 28,470 m 2 and is partially covered by a seagrass meadow of C. nodosa and Zostera noltii (with an extremely limited presence of Z. marina along the south channel). Records and maps (Atlante della Laguna, 2021; SOLVe, 2021) report that both species were already present in 2002; however, at that time, the population was mixed and the portion of the tidal flat covered by seagrasses was smaller ( Figure  S1 in Supporting Information S1). Both the cover extent and the density of the two species have increased over time (see Table S1 in Supporting Information S1). This is consistent with the information available from the entire central and south basins of the Venice lagoon, where tidal flats are mainly covered by seagrass meadows, and we find specifically the three species C. nodosa, Z. marina, and Z. noltii. Nowadays the meadows are mainly composed by monospecific patches, even though mixed patches can be found far from the inlets ( Figure S2 in Supporting Information S1). The composition of seagrass meadows has changed over time, shifting in the period 2002-2017 from a multispecies population to an almost monospecific one, dominated by C. nodosa ( Figure  S2 in Supporting Information S1), a species that is still the most abundant across the tidal flats (Atlante della Laguna, 2021;SOLVe, 2021). Both the seagrass cover extent and density have also increased over time across the central and south basins (Table S2 in Supporting Information S1).
The expansion of C. nodosa can probably be ascribed to favorable conditions such as a constant increase in water temperature as will be further discussed below.
It is interesting to note that in the studied mudflat C. nodosa and Z. noltii cover two separate portions, with C. nodosa colonizing the half to the West and Z. noltii the half to the East ( Figure S1 in Supporting Information S1).

Unmanned Aerial Vehicle and GNSS Surveys
On the 23 June 2020 at 06:05 a.m. (GMT; 07:05 solar local time) unmanned aerial vehicle (UAV) images were collected at low tide with a DJI Zenmuse X4S camera (20 Mpixels, focal length 8.8 mm, 1-inch CMOS Sensor) mounted on a professional quadcopter DJI Matrice210v2. At the time of the flight, about half of the tidal flat was flooded while half was above the water level. The drone flight was planned to provide a side overlap of 70% and a front overlap of 85% among pictures, which is essential for the image matching algorithms used in Structure from Motion (SfM) photogrammetric technique (Carrivick et al., 2016). In order to ensure such overlap as well as a high resolution (image footprint with a Ground Sampling Distance-GSD, i.e., distance between two consecutive pixel centers measured on the ground, of 0.008 m), the flight altitude was ∼30 m. A total of 228 images were collected and the Agisoft Metashape Pro v 1.6.2 software was then used to produce an ortho-photo mosaic at 0.01-m resolution (i.e., pixels of 1 cm × 1 cm) through SfM technique. Before the acquisition of the images, 18 Ground Control Points (GCP) were selected across the emerged portion of the tidal flat, retrieving the coordinates with a Geomax Zenith40 Global Navigation Satellite System (GNSS) in Real-Time Kinematic (RTK) mode. The planimetric positional accuracy (RMSE) of the GNSS survey was 0.011 m and the vertical was 0.016 m in real-time mode. All point coordinates were referred to WGS84/UTM zone 32N (EPSG: 32632) reference system. Two-thirds of the GCP were used in the orthorectification process while the remaining one-third of the GCP were used afterward as Check Points (CP). This step was essential to georeference the SfM ortho-photo mosaic and evaluate the photogrammetric errors (James et al., 2019). The assessment of the SfM process uncertainty was evaluated on GCPs and CPs, considering the residuals of each point (i.e., the difference between the real coordinates of this point measured with the GNSS in the field and the modeled SfM values). The RMSE 3D (3D Root-Mean-Square Error; Remondino et al., 2017) of GCPs and CPs computed along the x, y, and z directions are 0.118 and 0.554 m, respectively. The mean of the residuals provides an indication of the accuracy of the registration process (0.031, 0.030, and 0.089 m for x, y, and z components) and the point cloud/ortho-photo mosaic (0.037, 0.039, and 0.481 m for x, y, and z components) when the GCP and CP residuals are considered, respectively. Instead, the standard deviation (StDev) of the CP residuals yields (0.027, 0.031, and 0.226 m for x, y, and z components) is an indication of the precision (Cucchiaro et al., 2018). We underline that the GCPs were placed exclusively on the emerged mudflat portion; hence, the elevation of the submersed area (which is the area of interest in this study) could not be determined due to the presence of water. For this reason, in order to retrieve the Digital Elevation Model (DEM) of the entire mudflat, we performed a second GNSS RTK survey on 26 February 2021, during a pronounced low tide that left the entire tidal flat exposed. During this second survey, a total of 313 measurements were recorded with a Topcon Hyper V GNSS system (planimetric positional accuracy around 0.017 m; vertical uncertainties around 0.024 m in real-time). The geographic coordinates were converted into WGS84 UTM 32N using the CONVERGO model provided by CISIS (Centro Interregionale per i Sistemi Informatici, geografici, Statistici). All the 331 points were then used to calculate the variogram and perform an ordinary kriging, as explained below.

Water Quality Parameters
Data collected in the period 2003-2020 with a multiparametric probe located 2.5 km away from the study site ( Figure 1; probe VE3: latitude = 45.32556796 N; longitude = 12.28574527 E) were analyzed (data provided by the Ministero delle Infrastrutture e dei Trasporti-Provveditorato Interregionale alle OO.PP. del Veneto-Trentino Alto Adige-Friuli Venezia Giulia). The study site and the multiparametric probe VE3 are located within the Malamocco subbasin, central lagoon (Cucco & Umgiesser, 2006;D'Alpaos, 2010;Matticchio et al., 2017;Viero & Defina, 2016). Together with the Chioggia subbasin (southern subbasin), the Malamocco subbasin is characterized by the greatest depths in the Lagoon, and thus by limited effects on tidal propagation (Carniello et al., 2009). Given the proximity to the Malamocco inlet, the amplitude, and phase differences of the tidal oscillation at the study site and at the multiparametric probe are negligible. Moreover, the two sites are within the same physically homogeneous area for hydrodynamic aspects (Solidoro et al., 2004a) and for other water quality parameters such as water temperature, salinity, chlorophyll-a concentrations, etc. (Ghezzo et al., 2011;Solidoro et al., 2004b). Observations and models described in Solidoro et al. (2004aSolidoro et al. ( , 2004b confirm that long-term observations of water quality performed by the probe VE3 can be reliably used to study the tidal flat selected for our study. Data were collected with a half-hourly temporal resolution, and daily mean, maximum, and minimum values were analyzed. Temperature, oxygen concentration, salinity, conductivity, and pH were measured by an Idronaut OCEAN SEVEN 316 multiparametric probe. Water turbidity is measured through a backscattering optical probe (Seapoint turbidity meter, operating at 880 nm), and is expressed in Formazine Turbidity Units (FTU), which is directly related to the suspended sediment concentration (Old et al., 2001). Chlorophyll concentration was measured with a Seapoint Chlorophyll Fluorometer. Finally, wind speed measurements recorded hourly at the station "Stazione Malamocco Porto" by the Venice municipality (see Figure 1; latitude = 45.339883 N; longitude = 12.295017 E; https://www.comune.venezia.it/it/archivio/1653) were used to calculate yearly averaged wind speed from 2004 to 2019.

Digital Elevation Model of the Tidal Flat, Emersion Times, and Frequencies
The 331 GNSS points were used in SAGA GIS to perform an ordinary kriging, which allowed us to infer the mudflat DEM at 1-m spatial resolution. We fitted a linear variogram model with no nugget in the 0-37-m distance range. In order to calculate the accuracy of our results, we performed a leave-one-out analysis that yielded an RMSE just under 0.05 m. Once the error StDev across the map was calculated, we clipped the borders where the StDev values were higher than 0.08 m in order to retain only bottom elevation estimates characterized by high estimation accuracy.
Water levels in the Venice lagoon are measured every 10 min by a network of 14 tide gauges managed by the Venice municipality. After data checking, the municipality makes hourly data publicly available through their website (https://www.comune.venezia.it/it/archivio/1653). Tidal levels refer to a local datum which is anchored to a benchmark located downtown, at Punta della Salute, defined as ZMPS (https://www.comune.venezia.it/it/ content/riferimenti-altimetrici). The ZMPS was referred to the Italian unique zero benchmark, which was established by measuring the sea level in the period 1937-1946 by the Italian "Istituto Geografico Militare" (IGM) in Genova (Hayford Ellipsoid oriented at Roma Monte Mario; Datum: Roma40) and then connected to the Venice ZMPS in 1968 through altimetric triangulations (Ferla et al., 2006). Specifically, at that time, there was a difference of 23.56 cm between the Genova and the Venice benchmarks. Unfortunately, the ZMPS elevation is not constant in time due to natural and human-induced subsidence. Specifically, based on several studies (Carbognin et al., 2004;Ferla et al., 2006;Tosi et al., 2013;Trincardi et al., 2016), it is found that only a natural subsidence of 0.5 mm/y affected the Venice benchmark after 1968 (i.e., no anthropogenic subsidence), for a total of 2.5 cm to be added to the 23.56 cm difference measured in 1968 (Da Lio et al., 2017;Ferla et al., 2006). Therefore, a total difference of 26.06 cm has been subtracted to the tide-gauge sea-level measurements in order to account for the difference with the national benchmark and for the subsidence experienced by the ZMPS datum.
As for the selection of the most appropriate tide gauge among all those available, the sea level at our study site is considered equivalent (i.e., same level and synchronous) to the level measured at the "Stazione Malamocco Porto" (see Figure 1; latitude = 45.339883 N; longitude = 12.295017 E). The correspondence has been checked on 2019 data using the finite-element hydrodynamic model developed by D'Alpaos and Defina (2007), finding a mean difference of 1.6 cm between the levels measured at "Stazione Malamocco Porto" and those simulated for the study site. Given such a close correspondence, we decided to use for our analyses the levels measured at "Stazione Malamocco Porto" instead of the simulated values.
Water levels recorded in 2019 were used to determine the amount of time that the surface is emerged during a year as well as the number of times the emersions occur. The calculated emersion times and frequencies were combined with the topography of the studied tidal flat to produce the corresponding maps and statistics. We used 2019 data instead of 2020 because in October-December 2020 the Mo.S.E. system became operational (Mo.S.E. stands for "Modulo Sperimentale Elettromeccanico," i.e., Experimental Electromechanical Module, the gates that prevent high tide in the lagoon during extreme surge events); hence, the measurements of the 2020 tidal levels are affected by the closures and are not representative of the submergence regime experienced by the P. nobilis population in the past. We exclusively used water levels recorded in 2019 to be consistent with the fact that the spatial distribution of P. nobilis shells and the topographic data were recorded one single time, in June 2020 and in February 2021, respectively. Therefore, considering that the topography of the tidal flat may have changed overtime, we limited our analysis to the water levels recorded in the prior year. However, for completeness, we show in Figure S3 in Supporting Information S1 that the variation of the mean tide levels calculated in the period 2005-2020 is negligible, with an annual variation of 0-0.1 m above mean sea level.
The number of emersions of the tidal flat in 2019 and the total emersion time (sum of the time of each emersion occurred in 2019) were calculated every 5 cm of depth, between −0.80 and −0.10 m. A fourth-order polynomial function was then fitted in order to assign the number of emersions and the emersion time corresponding to each pixel of the DEM.

Detection of P. nobilis Individuals and Seagrass Cover Assessment
Visual inspection of the ortho-photo mosaic showed that in the emerged portion of the tidal flat the seagrass partially or completely covers most of the P. nobilis individuals, preventing an accurate assessment of the bivalve presence. On the contrary, individuals are clearly detectable in the submerged areas (see Figure S4 in Supporting Information S1) and, since we could calculate the GSD value (i.e., GSD = 1 cm), we can estimate the size of the smallest detectable object in the image. Specifically, the Nyquist Sampling Criteria dictates that the pixel size of an image must be at least 2.3 times smaller than the object that is being resolved, and 2.8 times smaller to be safe (Pawley, 2006). A similar approach has been successfully tested in other contexts, e.g., landslide feature automatic detection (Tarolli et al., 2012) and ditches detection (Cucchiaro et al., 2021). Applying this criterion to our ortho-photo mosaic, we estimate that the smallest detectable shell in our ortho-photo had a maximum shell width of 3 pixels, i.e., 3 cm. Considering the possible uncertainties related to other factors (i.e., the optics of the camera, slight changes in the drone altitude, distortions during the orthorectification process, etc.) we prefer to use an even more conservative criterion, hence considering a pixel 5 times smaller than the smallest detectable object. Thus, estimating that the minimum detected shell width is 5 cm. According to measurements performed by Alomar et al. (2015) in eutrophic environments and to the growth models developed by García-March et al. (2020), we speculate that the method allows us to detect individuals that are at least 1 year of age and older. Even though the visual interpretation of ortho-photos is a common process widely applied in photogrammetry, we note that the presence of seagrass may interfere with the visual detection process, especially for smaller shells and during the summer season, resulting in an underestimation of the number of small shells with respect to larger shells. Moreover, we underline that the method does not allow the operator to distinguish between alive and dead individuals, so all standing shells were visually recognized and recorded.
On the ortho-photo mosaic, we first detected the line that separates the emerged from the submerged portion of the mudflat and excluded the emerged portion from our analyses. A 5 m × 5 m grid was overlaid on the orthophoto mosaic in the submerged area using QGIS software ( Figure S5 in Supporting Information S1). The cells that entirely or partially fell outside the ortho-photo mosaic were discarded from the analysis, as well as cells where the certain recognition of individuals was prevented by the presence of rocks or artifacts. A total of 484 cells were retained. For each 25 m 2 cell, we determined, through visual inspection, the number of P. nobilis individuals, and calculated the density as number N per m 2 (i.e., N/m 2 ).
Seagrass percentage cover was evaluated in two steps. An initial evaluation was performed based on the visual assessment of the seagrass distribution across of the entire mudflat, assigning each 25 m 2 cell to one of the four following cover classes: (a) high-cover, where plants fully cover the sediment, (b) low cover, where a mixture of plants and sediment was detectable, (c) sparse vegetation where sparse plants could be detected over a bare sediment background, and (d) bare sediment (see Figure S6 in Supporting Information S1). The high spatial resolution of the ortho-photo mosaic allowed a clear assignment of each pixel to one of the four classes by the operator also thanks to the homogeneity of the vegetation cover within each grid element. We subsequently assigned a quantitative range of vegetation percentage cover to each class by randomly selecting 10% of the 25 m 2 cells. A 1 m × 1 m grid was then defined within each selected cell thus obtaining 25 subcells of 1 m 2 ( Figure S5 in Supporting Information S1). Finally, we determined the proportion of gridline intersections falling on vegetated versus nonvegetated points. This provided an estimate of the percentage cover class within the selected 25 m 2 cells. Specifically, we found that the high-cover class corresponds to a percentage cover range of 85-100% (i.e., 85-100% of the gridline intersections falling on vegetated pixels); the low cover class corresponds to a percentage cover in the 40-85% range; the sparse vegetation class corresponds to a 10-40% percentage cover; bare sediment corresponds to a percentage cover <10%. The correlation between the P. nobilis density and the seagrass percentage cover was performed using the mean percentage cover of each class (i.e., 0.93 = high-cover; 0.63 = low cover; 0.25 = sparse vegetation; 0.05 = bare sediment).
The correlation between the P. nobilis density and the mudflat elevation was performed based on the average elevation values calculated for the 25 m 2 cells falling within the DEM. A multiple regression analysis was also performed considering both the mudflat average elevation of the 25 m 2 cells and the seagrass percentage cover determined for the cells as independent variables. The number of emersions/yr and the emersion time (hr/yr) were calculated using the fourth-order polynomial functions retrieved as explained above and then correlated to the number of P. nobilis individuals. Finally, we calculated the probability density function (PDF) for different classes of emersions.

Results
The distribution of P. nobilis shells across the studied tidal flat is shown in Figure 2a. The blue line separates submerged area from the emerged area, which are the west and the east portions, respectively. The very high spatial resolution of the ortho-photo mosaic allows the visual detection of the shells as well as the characterization of the distribution of the seagrass plants (see Figures S4 and S6 in Supporting Information S1). In the submersed portion of the tidal flat, we detected 5,252 standing shells. The population has a variable density of individuals (Figure 2a): the central area is colonized by a very dense population, with up to 2.68 N/m 2 , the shell density decreases moving away from the central part, with areas near the tidal flat boundary having densities of 0-0.2 N/m 2 . The central portion of the mudflat appears to be also the favorite area for C. nodosa development (Figure 2b). The percentage cover of the meadow, in fact, is high in the central and in the south portions of the mudflat (along the main channel), while all other areas have low cover seagrass meadow or sparse plants. The two maps in Figure 2 suggest that the presence and density of P. nobilis is not independent of the percentage cover of C. nodosa even though the high density P. nobilis portion is smaller and concentrated in a band-shaped area that goes from NW to SE. Figure 3 also shows that the dense P. nobilis population grows within the high-cover C. nodosa meadow, and the density of the shells increases with increasing seagrass percentage cover. The ANOVA test gives a significant linear correlation (P-value < 0.01). It is interesting to note that in Figure 3 the box-plot corresponding to the high-cover seagrass meadow exhibits the highest difference between the third and the first quartiles, indicating that there is a high variability in the number of P. nobilis individuals that colonize this vegetation class. This is even more evident in the maps of Figure 2: even though the dense P. nobilis population (i.e., >0.8 N/m 2 ) tends to colonize the C. nodosa high-cover class, its extent is limited to the south-west portion, and creates an elongated band-shaped structure along the edge of the meadow, while the rest of the high-cover seagrass meadow is colonized by a lower density P. nobilis population, with 0.2-0.8 N/m 2 .
The distribution of both P. nobilis and C. nodosa suggests the possible combined influence of tidal and geomorphic factors. The mudflat, may emerge during very low tides, exposing the biota that grows on its surface to possible desiccation. In order to understand the importance of this factor for  Surfaces that lie below z = −0.8 m never emerge, not even during extremely low tides, while the number of emersions notably increases for z > −0.8 m. However, this increase is nonlinear and is small for surfaces between −0.8 m < z < −0.5 m, while the number of emersions increases more rapidly moving from −0.5 m to higher elevations. A similar behavior is found for the number of hours in which the mudflat was emerged in 2019 ( Figure 4b) and for the seasonal variation of the number of emersions (see Table S4 in Supporting Information S1), especially for the summer months when the mudflat tends to be less frequently emerged than in winter months (i.e., in January and February the Venice lagoon typically experiences the major low tides).
These results were used to compute the number of emersions and the total emersion time experienced spatially by the mudflat in 2019. In order to do this, we considered the DEM computed from kriging the 331 GNSS points collected across the tidal flat and cut based on the area surveyed with the UAV system ( Figure 5a). The elevation of the tidal flat varies between z = −0.86 m and z = −0.21 m, with a depression in the central area of the permanently submersed zone. Comparing Figure 5a with the distribution of P. nobilis density classes (Figure 2a), we notice that shell density tends to be high within the depression, while areas with higher elevations are characterized by lower shell densities. A similar behavior can be noticed when comparing the P. nobilis density classes with the map of the number of emersions/year ( Figure 5b). Also, in this case, the dense P. nobilis population tends to develop where the number of emersions is low. Finally, when comparing the distribution of the bivalve with the map of total emersion time of surfaces in 2019 (Figure 5c), we notice the same behavior. Figure 6a shows the relation between the density of P. nobilis and the mudflat elevation. There is a significant linear correlation between the bivalve presence and the mudflat elevation (R 2 = 0.50, P-value < 0.01). When both the mudflat elevation and the seagrass density are considered as independent variables, the results of the multiple regression analysis show that R 2 increases from 0.50 to 0.52 (P-values < 0.01), suggesting that there must be other factors that influence the distribution and density of the P. nobilis population.
The distribution of P. nobilis density as a function of the number of emersions (Figure 6b) or of emersion time ( Figure 6c) displays a high variability in the number of individuals that develop on the mudflat portions that experience 0-40 emersions/yr, which corresponds to the emersion time range 0-100 hr/yr. This means that in this range we find sites with both dense P. nobilis populations and low-density populations. On the contrary, if we examine areas of the tidal flat with z > −0.52 m that experience a higher number of emersions (>60 emersions/ yr) and emersion times (>142 hr/yr) we never find >0.5 N/m 2 (Figures 6b and 6c). The empirical PDFs of shell density (Figure 7) computed for different intervals of topographic elevation time clarify this behavior. Mudflat areas that experience <20 emersions/yr (where z < −0.62 m) host a large variety of possible shell density values, indicating that emersion stress is not a limiting factor in this range of elevation. Areas at higher elevations, thus experiencing higher emersion times and frequency, are associated with progressively narrower ranges of possible shell densities, indicating that where there are >60 emersions/yr (areas with z > −0.52 m) the emersion stress becomes the main limiting factor for P. nobilis development and survival. Among the water quality parameters, water temperature records show a small but constant increase since 2003 for the daily mean values, as shown in Figure 8. The analysis of the other water quality parameters do not show particular trends for the period 2003-2020 ( Figure S7 in Supporting Information S1), apart from a small decrease of suspended sediment concentrations.

Discussion
Coastal lagoons and estuaries across the Mediterranean host several P. nobilis populations and so far are the ecosystems where the species has been less attacked by harmful pathogens . The hydrodynamic conditions typical of Mediterranean microtidal systems are ideal for the development of young P. nobilis individuals, which within the quiet lagoon waters can grow protected from strong currents and extreme events (García-March et al., 2007. These tidal systems provide nutrients and act as nurseries (Marrocco et al., 2018). However, a number of environmental forcings typical of tidal systems appear to limit the spread of P. nobilis population and the growth of individuals. Among the limiting factors, the influence of the tide combined with the bottom topography is one of the most important because it regulates the frequency of emersion periods and the duration of exposure to desiccation. Our analyses show that the density of P. nobilis population is highly correlated to mudflat elevation (Figure 6a). This result suggests that there is a range of elevations at which P. nobilis preferably grows and there is a higher elevation limit to the colonization of P. nobilis, which is likely related to the number of emersions and the exposure time that the species can tolerate. In fact, the number of individuals decreases rapidly as the number of emersions/year (Figure 6b) and the duration of emersion periods (Figure 6c) increase. This is likely because at higher elevations the mollusk faces the risk of desiccation. Seasonality may also have an effect in terms of desiccation. This study did not explore the consequences of emersions on the physiology of individuals; however, we notice that, given a specific elevation of the tidal flat, during spring-summer months (April to September) the number of emersions is lower than in autumn/winter months (see Table S4 in Supporting Information S1). This is because at our study site major spring-tides occur between December and February. However, we suggest that air temperature and incident sunlight may be important factors in the desiccation process of the mollusk. Thus, prolonged emersions in winter may be less dangerous to P. nobilis individuals than in spring and summer, when higher irradiance may put at risk the survival of individuals emerged for long periods. Another important factor that may influence desiccation is the presence of seagrass. We have seen that during low tides, seagrass leaves tend to cover the mollusk preventing its identification through UAV images. This circumstance may also create a sort of favorable microclimate for the shells that are protected from the direct incidence of sunlight by being covered by leaves and are less subject to continuous shell fractures and reconstructions (Prado et al., 2021). Looking at the distribution of P. nobilis individuals over the mudflat we see that areas experiencing <20 emersions in a year host the majority of the population, but interestingly in these areas we can find a wide range of population densities (Figure 7). Therefore, in these areas, the number of emersions and the emersion time are not the only factors driving the distribution of P. nobilis. On the contrary, there must be also other factors affecting the presence and distribution of the mollusk population, as, e.g., the food availability. As we move toward shallower areas, where the emersions increase, the number of individuals quickly decreases, and we never find >0.5 N/m 2 in areas that experience >60 emersions in 1 year (Figure 7). The desiccation conditions related to frequent emersions are clearly the driving factor here, limiting the colonization of the bivalve.
The distribution of P. nobilis over the mudflat confirms our hypothesis: the large majority of the population tends to grow in the depression in the central part of the mudflat that was flooded during the UAV survey (Figure 5a), where the frequency of emersions and the total emersion time are very low (Figures 5b and 5c). Few individuals survive at higher elevations, which are characterized by complete absence of individuals over large surfaces and are often emerged (Figure 5b) for long periods (Figure 5c). Our results also suggest that P. nobilis individuals tend to survive at depth lower than −0.5 m. This appears also as the threshold depth at which we observe a clear change in the number of emersions (Figure 4a): the increase in the number of emersions of the mudflat with elevations between −0.8 m and −0.5 m is much lower than the increase of emersions of the surfaces higher than −0.5 m. A similar relation is evident when looking at the amount of emersion time of the mudflat with respect to increasing elevation (Figure 4b). The elevation −0.5 m corresponds to 63 emersions, for a total of 171 hr in a year ( Figures 4a and 4b); hence, it is in good agreement with the results obtained by calculating the PDF (Figure 7).
Our results confirm that, as in open water ecosystems, in the Venice lagoon P. nobilis colonizes preferably dense C. nodosa meadows (Figure 2), and there is a significant positive correlation between the density of shells and the seagrass percentage cover (Figure 3). We speculate that, given the absence of harmful hydrodynamic conditions in lagoon environments, one of the main reasons that lead P. nobilis to colonize seagrass meadows in tidal systems may be related to protection from desiccation during low tides and emersion periods. What is evident from Figure 2 is that, within the dense seagrass meadow, P. nobilis prefers to grow along the edge toward the main channel to the west, which is directly connected to the lagoon inlet, where nutrients are likely to be more abundant and easier to trap since the tide rises from that direction. Other studies have reported a similar pattern (Coppa et al., 2013(Coppa et al., , 2019Koch et al., 2006;Manca, 2010), even though in open water environments, confirming the trade-off strategy between living protected by the seagrass leaves and the maximization of nutrient trapping. Even though P. nobilis colonizes preferably dense seagrass meadows, this does not mean that all dense seagrass meadows host P. nobilis. This is evident looking at the results of the multivariate analysis: the small increment of R 2 shows that the seagrass cover gives a small contribution in explaining the density of P. nobilis and there must be also other factors affecting the density of the bivalve, as discussed above. As previously discussed, starting in 2006/2007, P. nobilis has recolonized the Venice lagoon and such behavior might indicate that some environmental factors have changed during those years, favoring the spread of this bivalve. As we have shown in Figures 2 and 3, one of the main factors that created a favorable environment for P. nobilis has probably been the expansion of the C. nodosa meadow. In the last two decades, C. nodosa has substituted other seagrass species in the Venice lagoon, as confirmed by maps and records provided by local agencies ( Figure S2 and Table S2 in Supporting Information S1). In general, the area covered by the seagrass meadows and the density of plants have increased at the lagoon scale ( Figure S2 in Supporting Information S1) as well as at the study site scale ( Figure S1 in Supporting Information S1). This behavior of C. nodosa is in line with other studies that describe it as a fast colonizer species (Cancemi et al., 2002). Being a subtropical species (Lee et al., 2007), the spread of C. nodosa may have been facilitated by a small but constant increase of water temperature that has been recorded by the local probe VE3 through continuous measurements performed in the last 18 years (Figure 8). C. nodosa has not only colonized greater surface areas but also increased the percentage cover, developing a dense meadow that now covers the large majority of the tidal flats of the southern lagoon. There has been a decrease of turbidity recorded by the local probe ( Figure S7g in Supporting Information S1), especially during high resuspension events. Wind speed records do not show any specific decrease in the last 18 years ( Figure S7h in Supporting Information S1), sustaining the hypothesis that another factor must have influenced the lower maxima turbidity values recorded over time. Seagrass is known for its ability of trapping sediments and decreasing resuspension (Carniello et al., 2014;Carr et al., 2010;Folkard, 2005;McGlathery et al., 2007;Moriarty & Boon, 1989;Sfriso et al., 2005) and thus lower turbidity values confirm the presence of larger and denser C. nodosa meadows, as reported by maps and records (Atlante della Laguna, 2021;SOLVe, 2021). We conclude by noting that even though the constant water temperature increase may be a favorable condition for C. nodosa, it may become a threat for P. nobilis and especially for juveniles, by increasing mortality rates (Basso et al., 2015). This is not just a threat for P. nobilis: high summer temperatures and low food availability will probably characterize the Mediterranean Sea in the near future, bringing more frequent mass mortality events of invertebrates in general (Coma et al., 2009). Beside water temperature variations, other consequences of climate change and increasing human pressure put P. nobilis populations at risk since the mollusk has a strong response to geophysical and biogeochemical changes. These include  but are not limited to sea-level rise, sedimentation and erosion, pollution, invasive species, and many other factors that will particularly impact estuaries and lagoons. Even though, apart from water temperature and turbidity, no specific trends of other water quality parameters were observed in the 2003-2020 records of the multiparametric probe VE3 (see Figure S6 in Supporting Information S1), in future scenarios this could change. In fact, the data and the observations presented in this study have been collected before the Mo.S.E. system became operational; therefore, the impact of the gates has not been assessed. The closures of the Mo.S.E. system may alter the emersion frequencies and duration of the studied area. However, the gates are operated only a few times every year to prevent high tides so the influence on low tide events will be limited. There is another interference of the Mo.S.E. system which may influence tidal flat ecosystems: a recent study performed by Tognin et al. (2021) shows that the Mo.S.E. system is likely to interfere with sedimentation processes on the Venice lagoon salt marshes. Reducing the tide level during storms, the gates closure prevents the flooding of marshes with sediment-rich waters, hence decreasing the sediment contribution vital to the marsh vertical accretion. The main sources of such sediments are tidal flats, from which wind waves resuspend large quantities of sediments during storms (Green & Coco, 2014), therefore we speculate that closures of the Mo.S.E. system gates may prevent the export of suspended sediments to the North Adriatic Sea, favoring the resedimentation within the lagoon and therefore also on tidal flats. Specific studies are needed in order to understand if such processes may alter the sedimentation/erosion processes of the studied tidal flat and therefore the fate of the P. nobilis population. Such effects may however be minor because the studied mudflat is located behind the Pellestrina barrier island, which protects the mudflat from the influence of the dominant winds, coming from NE and SE.
The use of a UAV system to detect P. nobils shells over a mudflat presents advantages and disadvantages with respect to field surveys. In this study, thanks to the ortho-photo mosaic produced from the UAV images, we could detect 5,252 shells, accurately locating them in a georeferenced system without provoking any perturbation of the system or damage to the shells. This result would have been impossible to obtain with a regular field survey. However, the method presents also important limits. Depending on the resolution of the ortho-photo, and in turn on the characteristics of the camera and the height of the flight, only shells larger than a certain size can be detected. The underestimation of juveniles' number and distribution may affect the study of population dynamics. However, we underline how in tidal systems the underestimation issue may, in general, also occur when performing field surveys. In fact, even when submerged by the tide, mudflats are characterized by very low water levels that make it difficult to swim or dive, especially when there is a highly dense P. nobilis population that colonizes the mudflat surface. Walking or diving in these conditions is dangerous and risks damaging the shells. Moreover, both aforementioned activities disturb the mudflat surface and suspend large amounts of sediments, thus preventing clear recognition of shells, particularly small ones. Even when the tide is low, walking between shells is difficult and produces perturbation and potentially damages the shells. In the case that the P. nobilis population develops within a dense seagrass meadow, the detection of small shells may also be impeded by leaves that cover them. Finally, if the purpose were the detection of the exact geographic location of each shell, as in our study, the use of a GNSS receiver would pose further difficulties considering the need to perform thousands of precise measurements. A second important limitation of the method based on UAV systems is the impossibility of distinguishing dead from alive/moribund P. nobilis individuals. This is a relatively minor problem when the large majority of individuals are alive and healthy, however, it may become an important issue if we consider the recent mortality events that have decimated the bivalve populations. In our case, pathogens have attacked the Venice lagoon populations several months after our survey, therefore our data describe the pre-epidemic condition of the studied P. nobilis population. Nevertheless, a reliable method for monitoring the mortality events would be highly desirable. Given the scale of the populations developing on tidal flats, as in our study site, field observations may not be the solution. A recent paper published by Cabanellas-Reboredo et al. (2019) found that shells might be affected by pathogens several weeks before showing the effects; hence, a simple observation of their state is not sufficient to determine their health condition. The uncertainty around the amount of time that the pathogens take to kill the individual does not allow a clear distinction between alive and healthy, and alive but moribund shells. Thus, posing a number of questions about the efficacy of field surveys for tracking the mortality events. Surveys performed using UAV systems may be able to provide a solution by using, e.g., thermal cameras as suggested by Franceschini et al. (2017) or other sensors/methods not yet explored.
Finally, our results are based on one single UAV flight, and no previous observations of the population distribution are available, therefore no conclusions on the temporal evolution of the system can be inferred. A comprehensive monitoring effort will be needed, at least yearly, to discriminate between dead and alive shells, monitor the variations of the water quality parameters, record the morphological variations of the bottom topography, monitor the variation of the seagrass percentage cover and assess the possible influence of the Mo.S.E. system on the P. nobilis population.

Conclusions
The limits imposed on the survival of P. nobilis individuals by the number and duration of emersions is a key factor in the colonization of tidal systems. The population density patterns analyzed in our study show that above a certain topographic level, the frequency and duration of emersions are the main limiting factors for P. nobilis development and survival. Fewer individuals are recorded above this topographic level, where the risk of desiccation is high, while the majority of the population tends to colonize the mudflat depressions, where the number and frequency of emersions are low. In these deeper zones, the population densities range from low to high, which indicates that within a favorable environment there are other factors driving the observed density variability. Among such factors, the presence and percentage cover of the seagrass meadow is important, since we show that the dense P. nobilis population, with 0.8-2.68 N/m 2 , grows within areas characterized by high-C. nodosa meadow cover (seagrass percentage cover >85%), and the density of the shells decreases with decreasing seagrass percentage cover. However, even considering both topography and seagrass cover as independent variables, we explain the majority but not the whole variability of the P. nobilis density. Therefore, there must be other factors driving the P. nobilis density patterns, probably linked to the availability of nutrients. In fact, we find that the bivalve tends to grow along the edge of the high-cover seagrass meadow that develops toward the main channel which connects to the inlet. We speculate that there is a trade-off between the advantage of being protected by the dense seagrass meadow and the interference that plants, trapping particles with their leaves and stems, may have with the P. nobilis filtering activity.
The presence and state of the seagrass meadows may also have had a role in the recolonization of P. nobilis of the Venice lagoon starting from 2006/2007. Maps and records show that the seagrass density and the extension of the meadows have increased over time, and the slight but constant increase of the water temperature during the last 18 years probably triggered this shift toward monospecific C. nodosa meadows. The presence of denser seagrass meadows has likely had a role in the decreased water turbidity that we have observed in the area, with positive feedback for the seagrasses and consequently for the colonization of P. nobilis.
The P. nobilis population of the Venice lagoon has been recently attacked by the pathogen H. pinnae as reported by the Centro Specialistico Ittico dell'Istituto Zooprofilattico Sperimentale delle Venezie (IZSVe; http://www. ismar.cnr.it/eventi-e-notizie/notizie/pinna-nobilis-si-aggrava-lemergenza). Given the current situation, our study can be considered as a pre-epidemic reference since the UAV survey was performed several months before the pathogens attacked the population. Since our research has detected some important geomorphological and environmental factors that drove the population growth over the last two decades, it provides management information for protection and/or restoration actions.
As for the use of UAV systems for monitoring P. nobilis within lagoon waters, we underline how this methodology allowed us to detect >5,000 individuals, thus building a georeferenced data set impossible to create with direct field observation methods. The method is extremely effective for the detection of standing shells, especially when a layer of water of 20-30 cm covers the seagrass meadow colonized by P. nobilis. On the contrary, when the tide is very low and the mudflat surface is emerged, the seagrass may cover the shells, especially in summer, producing errors in the estimation of their number and distribution. Therefore, very low tide conditions should be avoided to prevent the underestimation of the population. The two most important limitations of the method are the impossibility of distinguishing between alive and dead individuals and the underestimation of recruits and small individuals because, depending on the image resolution, only shells larger than a certain size are detectable. Such limits must be carefully considered when using UAV observations for modeling population dynamics.

Data Availability Statement
Data set for the article is freely available online in the Zenodo repository (DOI: https://doi.org/10.5281/ zenodo.4812987). The UAV ortho-photo mosaic is freely available online in the Zenodo repository (DOI: https:// doi.org/10.5281/zenodo.4812727). GIS databases and shape files on the distribution and abundance of seagrass species are freely available online in: (a) Atlante della Laguna (2021) del Veneto-Trentino Alto Adige-Friuli Venezia Giulia. We acknowledge the support of the CIRGEO center (University of Padova) and specifically we thank Antonio Vettore for having provided the Topcon Hiper V GNSS system and Alberto Guarnieri for his precious help with the GNSS measurements. We thank Tegan Blount for the English language review and her useful suggestions. We also thank five anonymous reviewers and the associate editor who contributed key suggestions that allowed us to improve the quality and significance of the paper.