Mating system variation in the eelgrass, Zostera marina

Genetic diversity can modulate a population’s response to a changing environment and plays a critical role in its ecological function. While multiple processes act to maintain genetic diversity, sexual reproduction remains the primary driving force. The eelgrass (Zostera marina) is an important habitat forming species found in temperate coastal ecosystems across the globe. Recent increases in sea surface temperatures have resulted in shifts from a perennial to a largely annual life-history strategy at its southern edge-of-range. Given that mating systems are intimately linked to standing levels of genetic variation, understanding the scope of sexual reproduction can illuminate the processes that shape genetic diversity. To describe and characterize edge-of-range eelgrass mating systems, developing seeds on flowering Z. marina shoots were genotyped from three meadows in Topsail, North Carolina. In all meadows, levels of multiple mating were high, with shoots pollinated by an average of 8 sires (range: 3 – 16). The number of fertilized seeds (i.e., reproductive success) varied significantly across sires (range: 1 – 25) and was positively correlated with both individual heterozygosity and self-fertilization. Outcrossing rates were high (approx. 70%) and varied across spathes. No clones were detected and kinship among sampled flowering shoots was low, supporting observed patterns of reproductive output. Given the role that genetic diversity plays in enhancing resistance to and resilience from ecological disturbance, disentangling the links between life-history, sexual reproduction and genetic variation will aid in informing the management and conservation of this key foundation species. 

Methods

Study Organism  

During Zostera marina’s reproductive season, reproductive shoots (i.e., flowering shoots) extend from a basal node and form multiple flowering branches (rhipidia). Each rhipidium contains additional branches with several spathes (i.e., inflorescences) containing flowers (De Cock 1981; Thayer et al. 1984, Kuo & Den Hartog 2006). Spathes contain a spadix with rows of both male and female flowers clustered in ratios of 2:1. Flowering occurs in two stages: (1) the flowering and exposure of ovaries and (2) the opening of anthers and release of pollen, each separated by approximately 5 days (De Cock 1980). Past studies describe sequential flowering of spathes (i.e., temporal dichogamy) within Z. marina reproductive shoots beginning with basal structures (De Cock 1980, 1981). Specifically, one spathe per rhipidium flowers at a time, during which the next spathe of the same rhipidium is not yet entirely developed. The second spathe starts flowering a few days following the first spathe opening its anthers. Temporal dichogamy seemingly does not exist across rhipidia (De Cock 1980, J. Jarvis unpubl. data). In the present study, spathe position is therefore used as a proxy for reproductive timing, and rhipidium position is used as a proxy for height within the water column, consistent with previous literature (e.g., Furman et al. 2015, Follett et al. 2019). 

2.2. Study Sites & Sample Collection  

Three seagrass meadows located in Topsail Sound, North Carolina were sampled to characterize Zostera marina mating system variation. The meadows were in a shallow coastal lagoon (34.22 N, 77.37 W) protected from the Atlantic Ocean by Topsail Island. This lagoon contains Z. marina in both monospecific and mixed assemblages (i.e., co-occurring with Halodule wrightii and Ruppia maritima) (NCDEQ 2021). Sites were classified as shallow subtidal (depth < 2m MLLW) and were on a narrow shelf between the Intracoastal Waterway and the adjacent shoreline. Twelve flowering shoots were haphazardly collected roughly 5m apart from each site (n = 36 shoots total) on 4 May 2021 and transported on ice to the University of North Carolina Wilmington’s Center for Marine Science. Shoots were then cleaned and blotted dry, and the morphological position of each seed was recorded. Each seed was given a position within a spathe, assigned to a given spathe and to a given rhipidium. Labels were assigned in ascending numeric order from increasing proximity to the rhizome (e.g., basal positions were given a value of 1).  

Seeds were removed from spathes, cleaned, blotted dry, and tested for viability using the “squeeze test” by gently compressing individual seeds with a pair of tweezers (Marion & Orth 2010). Those with a seed coat that compressed were considered nonviable. Viable seeds were removed from the seed coat and genotyped. Historically, the “squeeze test” is performed on seeds that have been released from the flowering shoot; however, if this had been done, the position of the seed within and the identity of the parent plant would not have been known. Moreover, seeds were sampled in May at the peak of the eelgrass reproductive season in North Carolina where most sampled seeds were mature (i.e., at stages 4 and 5; Combs et al. 2021). Indeed, when the sites were re-visited one-week post-sampling, all seeds had been released from their spadices (Jarvis, pers. obs.) indicating that seeds collected as part of this study were mature at the time of collection. Moreover, any immature seeds released from the spadix within 1 week were effectively nonviable; they would not have successfully germinated.  

Seed Genotyping 

DNA was extracted from viable seed samples using a PowerPlant® Pro DNA Isolation Kit. Ten microsatellite loci previously characterized for Z. marina (Reusch et al. 1999, Reusch 2000a, Oetjen & Reusch 2007, Oetjen et al. 2010) were amplified in two multiplex Polymerase Chain Reactions (PCR). Individual primer working stocks contained 1μL of 10μM fluorescently-labeled forward primer and 10μL each of 50μM unlabeled forward and reverse primers diluted in 80μL of ddH2O. Primers were then combined into two primer mixes – each containing five different primers. PCR conditions for all multiplex conditions were as follows: 95.0°C for 15 minutes; 2 cycles of 94.0°C for 15 seconds, 60.0°C for 30 seconds, 72.0°C for 45 seconds; 2 cycles of 94.0°C for 15 seconds, 59.0°C for 30 seconds, 72.0°C for 45 seconds; 2 cycles of 94.0°C for 15 seconds, 58.0°C for 30 seconds, 72.0°C for 45 seconds; 2 cycles of 84.0°C for 15 seconds, 57.0°C for 30 seconds, 72.0°C for 45 seconds; 28 cycles of 94.0°C for 15 seconds, 56.0°C for 30 seconds, 72.0°C for 45 seconds; and a final 2 minute extension at 72.0°C. Following PCR, two reactions were prepared: one containing 0.5μL of each PCR product from each of the multiplex mixes. PCR products were added to 9μL of highly deionized formamide (HiDi) and 0.4µL of GeneScan-600 (LIZ) size standard (Applied Biosystems, Foster City, CA, USA) for capillary sequencing on an ABI Prism 3130XL Genetic Analyzer. Fragments were scored using Applied Biosystems Microsatellite Analysis Software (ThermoFisher Scientific Inc). 

Paternity Analyses 

Known maternal half-siblings were used for sibship reconstruction and paternity assignment with the maximum likelihood approach in COLONY v2.0.7.0. (Jones & Wang 2010). COLONY parameters included a polygamous mating system for both sexes; inbreeding; and a monecious, diploid species. A long run with medium-likelihood precision and a genotyping error rate of 1% was performed. Maternal and paternal genotypes were reconstructed using an allele probability threshold of 0.925 for allele calls at each locus. Seeds were categorized as selfed if the putative father had the same genotype as the putative mother. Outcrossing rates were calculated as the proportion of outcrossed offspring per meadow, shoot, rhipidium, and spathe. The effective number of sires and paternity skew per spathe, rhipidium, and shoot were calculated after Neff et al. (2008) in which effective sires = 1/Σ(rsi /seeds)2 where rsi = the number of offspring assigned to sire i, and the summation is over all sires contributing to a maternal spathe, rhipidium, or reproductive shoot. Skew was then expressed as 1 – (effective number of sires/actual number of sires). As such, a value of 0 implies no skew in which case all sires contribute equally to a seed set, and a value approaching 1 implies maximal skew in which case nearly all offspring are assigned to a single father. The paternity skew of each sire was calculated as the proportion of genotyped seeds per shoot sired by a particular sire.   

Using the reconstructed genotypes, average kinship (k), observed and expected heterozygosity (HO and HE, respectively) and clonality were calculated for the parent plants in each meadow using GENODIVE (Meirmans & Van Tienderen 2004). Following Iacchei et al. (2013), individuals were categorized by levels of kinship (k): ‘nearly identical’, 0.57 > k > 0.375; ‘full siblings’, 0.375 > k > 0.1875; ‘half siblings’, 0.1875 > k > 0.09375; and ‘quarter siblings’, 0.09375 > k > 0.047 (Loiselle et al. 1995). In addition, the heterozygosity of each paternal genotype was calculated as the proportion of heterozygous loci.  

Statistical Analyses  

Statistical analyses were conducted in RStudio with R v4.2.1 (R Core Team 2022, Posit Team 2023), and figures were generated with “ggplot2” and “lattice” (Sarkar 2008, Wickham et al. 2016). Data were tested for outliers, collinearity, even sample size, and normal distribution (Zuur et al. 2007). To assess patterns in seed viability, a generalized linear mixed model (GLMM) was fit to test the fixed effects of meadow, rhipidium position, and spathe position on the number of viable seeds per spathe (Poisson distribution, offset by the number of seeds per spathe). A random effect of maternal identity was added to account for differences among mothers. To assess patterns in mating system variation, GLMMs with the appropriate distribution were fit to test the fixed effects of meadow, rhipidium position, and spathe position on the response variables of number of outcrossed seeds (Poisson distribution), paternity skew (log+1 transformed, Gaussian distribution), and number of sires (Poisson distribution) per spathe using the package “lme4” (Bates et al. 2015). A random effect of seeds per spathe was added to control for expected variance due to the fair raffle process in sperm competition (Parker 1990, Zuur et al. 2007), and random effect of maternal identity was added to account for differences among mothers. Model residuals were visually inspected for normality and homogeneity using the package “DHARMa” (Hartig & Lohse 2022). Global models were used to perform ANOVAs (α = 0.05) and post hoc pairwise comparisons (α = 0.05) using the package “car” (Fox & Weisberg 2011).  

To explore whether the number of sires was influenced by (a) the number of genotyped seeds, (b) the total number of seeds, and (c) the proportion of selfed offspring, GLMMs (Poisson distribution) were also fit on both a ‘per spathe’ and ‘per shoot’ basis, with a random effect of maternal identity. To assess the impact of paternal genotype on reproductive success, GLMMs (Poisson distribution) were fit to test paternal heterozygosity, whether a sire selfed or outcrossed, and their interaction on (a) the number of seeds sired per male and (b) paternity skew per male, with a random effect of the number of successfully reconstructed loci per sire. Model residuals were visually inspected for normality and homogeneity using the package “DHARMa” (Hartig & Lohse 2022).