Maternal quality, paternal effects, and sibling interactions influence seed size in the eelgrass, Zostera marina

Seed size is an essential determinant of germination and survival in angiosperms. Zostera marina, one of few marine angiosperms, is a key foundation species present in coastal marine ecosystems, and edge-of-range population persistence is increasingly reliant on seed production. While environmental conditions have been invoked to explain regional patterns of seed size variation, far less is known about seed size differences within a plant. In this study, we genotype and measure individual seeds across parent plants to investigate the relative contributions of maternal, paternal and offspring traits on seed weight. Specifically, we investigate how parent heterozygosity, paternal siring success, outcrossing, a size-number tradeoff, and sibling interactions influence both seed weight (mg) and cross-sectional area (mm2). Here, we provide all necessary files to conduct the entirety of analyses used in the manuscript. 

Methods

1. Study Sites and Sample Collection  

Three seagrass meadows in the intracoastal waters of Topsail Sound, North Carolina were sampled to characterize Zostera marina seed size variation. Meadows were on a narrow shelf classified as shallow subtidal (depth < 2m MLLW) between the Intracoastal Waterway and the adjacent shoreline (34.22 N, 77.37 W; Fig. S1). On 4 May 2021, near the end of the eelgrass reproductive season, flowering shoots were haphazardly collected at least 5m apart, yielding 35 shoots total. Samples were transported on ice to UNCW’s Center for Marine Science, and the morphological characteristics of each shoot were recorded. Specifically, each rhipidium (branching reproductive structure), spathe (seed-containing branch within a rhipidium), and seed position was labeled in order of decreasing proximity to the rhizome (i.e., basal positions were given a value of 1; Fig. S2). 

2. Seed Measurements and Genotyping 

Seeds were removed from spathes, blotted dry, and tested for viability using the “squeeze test” by gently compressing individual seeds with a pair of tweezers (Marion and Orth 2010). Those with a seed coat that compressed were considered nonviable. Viable seeds were then weighed (mg) with a microbalance and photographed at ×10 magnification on a stereo microscope (Fisher Scientific, MA, USA). Length (mm), width (mm), and cross-sectional area (mm2) were measured for each seed using ImageJ software (Schneider et al. 2012). Cross-sectional area of seeds was then calculated from the following equation: v = πab, where a = radius of length and b = radius of width (Wyllie-Echeverria et al. 2003). DNA was extracted from viable seed samples using a PowerPlant® Pro DNA Isolation Kit (QIAGEN, Hilden, GER). Ten microsatellite loci previously described for Z. marina (Reusch et al. 1999, Reusch 2000, Oetjen and Reusch 2007, Oetjen et al. 2010; Table S1) 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 (Table A1). PCR conditions for all multiplex conditions were as follows: 95.0°C for 15 min; 2 cycles of 94.0°C for 15 s, 60.0°C for 30 s, 72.0°C for 45 s; 2 cycles of 94.0°C for 15 s, 59.0°C for 30 s, 72.0°C for 45 s; 2 cycles of 94.0°C for 15 s, 58.0°C for 30 s, 72.0°C for 45 s; 2 cycles of 84.0°C for 15 s, 57.0°C for 30 s, 72.0°C for 45 s; 28 cycles of 94.0°C for 15 s, 56.0°C for 30 s, 72.0°C for 45 s; and a final 2 min 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.)

3. Paternity Analyses 

Paternity assignment and parent genotype reconstruction were performed using all offspring genotypes in COLONY v2.0.7.0 (Jones & Wang, 2010; Wang, 2019). COLONY input parameters included a polygamous mating system for both sexes, possible inbreeding, and a monecious, diploid species. The maternal plant of each seed (i.e., in which seeds are at least maternal half-siblings) was specified. The maximum likelihood approach with a long run, medium-likelihood precision, and a genotyping error rate of 1% was performed. Both maternal and paternal genotypes were reconstructed using the maximum likelihood approach and a probability threshold of 0.925 to infer the most likely parental genotypes based on observed sibship patterns. Seeds were categorized as selfed if the putative father had the same multi-locus genotype as the known mother. It was also noted whether a seed's basal neighbor (i.e., the seed closer to the base of the plant) was fertilized via outcrossing or selfing. Paternity skew was used as a proxy for within-spathe seed relatedness and calculated following Neff et al. (2008). A value of 0 indicates no skew in which all sires contribute equally to seeds within a spathe (seeds would be some combination of half- and full-siblings), and a value approaching 0.5 indicates maximum skew in which all seeds were sired by a single father (all seeds would be full-siblings). Paternal siring success was defined as the extent to which a sire monopolized available ovules on a plant and was measured as the proportion of viable seeds per shoot sired by a particular father. Heterozygosity of the putative sires (as well as the maternal shoots and individual seeds) was calculated as the proportion of heterozygous alleles across successfully reconstructed loci.  

2.5 Statistical Analyses  

Statistical analyses were conducted in RStudio with R v4.2.1 (R Core Team, 2022; Posit Team, 2023). Data were tested for outliers, collinearity, and normal distribution (Zuur et al. 2007; see Figure S3 in the supporting information for the relationships among fixed effects). All model residuals were visually inspected for normal distribution (DHARMa; Hartig & Lohse, 2022), and all figures were generated using the package “ggplot2” (Wickham et al., 2016).  

To test the fixed effect of maternal heterozygosity on reproductive shoot characteristics, generalized linear models (GLMs) were fit to the total number of seeds, viable seeds (offset by the number of seeds), spathes, and rhipidia per shoot using a Poisson distribution as well as to the mean and CV of seed weight and area per shoot using a gamma distribution and a log link function (lme4; Bates et al., 2015). To test the size-number tradeoff, GLMs were fit between the proportion of viable seeds, the total number of seeds, and mean seed weight per shoot using a gamma distribution and a log link function. 

To first test the fixed effects of rhipidium position, spathe position, and seed position on individual seed size, a GLM was fit to seed weight and area using a gamma distribution and a log link function (lme4; Bates et al., 2015). This model was also fit to reproductive mode (selfed vs. outcrossed), as previous work has shown that outcrossing decreases in more distal spathes (Sgambelluri et al. 2024). An ANOVA and Tukey's Honest Significant Difference test was performed on the model outputs to assess the significance of each fixed effect (car; Fox & Weisberg, 2019). 

To then test the fixed effects of maternal heterozygosity, paternal heterozygosity, offspring heterozygosity, outcrossing, paternal siring success, within-spathe seed relatedness, basal neighbor effects, and the proportion of viable seeds per shoot on individual seed size, generalized linear mixed models (GLMMs) were fit to weight and area using a gamma distribution and a log link function (lme4; Bates et al., 2015). Here, the random effects of rhipidium and spathe position were included to account for observed positional differences (see Results) and a random effect of maternal identity was added to account for possible differences among mothers. Models were run with different combinations of fixed and random effects, and those garnering the highest Akaike weight were considered top models for inference (AICcmodavg; Mazerolle, 2023). To then evaluate the relative impact of maternity, paternity, and developmental timing on seed size, an REML variance components analysis (VCA) was performed with spathe position, maternal identity, and paternal identity as random effects and with individual seed weight as the response variable.  

Finally, to test the fixed effects of maternal heterozygosity, the number of seeds per spathe, the proportion of selfed seeds per spathe, and within-spathe seed relatedness on variance in seed size, GLMMs were fit to the mean and the coefficient of variation (CV = σ/µ, where σ = standard deviation and µ = sample mean) in seed weight and area within spathes.