Sexual reproduction vs . clonal propagation in the recovery of a seagrass meadow after an extreme weather event

1 Centre of Marine Sciences (CCMAR), University of Algarve-Campus de Gambelas, 8005-139 Faro, Portugal. (DP) (Corresponding author) E-mail: diogosoarespaulo@gmail.com. ORCID-iD: https://orcid.org/0000-0001-9061-6850 (OD) E-mail: onno.diekmann@gmail.com. ORCID-iD: https://orcid.org/0000-0002-6141-3429 (AAR) E-mail: xanaramos@gmail.com. ORCID-iD: https://orcid.org/0000-0003-1099-4005 (EAS) E-mail: eserrao@ualg.pt. ORCID-iD: https://orcid.org/0000-0003-1316-658X 2 Department of Biological Sciences, University of Wisconsin Milwaukee, WI, 53201 0413, USA. (FA) E-mail: albertof@uwm.edu. ORCID-iD: https://orcid.org/0000-0003-0593-3240


INTRODUCTION
Seagrasses are declining worldwide due to human and natural causes (Waycott et al. 2009), leading to loss of their important ecosystem functions (Hemminga and Duarte 2000). Disturbance and stress are well known environmental factors affecting successful plant establishment (Craine 2009). Catastrophic events in particular can change seagrass cover drastically (den Hartog 1987, Williams 1988, Larkum and den Hartog 1989 and may trigger energy allocation into sexual reproduction (e.g. Reusch 2006, Alexandre et al. 2005, Hammerstrom et al. 2006. Although clonal propagation is the major mode of space occupation for many seagrasses (Marbà and Duarte 1998), under some environmental conditions sexual reproduction can be very important in population persistence and in particular for recovery from disturbances (Greve et al. 2005, Bell et al. 2008, Park et al. 2011. Several studies have focused on the relationship between genotypic diversity and resilience and resistance to perturbations and extreme stress events in seagrasses (e.g. Hughes and Stachowicz 2004, Massa et al. 2013, Evans et al. 2017. These studies report that higher genotypic diversity can have several positive effects on the response to perturbations, including increased resistance to loss and faster recovery. A better understanding of how seagrass meadows recover from extreme events and of the relative role of sexual versus clonal propagation (i.e. genotypic diversity) may thus be an important concern in conservation and management decisions (Reusch et al. 2005, Qin et al. 2014). However, most of the evidence is based on short-term experiments or at most a year of recovery. Very little is known about the relationship between genotypic diversity and disturbances in the long term, analysed at scales encompassing several years. Therefore, this paper aims to focus on understanding the relationship between sexual versus clonal propagation and disturbance over a multi-year time scale.
Zostera marina (L.), a monoecious seagrass species, is the most widely distributed seagrass species in the northern hemisphere (den Hartog 1970), forming perennial and annual populations (Bos et al. 2007). Annual populations are characterized by high sexual allocation, where all vegetative shoots turn into flowering shoots, creating seed banks which can germinate when conditions are favourable (Keddy and Patriquin 1978, Santamarı́a-Gallegos et al. 2000, Jarvis et al. 2014. The genotypic (i.e. clonal) diversity of annual Z. marina meadows is therefore expected to be very high. Perennial populations are characterized by lower sexual allocation, and therefore lower genotypic diversity, with clonal growth being the favoured strategy for area occupation (Olesen and Sand-Jensen 1994, Olesen 1999, Kim et al. 2008. Despite these general predictions, high sexual allocation has been observed in some perennial populations (van Lent and Verschuure 1994, Meling-López and Ibarra-Obando 1999, Olesen 1999. Such unpredicted behaviour in persistent populations may be linked to disturbance events (Cabaço andSantos 2012, Qin et al. 2014) such as storms.
The aim of this study was to assess how the importance of sexual and clonal reproduction varies before, during and after recovery from natural storm disturbance in a natural Zostera marina population followed over a five-year period. We proposed to test two alternative hypotheses that resulted in contrasting genotypic diversity in the recovered meadow; we hypothesized that, in the long term, seagrass recovery may lead to either 1) clone dominance due to vegetative propagation being more effective than sexual reproduction, resulting in low genotypic diversity; or 2) high genotypic diversity due to recurrent seedling recruitment.

Location
The studied meadow is located at the tip of the Troia Peninsula (Ponta do Adoxe) at the mouth of the Sado estuary, Portugal (38°29.525'N 08°54.507'W, Fig. 1). The area is a sand point characterized by very strong tidal currents. The perennial Z. marina bed is located at the beach, starting at 1.5 m depth and extending to 3 m depth.
The meadow had been monitored and sampled for population genetic analysis before its disappearance in 2008 (Diekmann and Serrão 2012). In March 2010, after intense winter storms, we observed that all 1.12 ha were lost (Cunha et al. 2013). Z. marina recovery was monitored during the spring and summer of 2010. By 2013, recovery of the meadow was observed (continuous seagrass cover of 0.57 ha) and it was resampled for genetics.

Field work
The total area of the meadow was estimated in the spring of 2009 using GPS in tracking mode while fol- lowing its contour from the boat on a day of good visibility. A comparison between the GPS tracking data and Google Earth images from the same year confirmed that the areas were the same. In 2007, 2013 and 2015, the area was estimated using Google Earth images. Before the disturbance, shoot density had been recorded in 2009. Immediately after the winter disturbance, we recorded seedling and shoot densities in April, June and September 2010. The meadow was resampled in July 2013 (three years after the disturbance), when it was visually estimated to have attained 100% of the cover it had had before its disappearance. SCUBA divers monitored along two 50 m transects perpendicular to the shore. Every 2 m along the transects, the divers systematically counted seedling and shoot density using a 25×25 cm quadrat at the interception point, and measured leaf lengths of the longest leaf within each shoot using a ruler, in a total of 25 quadrats per transect.
We performed ANOVA and, when necessary, Tukey tests (R Core Team 2017) to test mean differences before and after disturbance for the following variables: seedling density (seedling m -2 ), shoot density (shoots m -2 ) and leaf length (cm).

Genetics
Plant leaves were collected haphazardly by SCUBA diving, maintaining a distance of 1.5 m between sampled shoots. The pre-disturbance samples were collected in May 2008 (n=83) and the post-recovery samples were collected in May 2013 (n=31). The number was not equal between pre-and post-disturbance because the sampling was opportunistic, with SCUBA divers collecting shoots haphazardly and only counting them when they arrived at the boat. Samples were conserved in silica gel desiccant.
Plant DNA was extracted following the CTAB method (Doyle and Doyle 1988). After extraction, two chloroform⁄isoamyl alcohol (24:1) extractions and one ethanol (100%) precipitation were used to purify the DNA. Samples were genotyped for eight microsatellite loci (Stam et al. 1999, Reusch 2000 and followed the protocols of Reusch et al. (2000) for multiplexing. Fluorescently labelled PCR fragments were analysed on an ABI PRISM 3130 Genetic Analyser (Applied Biosystems) at CCMAR. The STRAND software (http://www.vgl.ucdavis.edu/informatics/strand.php) was used to score raw allele sizes. The R package Msa-tAllele (Alberto 2009) was used to bin the allele sizes, and ambiguities were manually reviewed.
The genotype data were analysed to estimate the number of sampled genets (G) for comparison with the number of sampled ramets (N= sample size). The probability that identical multilocus genotypes (MLG) were produced by sexual recombination (Psex) rather than being clones was estimated using the GenClone software (Arnaud-Haond and Belkhir 2007). Genotypic diversity, the proportion of genets per population, was estimated following Dorken and Eckert (2001): R = (G-1) ⁄ (N-1). We also analysed samples from before and after disturbance together in GenClone to infer whether we could find the same clones in both populations from different years.
Population genetic parameters were calculated after the removal of clone replicates for pre-disturbance and post-disturbance samples. The observed and expected heterozygosity (H O ; H E ) and the inbreeding coefficient (G IS ) were estimated using the GenoDive software (Meirmans and van Tienderen 2004) and G IS were tested for significant difference from zero using a permutation test (Goudet 1995) to detect departure from Hardy-Weinberg equilibrium. Pairwise population differentiation (F ST ) was estimated and significance (difference from zero) was tested by permutations (AMOVA, Excoffier 1992, Michalakis andExcoffier 1996) using the same software (Meirmans and van Tienderen 2004).
Using the allele frequencies in the two samples (pre-and post-disturbance), we calculated the likelihood of an individual genotype being found in each of the samples (Paetkau et al. 1995), allowing us to assign the sample from which each post-recovery individual is most likely to come. Individual population assignment was inferred in Genodive (Meirmans and van Tienderen 2004) based on a Monte Carlo test with an alpha of 0.002 on 100 replicated datasets and 11400 resampled individuals. to 0.07 seedlings m -2 (±0.04 SE) in June 2010 and no seedlings were detected in September 2010 and July 2013 (Fig. 2). Shoot density also varied over time, with an opposite trend relative to that of seedling density. It changed from 288.48 shoots m -2 (±19.30 SE) in July 2009 (before the disturbance), to total loss of shoots in February 2010. In that same year, in April we only observed seedlings, none of which had initiated clonal growth. By June 2010, shoot density was 12.62 (±3.91 SE) and by September it had nearly duplicated to 22.19 (±4.05 SE) shoots m -2 . Finally, the population reached densities similar to the pre-disturbance values in July 2013 (314.80 shoots m -2 ±9.35 SE), representing 100% plant cover (Fig. 3).

Before and after meadow comparison
There were significant shoot density differences between all dates (Tukey HSD P<0.05) with the exception of June and September 2010 (Tukey HSD P=0.73). Leaf length also increased significantly from April 2010 to July 2013, as shown by significant differences between all dates during the monitoring period (Tukey HSD P<0.05) (Fig. 4).

Population genetic structure
Different clones were identified in the samples collected before and after disturbance. The genotypic diversity estimated decreased from R=0.89 in 2008 to R=0.60 in 2013 (Table 1). Based on allelic frequencies, the individual population assignment test tagged only a single individual from the post-disturbance population as belonging to the pre-disturbance population.
Expected and observed heterozygosity were similar in both samples (permutation test, P=1) ( Table 1). Although the inbreeding coefficient was similar between populations (permutation test, P=1) ( Table 1) the G IS of the pre-disturbance population revealed heterozygote deficiency (P=0.001), whereas the post-disturbance showed no deviation from Hardy-Weinberg equilibrium (P=0.134). The pairwise population differentiation fixation index, F ST , between samples was 0.091 and not significantly different from zero (permutation test P=0.001).

DISCUSSION
Our results show that both sexual and clonal reproduction of Z. marina played an important role, but contributed differentially, in different stages of recovery of a meadow from a severe natural disturbance event that had caused its disappearance. At an early stage, the seed bank played a crucial role in the meadow recovery, shown by seedling recruitment. It has been hypothesized (Orth et al. 2006) that the natural recover of Z. marina banks is likely to be possible due to small remnant stands. We made a considerable effort to cover the total former area covered by the seagrass meadow after its decay and found no evidence of remaining shoots. Therefore, our results demonstrate the important role of the seed bank in recovering a population in which all plants were lost.
On the long term, clonal growth increased shoot density, increasing the vegetated area. This perennial meadow behaved as an annual meadow during extreme climate conditions, demonstrating the species's capacity to recover from large scale losses, via seedling recruitment, as similarly observed in other locations (Jarvis and Moore 2010, Kim et al. 2014, Qin et al. 2016. Most of all, it shows that the contribution of sexual and clonal propagation is a population trait that is very variable (Rafajlović et al. 2017), and this plasticity plays a key role in population persistence despite catastrophic disturbance events.
The new meadow originated exclusively from sexual recombination from the previous population and not from seed dispersal from a differentiated population, as shown by the absence of genetic differentiation between the pre-and post-disturbance samples. All   other populations in the region are very differentiated and separated by >100 km (Diekmann and Serrão 2012). Closer to this study site there were small patches within a distance viable for seed dispersal but only one survived the storm. This patch was less than 10 m 2 , a very small patch inside the Sado estuary, most likely not containing sufficient plant material to have given rise to this new population. Also, the small surviving patch in the Sado River Estuary was studied in previous works and had very low genotypic richness, leading to a prediction of low sexual allocation (Diekmann and Serrão 2012).
Our study site is located in an extremely dynamic area, mainly due to strong tidal currents and storms during winter months (Paulo et al. 2019 -Supplementary materials). In such environments, high sexual allocation by seagrasses is well documented (van Lent and Verschuure 1994, Meling-López and Ibarra-Obando 1999, Santamaría-Gallegos et al. 2000. The high genotypic diversity found in the pre-disturbance population demonstrates that this meadow allocated resources into sexual reproduction in the years previous to disturbance. As a result, a seed bank allowed the population to recover. The observed reduction in genotypic diversity in 2013 shows that shoot recruitment from sexual propagules, which started as 100% seedlings (R=1), was diluted with clonal growth over time. Nevertheless, the post-disturbance population genetic diversity was still high.
The pre-disturbance population showed heterozygosity deficiency, which may be explained by inbreeding (Diekmann and Serrão 2012), suggesting that in the past, sexual reproduction occurred between closely related individuals or selfing. The post-disturbance population, sampled after three years of recovery, had no heterozygote deficiency, indicating random mating. This change in mating system might hypothetically be a consequence of the different clonal structure of the population. Although male and female flowers in the same inflorescence are not mature in synchrony, those in different flowering shoots from the same clone can reproduce with each other (Reusch 2001). In contrast, in a recent population formed entirely by smaller clones, each clone is not likely to have many more flowering shoots than any other one.
Southern edge populations, such as the one studied here, are particularly important for the persistence and evolution of Z. marina as a species. These populations are small and reproductively isolated from others (Diekmann and Serrão 2012). Although at the worldwide scale, eelgrass conservation status is of least concern, the current trend is for populations to decline (Short el al. 2010). The southern populations in particular are considered highly endangered and many have recently disappeared (Cunha et al. 2013). As an example of the species vulnerability in this region of the world, 5 km northwest from our study area, 10 ha of eelgrass meadow were recently lost and did not recover naturally (Cunha et al. 2014, Paulo et al. 2019). Due to the opportunistic nature of this research, we did not have the opportunity to study sexual and clonal allocation. As future recommendations, we suggest the study of the sexual investment in these populations to better understand the chances of recolonization after localized extinction. Seedbank size should be assessed yearly, together with the meadow flowering rate, information that is vital to understand the chances of survival in case of another mass mortality event. Such knowledge can be used to plan seed-based restoration programmes (Tanner andParham 2010, Marion andOrth 2010). Because Z. marina seed banks have a transient nature, the window of opportunity for recolonization might be short. Seed banks' highest germination success is limited to 6 months and falls to less than 5% in 15 months (Jarvis et al. 2014). It is still unknown how long it takes for a seagrass population to undertake another extreme event and still recover. Theory predicts that after colonization is completed, sexual reproduction can be favoured, decreasing clonal dominance (Rafajloví et al. 2017). Our study period was too short to confirm this hypothesis. Therefore, in future research, it is very important to quantify sexual investment (flower and seed density) and seedling success vs. clonal growth to determine long term trade-offs and vulnerability status to seagrass meadows.
This study demonstrated that the relative roles of sexual versus clonal reproduction in partially-clonal organisms, such as seagrasses, vary with the temporal scales of population disturbance and are best understood when analysed on multi-year scales. The study also revealed that sexual reproduction is of extreme importance in population maintenance and resilience in the face of unpredicted future climate regimes.