DOI: 10.5061/dryad.dfn2z358x
title: Mating system variation in the eelgrass, Zostera marina

description: 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. 
description: <em>Study Organism
 </em>  During
 <em>Zostera marina</em>’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
 <em>Z. marina </em>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).  
 <em>2.2. Study Sites & Sample
 Collection </em>  Three
 seagrass meadows located in Topsail Sound, North Carolina were sampled to
 characterize <em>Zostera marina </em>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
 <em>Z. marina</em> in both monospecific and mixed assemblages
 (i.e., co-occurring with <em>Halodule wrightii</em> and
 <em>Ruppia maritima</em>) (NCDEQ 2021). Sites were classified
 as shallow subtidal (depth &lt; 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.  
 <em>Seed
 Genotyping</em>  DNA was
 extracted from viable seed samples using a PowerPlant® Pro DNA Isolation
 Kit. Ten microsatellite loci previously characterized for <em>Z.
 marina </em>(Reusch et al. 1999, Reusch 2000a, Oetjen &
 Reusch 2007, Oetjen et al. 2010)<em> </em>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
 ddH<sub>2</sub>O. 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).  
 <em>Paternity Analyses</em>  
 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/Σ(<em>rs</em><em><sub>i</sub></em>
 /seeds)<sup>2</sup> where
 <em>rs</em><em><sub>i</sub></em> = the
 number of offspring assigned to sire <em>i</em>, 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 (<em>k</em>),
 observed and expected heterozygosity (H<sub>O </sub>and
 H<sub>E</sub>, 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 (<em>k</em>): ‘nearly
 identical’, 0.57 &gt; <em>k </em>&gt; 0.375; ‘full
 siblings’, 0.375 &gt; <em>k </em>&gt; 0.1875; ‘half
 siblings’, 0.1875 &gt; <em>k </em>&gt; 0.09375; and
 ‘quarter siblings’, 0.09375 &gt; <em>k </em>&gt; 0.047
 (Loiselle et al. 1995). In addition, the heterozygosity of each paternal
 genotype was calculated as the proportion of heterozygous
 loci.   <em>Statistical
 Analyses </em>  
 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).  
description: # Mating system variation in the eelgrass, Zostera marina
 [https://doi.org/10.5061/dryad.dfn2z358x](https://doi.org/10.5061/dryad.dfn2z358x) ## Description of the data Microsatellite genotypes of sampled offspring are available in the `Seed Genotypes` file. Parent genotypes were reconstructed using the program COLONY (Jones & Wang 2010) and are available in the `Parent Genotypes` file. Offspring genotypes were then associated to their respective seed on each parent plant and organized in the `Paternity` file. This file contains data organized on spatial scales corresponding to varying degrees of seagrass morphological structures - on a 'per shoot,' 'per rhipidium,' 'per spathe,' and 'per sire' basis. Data exploration was performed, in which outliers in both dependent and independent variables were identified and the appropriate distribution was assigned to each dependent variable (`Mating System Data Exploration Results`). Data were then quality controlled to remove any outliers and control for uneven sample size, and quality controlled (slim) data were added to the `Paternity` file. The `Paternity` file was subset into various .csv files for use in RStudio; specifically, `by dad dataset`, `by shoot dataset`, `by shoot dataset slim`, `by spathe dataset`, `by spathe dataset slim`, and `Slim viability by spathe` .csv files. Statistical models and data visualization was performed using the `Mating System` R code. Model outputs are available in the `Mating System Model Outputs` file. A subset of data was formed by rerunning COLONY using only offspring with 9-10 successfully amplified and scored loci. This dataset is referred to as the "reduced" dataset throughout. Parent genotypes were again reconstructed using the program COLONY (Jones & Wang 2010) and have been added to the `Parent Genotypes` file. Offspring genotypes were then associated to their respective seed on each parent plant and organized in the `Reduced Paternity` file. This file contains data organized on spatial scales corresponding to varying degrees of seagrass morphological structures - on a 'per shoot,' 'per rhipidium,' 'per spathe,' and 'per sire' basis. The Reduced Paternity file was subset into the `reduced per spathe` .csv file for use in RStudio. Statistical models and data visualization was performed using the `Mating System` R code. Model outputs are available in the `Mating System Model Outputs` file. In the `Paternal Heterozygosity` and `Reduced Paternal Heterozygosity` .csv files, heterozygosity was calculated with only reconstructed alleles that received a probability &gt;0.925 in COLONY. We formed another subset of data when investigating male heterozygosity by using all alleles (regardless of probability) for only males that yielded an average allele probability &gt;0.925. ## File structure `Seed Genotypes` 1. *Site*: Meadow from which each seed was sampled. See Methods and Supplementary Material for additional details on sampling sites. 2. *Maternal ID:* Reproductive shoot from which each seed was sampled. 3. *Rhipidium*: Rhipidium (branched structure) from which each seed was sampled. Please see Methods and Supplementary Material for additional details regarding the labeling of reproductive shoot structures. 4. *Spathe*: Spathe (branched structure within a rhipidium) from which each seed was sampled. Please see Methods and Supplementary Material for additional details regarding the labeling of reproductive shoot structures. 5. *Offspring Number:* Number that was assigned to a given seed within a given reproductive shoot. 6. *ZM Alelle_*: ZM1-10 refer to 10 microsatellite markers used in this study (please see Table A1 in the supplementary material), and Allele1-2 refers to the allele call for each diploid marker. 7. ****Please note that a blank cell denotes a marker that was unsuccessfully amplified and scored.* `Parent Genotypes` 1. *Sire*: Number assigned by COLONY to each identified sire. 2. *Prob*: Probability assigned by COLONY to each reconstructed parent locus. 3. *#Reconstructed Loci*: The number of loci that were reconstructed by COLONY. Any locus with a probability less than 0.925 was not considered to be successfully reconstructed. 4. *%Reconstructed Loci*: Calculated as the number of reconstructed loci, divided by that number of loci used in this study (10), multiplied by 100. 5. *Heterozygous Loci*: The number of reconstructed loci that are heterozygous. 6. *Heterozygosity*: Calculated as the number of heterozygous loci divided by the number of reconstructed loci. 7. *Average Probability*: The average locus reconstruction probability for each parent. 8. ****Please note that there are two tabs: 'Original Data' are genotypes reconstructed from that original COLONY run using all offspring genotypes. 'Reduced Data' are genotypes reconstructed from the COLONY run using only offspring with 9-10 successfully reconstructed loci.* `Paternity` & `Reduced Paternity` 1. *Shoots*: The number of reproductive shoots sampled per site. 2. *Sires*: The number of sires detected by COLONY across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 3. *Seeds*: The number of seeds sampled across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 4. *Viable Seeds*: The number of viable seeds sampled across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 5. *Genotyped Seeds*: The number of seeds with at least one successfully amplified and scored microsatellite locus across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 6. *%Genotyped*: Calculated as the number of genotyped seeds divided by the number of sampled seeds across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 7. *Outcrossed Seeds*: The number of outcrossed (i.e., not self-fertilized) seeds across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 8. *Selfed Seeds*: The number of self-fertilized seeds across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 9. *Outcrossing Rate*: Calculated as the number of outcrossed seeds divided by the number of genotyped seeds. 10. *Skew*: Paternity skew, calculated after Neff et al. (2008) across various spatial structures including sites, reproductive shoots, rhipidia, and spathes. 11. *Shoots Fertilized*: The number of reproductive shoots a particular sire was detected as fertilizing. 12. *Rhipidia Fertilized*: The number of rhipidia a particular sire was detected as fertilizing. 13. *Spathes Fertilized*: The number of spathes a particular sire was detected as fertilizing. 14. *Seeds Fertilized*: The number of seeds a particular sire was detected as fertilizing. 15. *Skew per Sire*: Calculated as the number of seeds fertilized by a particular sire divided the number of genotyped seeds on the shoot he fertilized. `by spathe dataset`, `by spathe dataset slim` & `reduced by spathe` 1. *loggs*: The number of genotyped seeds per spathe, log transformed. 2. *logSkew*: Paternity skew, log transformed. 3. ****Please note that a blank cell denotes data that is unavailable, as paternity skew is only applicable when there are more than one offspring within a spathe.* `Slim viability by spathe` 1. *logseeds*: The number of seeds per sampled spathe, log transformed. 2. *Proportion*: Calculated as the number of viable seeds per spathe divided by the number of seeds per spathe. `Paternal Heterozygosity` & `Reduced Paternal Heterozygosity` 1. *Reproductive Strategy*: Whether a sire fertilized a seed via outcrossing or self-fertilizing. ## Code/Software COLONY v2.0.7.0 (Jones & Wang 2010) was used to assign parent genotypes to each offspring genotype, with which individual heterozygosity, outcrossing rates, paternity skew (Neff et al. 2008), and reproductive skew were calculated. GENODIVE (Meirmans & Van Tienderen 2004) was used to calculate population expected heterozygosity, observed heterozygosity, and kinship coefficients with reconstructed parent genotypes. R is required to run Mating System Data Exploration and Mating System script, which were created using R v4.2.1. Annotations are provided throughout the script through 1) library loading, 2) dataset loading, 3) statistical analysis, and 4) figure creation. ## Supplementary Material Supplementary material including Table A1 (Microsatellite forward and reverse primers used in this study), Figure 1 (map of sampling sites), and Figure 2 (labeling system for eelgrass reproductive shoots) have been included in the Figure1_Figure2_TableA1.pdf file. Please see Materials and Methods for additional context regarding the supplementary material. \****All .xlsx and .csv files can be found in the Dryad deposit; supplementary material and R code can be found in the linked Zenodo dataset.*