Repeated Genetic and Adaptive Phenotypic Divergence across Tidal Elevation in a Foundation Plant Species

Microgeographic genetic divergence can create fine-scale trait variation. When such divergence occurs within foundation species, then it might impact community structure and ecosystem function and cause other cascading ecological effects. We tested for parallel microgeographic trait and genetic divergence in Spartina alterniflora, a foundation species that dominates salt marshes of the US Atlantic and Gulf coasts. Spartina is characterized by tall-form (1–2 m) plants at lower tidal elevations and short-form (<0.5 m) plants at higher tidal elevations, yet whether this trait variation reflects plastic and/or genetically differentiated responses to these environmental conditions remains unclear. In the greenhouse, seedlings raised from tall-form plants grew taller than those from short-form plants, indicating a heritable difference in height. When we reciprocally transplanted seedlings back into the field for a growing season, composite fitness (survivorship and seed production) and key plant traits (plant height and biomass allocation) differed interactively across origin and transplant zones in a manner indicative of local adaptation. Further, a survey of single nucleotide polymorphisms revealed repeated, independent genetic differentiation between tall- and short-form Spartina at five of six tested marshes across the native range. The observed parallel, microgeographic genetic differentiation in Spartina likely underpins marsh health and functioning and provides an underappreciated mechanism that might increase capacity of marshes to adapt to rising sea levels.


Introduction
Local adaptation is widespread in nature (Leimu and Fischer 2008;Hereford 2009;Savolainen et al. 2013) and can underpin patterns in the distribution and abundance of organisms (Schoener 2011), generate high genetic variation among populations (Bolnick et al. 2011), and affect the responses of populations to both disturbance and climate change (Valladares et al. 2014). Adaptive phenotypic divergence can occur across spatial scales finer than a species' dispersal range (i.e., microgeographic adaptation; Kawecki and Ebert 2004;Roda et al. 2013;Richardson et al. 2014), a pattern that has been regularly documented within plant species. For instance, plants can adapt to variability in edaphic conditions (e.g., serpentine soils: Kruckeberg 1967;Brady et al. 2005;mine tailings: McNeilly 1968;Antonovics and Bradshaw 1970;Antonovics 2006) and physical factors (e.g., emersion time : Hays 2007; hydrodynamics: Roberson and Coyer 2004;soil moisture: Scotti et al. 2016) over spatial scales that are smaller than the scale of dispersal. Such microgeographic adaptation can evolve repeatedly across independent populations, each of which responds to similar local gradients in the environment, resulting in the evolution of similar traits (i.e., parallel adaptive phenotypic divergence; Rundle et al. 2000;Nosil et al. 2002;Johannesson et al. 2010).
Although many plant species that exhibit microgeographic adaptation are widely distributed, few studies have tested whether this divergence occurs repeatedly across sites (but see Roda et al. 2013;Bertel et al. 2018;Brousseau et al. 2020). Parallel adaptive phenotypic divergence can arise through two main processes in the presence of gene flow: (a) initial adaptive divergence occurs once and is followed by colonization of differently adapted forms into distinct habitats at multiple sites, or (b) divergence between habitats occurs repeatedly at multiple sites following colonization (Schluter and McPhail 1992;Johannesson et al. 2010;Butlin et al. 2013). This repeated evolution could follow many nonparallel genetic pathways and still result in similar phenotypic differentiation (Johannesson et al. 2010;Bolnick et al. 2018;Van Belleghem et al. 2018). Thus, similar phenotypes occurring in similar environments across independent populations strongly suggest natural selection, as nonadaptive processes are unlikely to yield the same evolutionary shifts that are correlated with the environment (Conte et al. 2012;Bertel et al 2018;Bolnick et al 2018). However, a recent synthesis of bacteria and fungi experiments illustrated that nonadaptive processes (e.g., mutation) can also contribute to parallel evolution (Bailey et al. 2017). Given that genetic divergence prior to colonization can also give rise to parallelism (Butlin et al. 2013), observations of similar phenotypeenvironment associations are not sufficient to infer parallel evolution. Instead, parallel adaptive divergence requires evidence of local adaptation in multiple locations (Bertel et al. 2018) that evolved after demographic separation (Langerhans and DeWitt 2004;Butlin et al. 2013;Savolainen et al. 2013).
One promising scenario for generating parallel adaptive phenotypic divergence arises when species persist across strong environmental gradients that occur repeatedly throughout their distribution (Bailey et al. 2015), such as wetland plant species that span an inundation gradient. Prior studies have documented small-scale genetic or trait divergence in such species (e.g., Puccinellia: Gray 1985; Phragmites: Gao et al. 2012;Schoenoplectus: Sweetman et al. 2013), including examples of genetic differentiation across the land-to-water gradient (e.g., Salicornia: Davy et al. 1990;Phragmites: Clevering 1999a, 1999b. In addition, there is evidence for both genetic variation among individuals and phenotypic plasticity in the response of wetland species to environmental factors including flooding (Lessmann et al. 1997;Clevering and Hundscheid 1998;Howard and Rafferty 2006), salinity (Hester et al. 1998(Hester et al. , 2001Howard 2010;Grewell et al. 2016), and nutrient availability (Clevering 1999a(Clevering , 1999b. For example, Spartina patens exhibited adaptive genetic and phenotypic differences across distinct sand dune, swale grassland, and marsh microhabitats (Silander 1979(Silander , 1984(Silander , 1985Silander and Antonovics 1979). In combination, these studies suggest strong potential for microgeographic adaptation in wetland plants to environmental differences across an inundation gradient, yet there have been few comprehensive tests of this hypothesis (but see Knight and Miller 2004;Lenssen et al. 2004) and no tests of parallel divergence to inundation across a species' range.
The plant species Spartina alterniflora (hereafter, Spartina; Poaceae) is arguably the most important foundation species (sensu Ellison 2019) in estuaries of the Atlantic and Gulf coasts of the United States. Monocultures of Spartina create millions of acres of salt marshes, critical estuarine ecosystems that provide shoreline protection, water purification, carbon sequestration, and habitat for economically valuable and ecologically important species (Barbier et al. 2011). Within most marshes, Spartina exhibits pronounced variation in stem height over tens of meters across a natural stress gradient in tidal inundation (Pennings and Bertness 2001), which is within the distance that Spartina pollen and seeds can disperse (Taylor et al. 2004;Travis et al. 2004;Blum et al. 2007). Tall-form Spartina (~1-2 m stem height) is found in low elevations along creek banks that experience daily inundation, and short-form Spartina (!0.5 m) is found at higher elevations that are less consistently flooded (fig. 1A, 1B; Shea et al. 1975;Valiela et al. 1978). Early experiments provided evidence both for phenotypic plasticity (Shea et al. 1975;Mendelssohn and McKee 1988) and genetic differentiation of plant height (Stalter and Baston 1969;Gallagher et al. 1988) in Spartina (reviewed in Anderson and Treshow 1980). For instance, long-term fertilization of field plots shifted plant morphology from short to tall (Valiela et al. 1978;Morris et al. 2002), which was interpreted as supporting environmental control, though nutrient availability could have also selected for these trait differences. Alternatively, both reciprocal transplant (Stalter and Baston 1969) and common garden (Gallagher et al. 1988) experiments conducted with adult transplants detected consistent morphological differences across the two growth forms over a single growing season, indicating a genetic basis. Yet these studies could not eliminate the potential confounding effects of irreversible development plasticity often exhibited in adult plants (Callaway et al. 2003). In addition to experiments, population genetic surveys can generate estimates of historical demography and thus provide context to understand the evolution of phenotypic divergence. Previous work using allozyme and amplified fragment length polymorphism technologies did not detect genetic differentiation between tall-and short-form Spartina (Shea et al. 1975;Valiela et al. 1978;Richards et al. 2004;Foust et al. 2016), suggesting that weak to no barriers to gene flow exist among these forms. To date, no analyses have been done with highly variable single nucleotide polymorphisms (SNPs), which are better suited to detect low levels of genetic differentiation.
We combined field and greenhouse experiments of seedreared individuals with a field survey of genetic differentiation using SNPs to examine genetic and trait differentiation across both biogeographic regions and marsh elevations. By doing so, we evaluated whether Spartina exhibits microgeographic adaptation to tidal elevation and examined whether this divergence occurred in parallel in Spartina populations across its native North American range. We focused on fitness metrics including survival and reproduction (e.g., flower production and seed set) as well as growth responses (e.g., density, height, and biomass). Biomass production and allocation are important components of long-term fitness in perennial clonal plants (Younginger et al. 2017). In addition, stem height can influence pollen dispersal in wind-pollinated grasses such as Spartina (Van Boheemen et al. 2019), which may impact fitness. Given that Spartina morphology can influence both abiotic (e.g., sediment accumulation, soil erodibility, and light availability; Morris et al. 2002;Seliskar et al. 2002;Bernik et al. 2018) and biotic (e.g., associated species composition; Seliskar et al. 2002) marsh characteristics, determining the degree to which these two Spartina phenotypes are genetically determined and reflect evolutionarily distinct entities has broad implications for predicting how salt marshes will respond to sea level rise (Kirwan et al. 2010), as well as for informing marsh restoration and management efforts (McKay et al. 2005).

Common Garden Experiment
We germinated seeds collected in the field from tall-and short-form Spartina at three sites in South Carolina separated by at least 3 km ( fig. A1; app. 1; figs. A1-A12 and apps. 1-4 are available online) and grew them in a greenhouse common garden for 11 months to assess trait differences. Although using first-generation seedlings cannot completely eliminate maternal effects (Roach and Wulff 1987;Bischoff and Muller-Scharer 2010), because few F 1 plants flowered in the subsequent 2 years of our study, breeding additional generations in the greenhouse was not feasible. Inflorescences containing seeds were collected from the field in October 2014 from 20 reproductive stems per marsh, 10 each from characteristically tall-and short-form Spartina plants (N p 60 total). We also characterized the surrounding natural marsh by measuring Spartina stem density and height and sediment redox (millivolts) in the tall and short zones at these three marshes in May 2015 (see app. 3 for methods and results). Seeds were stored in the dark at 47C to simulate the natural dormant period  and then transferred to the Northeastern University Marine Science Center (MSC) greenhouse in January 2015 to stimulate germination (app. 1). In April 2015, we planted 25 individual seedlings from each site (N p 3) by zone (N p 2) combination in 6" diameter pots (volume p 2.9 L) in a 50∶50 sand and potting soil mix. Pots were randomly assigned to one of two water tables in the MSC greenhouse. All pots were irrigated daily with freshwater and submerged in seawater for 8 h weekly; we note that these conditions more closely mimic the environmental conditions in the short-form Spartina zone. We measured survival and density as well as height and leaf number of the tallest stem in each pot monthly during the growing season (April-October 2015; density measurements began in June, and leaf number was not measured in June) and then bimonthly through the following winter and spring (November 2015-March 2016. At the end of the experiment, we also measured the height of four additional randomly selected live stems to estimate average stem height. We were unable to measure above-and belowground biomass, a destructive response variable, at the end of the common garden experiment, as many of these individuals were then used in the reciprocal transplant experiment (see below). As noted above, the use of first-generation seedlings reduces, but does not eliminate, maternal effects. To test for this potential confounding effect, we measured seed size, a commonly cited trait that can influence germination and early seedling survival and growth (Roach and Wulff 1987;Kalisz 1989). Using field-collected seeds that had not germinated by May 2015, we dried tins with five seeds for each site#zone combination (N p 7-8 per combination) in a drying oven at 607C for at least 48 h. Total dry weight was divided by the number of seeds measured to obtain an average dry weight (grams) per individual seed. While ungerminated seeds may be a biased subset of the total seed pool, as size could have contributed to why these seeds did not germinate, we have no reason to expect this would vary consistently between tall-and short-form Spartina.
In addition, we note that we cannot rule out the potential role of maternal effects that are unrelated to seed weight (e.g., nutrient provisioning, seed architecture, and epigenetic changes including DNA methylation ;Rossiter 1996;Richards et al. 2017).

Reciprocal Transplant Experiment
We ran an 8-month factorial reciprocal transplant experiment at the collection sites using greenhouse-reared individuals. We measured a suite of morphological and fitness metrics both during and at the end of the experiment (app. 2). In March 2016, we selected 10 genets (i.e., unique genotypes as confirmed by DNA microsatellites developed for Spartina; Blum et al. 2004;Sloop et al. 2005) from each site by zone combination in the common garden experiment for our reciprocal transplant experiment (N p 60 genets; fig. A1). Six clonal ramets from each genet (i.e., 360 transplants total) were isolated into separate 6" diameter plastic pots in a 50∶50 sand∶soil mix, with the pot bottom removed and lined with window screen to allow for water transfer and soil interactions in the field. This setup eliminated the ability of our transplants to directly interact with the surrounding plant neighborhood in the field except for any shading effects, thus isolating the effects of the environment in the absence of neighbor effects. Transplants were irrigated daily with freshwater, fertilized biweekly using Miracle-Gro water soluble all-purpose plant food (100 mL per pot; NPK at 24-8-16 concentration and trace elements), and submerged in seawater twice weekly for 8 h for 1 month prior to transplanting in the field.
In May 2016, we transported all transplants overnight by truck to Charleston and planted them into the field over the course of 2 days. We established six transplant locations adjacent to each of the seed collection areas. We planted each genotype (N p 10 per origin site#zone) in each transplant location (N p 60 individuals total per transplant location). This design allowed us to control for potential effects of genetic identity on the response variables of interest (Zerebecki et al. 2017). Within each transplant location, we set up two parallel 30-m transects separated by 2-5 m. We randomly assigned and planted 30 transplants flush with the sediment surface at 1-m intervals along each transect (N p 5 per origin site#zone per transect; see fig. A1). We left intact natural stands of Spartina surrounding our potted transplants within each transect; thus, in the tall zone, transplanted seedlings were planted under the canopy of naturally occurring tall-form Spartina (see app. 3 for additional details of the plant characteristics at our field sites).
We measured transplant survival, live stem density, and height of the tallest stem (i.e., maximum stem height) at approximately monthly intervals from June through December 2016. In addition, we quantified three metrics of Spartina reproductive output throughout the flowering period (September-December 2016, sampled approximately every 2 weeks): proportion of transplants that flowered, total number of flowering stems per pot, and average number of seeds per flowering stem. At the end of the experiment, we measured the height of all live stems, harvested the plants, separated the biomass into aboveground (stems and leaves) versus belowground (roots and rhizome) portions, and then dried samples at 607C for at least 48 h to calculate total dry biomass as well as examine variation in above-and belowground allocation.

Population Genetics Field Survey
Logistical constraints limited our common garden and reciprocal transplant experiments to evaluate plants from multiple sites at a single geographic location. We complemented these results with a broader genetic field survey to more completely examine adaptation to zonal location in the marsh in Spartina and to evaluate the generality of our experimental findings. To do so, we generated an SNP survey to assess genetic differentiation at two scales: across biogeographic regions and between tall and short forms within marshes (table A6; tables A1-A6 are available online), including two of the three marshes at which we conducted the experiment. At six marshes across Massachusetts (MAS, MAW), Rhode Island (RI), South Carolina (SCBI, SCFB), and Florida (FL), we collected plants from the tall zone (T) within approximately 2 m of creeks or open water and from the short zone (S) at least 20 m away. We also collected tissue within the short zone at only one marsh (SCFJ), the tall zone at only two marshes (SC2, SC3), and from unidentified zones at two marshes (North Carolina [NC1] and New Hampshire [NH1]; table A6). At all locations, ∼10 cm 2 of green leaf tissue was collected from 20 plants at least 1 m apart along a linear transect and placed into silica for drying. Given high genotypic diversity within marshes (Richards et al. 2004;Hughes and Lotterhos 2014;Foust et al. 2016), the probability that plants 1 m apart are clonal replicates is low.
To generate genotypes, we entered all sample extracts into a single genomic library by following the protocols of Parchman et al. (2012) and used Illumina HiSeq technology to generate 249:2#10 6 reads approximately 80-100 base pairs in length. We then created a de novo assembly using a subset of these reads with SeqMan NGen, aligned reads from individuals to this assembly using a Burrows-Wheeler Aligner tool (bwa, ver. 0.7.12; Li and Durbin 2009), and generated genotype likelihoods with samtools/bcftools 1.3.1 (-mpileup and -call protocols) with the multiallelic calling option . After filtering, we generated genotype likelihoods for 2,735 SNPs across 309 individuals (n p 13-21 per zone per marsh; table A6), with 5:2050:07 reads (mean5 standard error) per SNP per sample. More details on DNA extraction, library preparation, and in silico processing and filtering are found in appendix 4.

Data Analysis
Common Garden and Reciprocal Transplant. We included three sites in our study to evaluate the consistency of morphological differences across tidal elevations and the potential for parallel adaptive phenotypic divergence to zone. Because we were not interested in specific characteristics of each site, we treated origin site and transplant site (when applicable) as independent fixed block effects in our experimental analyses. We did not design our experiment to robustly test all possible interactions between origin site and transplant site.
Morphological responses in the common garden experiment (app. 1) were analyzed with a linear model using origin zone and origin site as fixed factors. Further, we also assessed whether seed weight, a proxy for maternal investment, was correlated with differences in maximum seedling height (see "Results") by conducting a linear regression on final average maximum height and average seed weight from each site # zone combination.
In the reciprocal transplant experiment, we assessed survival over time using a Cox proportional hazard regression with origin zone, transplant zone, and their interaction as fixed effects. Transplants that did not survive were then removed from further analyses. We conducted linear models with fixed effects of origin zone, transplant zone, and their interaction, as well as origin site and transplant site on all other response variables (height, density, etc.) at the end of the experiment. Flowering occurrence was analyzed using a generalized linear model with a binomial response. Local adaptation was indicated by a significant twoway interaction of origin and transplant zones with the local versus foreign criteria (i.e., higher fitness of local individuals than foreign) satisfied in post hoc comparisons (Kawecki and Ebert 2004). Statistical methods and results of morphological metrics measured repeatedly throughout the reciprocal transplant experiment are detailed in appendix 2.
To assess local adaptation in composite fitness (i.e., survivorship and total seed production), we used the R aster package (Geyer et al 2019;following Geyer et al. 2007 andShaw et al. 2008 and as implemented by Samis et al. 2016). We modeled survival with a Bernoulli distribution, presence or absence of seed set with a Poisson distribution, and number of seeds produced as a truncated Poisson distribution. We assessed fixed effects of transplant and origin site, origin zone (tall/short), transplant zone (tall/short), and the origin zone by transplant zone interaction using likelihood ratio tests with the aster function. For graphical purposes, we conducted univariate analyses of the effect of origin zone separately for tall and short transplant zones and generated a mean estimate of composite fitness.
All experimental analyses were conducted in R statistical software (ver. 3.6.3; R Core Team 2018). We used the lm function in the car package (Fox and Weisberg 2019) for linear models. Tukey's post hoc tests were performed using the lsmeans (Lenth 2016), pbkrtest (Halekoh and Højsgaard 2014), and multcompView (Graves et al. 2015) packages. Cox proportional hazards regressions were performed and visualized using the survival (Therneau 2015) and survminer (Kassambara and Kosinski 2018) packages.
Population Genetic Analysis. We examined genetic differentiation across and within biogeographic realms with several analytical methods. With noted exceptions, all analyses were based on genotype likelihoods to incorporate the uncertainty that arises from errors in sequencing and sampling (Nielsen et al. 2012). We estimated expected (H E ) and observed heterozygosity (H O ) using angsd (Korneliussen et al. 2013). We compared estimated patterns of heterozygosity across latitude using one-way ANOVAs and between tall and short zones using a paired t-test.
For the principal components analysis (PCA), we converted the phred-scale genotype likelihoods (from samtools/ bcftools) per SNP-sample combination into probabilities that summed to 1 and then converted these to a single value ranging from 0 to 2, where 0, 1, and 2 represent the highest probability of a homozygote, heterozygote, and alternative homozygote, respectively. We assessed differences among samples at two hierarchical spatial scales: across all marshes and between tall and short zones within each of six marshes, using the prcomp function in R. We used analysis of similarity (ANOSIM) implemented in vegan (Oksanen et al. 2019), a Bray distance matrix, and 1,000 permutations to assess differences between tall-and short-zone plants in marsh-specific analyses. We note that we repeated these PC analyses with a subset of loci in which loci within the upper 10% quantile of loadings on PC were removed.
We confirmed these PC analyses using hard-called genotypes as an input. Spartina alterniflora is hexaploid derived from an ancient allopolyloidization event (Ainoche et al. 2012). Previous work with Spartina microsatellites indicated inheritance patterns of alleles consistent with genome diploidization (Blum et al. 2007), and our genotype likelihood similarly treated genotypes as diploid. Recently, Clevenger et al. (2015) recommended that RADseq studies of polyploids utilize the approach of Stacks (Catchen et al. 2013). As suggested by Paris et al. (2017) for polyploids, we used a lower M (e.g., M p 1; number of mismatches allowed between two alleles of a heterozygote sample) and m p 3 (number of identical reads required to initiate a putative allele). After filtering (i.e., r80 p those loci found in 80% of samples or more), we generated 391 variant sites that we analyzed with PCA.
We implemented maximum-likelihood admixture methods to infer ancestry of individuals (NGSadmix; Meisner and Albrechtsen 2018) of eight data sets: all marshes, only marshes with both tall-and short-form Spartina samples (N p 6), and tall versus short form for each of the six marshes independently. For the first two data sets, genotype likelihoods (from samtools/bcftools) were run across an a priori range of genetic clusters (k p 2-20) and replicated five times. An analysis of delta k (as suggested by Evanno et al. 2005) indicates that the best-fit k to the data set with all marshes and with six marshes was k p 5 ( fig. A8) and k p 4 (fig. A12), respectively. For the marshspecific admixture analyses, genotype likelihoods were run across k p 1-10. For five of six marshes, the best-fit k was k p 2; when comparing MAWS and MAWT, the best-fit k was k p 3.
To assess how genetic variation is partitioned among populations (i.e., marshes) and among subpopulations (i.e., zones) within marshes, we performed hierarchical analysis of molecular variance (AMOVA; Excoffier et al. 1992) using pegas for three subsets of the data: all 12 subpopulations at six marshes in Florida, South Carolina, Rhode Island, and Massachusetts; four subpopulations at two South Carolina marshes; and four subpopulations at two Massachusetts marshes (for details, see app. 4). We also calculated pairwise F ST among and within marshes using ngsFST ) and compared pairwise F ST among populations categorized into three bins of geographic distance: short versus tall within a marsh (!0.2 km) and marshes separated by 0.5-5 km and 5-15 km (see app. 4).
We conducted a preliminary redundancy analysis (RDA; implemented in R package vegan; Oksanen et al. 2019) to identify whether there were loci that best correlated with genomic divergence between ecotypes. For comparison purposes, we also isolated the 20 loci that had greatest loadings in the correspondence analysis for each marsh independently, and we compared the identity of these loci across marshes.
To generate imputed genotypes for this RDA, we first compared the log likelihood of the three genotypes of each SNP by individual combination. We called a genotype when its log likelihood was greater than two times the log likelihood of any other genotype. If no genotype reached that threshold, then no call was made. For each marsh, we imputed missing genotypes by using the most common genotype.

Results
After 11 months in a greenhouse common garden, tallorigin seedlings had~10% greater maximum height compared with short-origin seedlings, despite similar heights at the start of the experiment (repeated measures ANOVA zone # date: P p :005; fig. 1C; table A1). Average stem height was also greater at the end of the common garden experiment (app. 1). Thus, the height difference between talland short-form Spartina likely has a heritable component.
We did not find a significant relationship between fieldcollected seed weight-a potential proxy for maternal investment-and plant height in the common garden experiment ( fig. A2A; F 1, 4 p 0:43, P p :54). Further, the weight of ungerminated seeds did not differ among origin zones (global mean p 1:79 mg; F 1, 43 p 0:007, P p :93; fig. A3). Survival, stem density, and leaf number also did not consistently vary between tall-and short-origin seedlings (for details, see app. 1).
At the end of our 8-month reciprocal transplant experiment, there was a significant interaction between origin zone and transplant zone on fitness (composite metric of survivorship and seed production) in the aster analysis ( fig. 2; table 1; deviance p 6:60; P p :010). Short-origin seedlings had higher fitness in the short zone than tall-origin seedlings (P p :005), whereas tall-origin seedlings tended to have slightly higher fitness in the tall zone than did short-origin seedlings (P p :101; fig. 2). This marginal effect likely occurred because we found only a single seed-producing stem in the tall transplant zone (i.e., a tall-origin plant from Bowens Island that flowered in the tall zone at Fort Johnson), but this individual transplant produced a large number of seeds (i.e., 190).
We also detected an interaction between origin zone and transplant zone on plant height (P p :006; fig. 3A; table 1) and the ratio of above-to-belowground biomass (P p :02; fig. 3B). In the tall zone, tall-origin transplants were significantly taller ( fig. 3A) and allocated more biomass aboveground ( fig. 3B) than short-origin transplants. In contrast, there were no significant trait differences between tall-and short-origin Spartina when transplanted to the short zone.
There were also clear transplant-zone effects on survival, growth, and reproduction when examined independently. For instance, survival of Spartina transplants was dramatically higher in the short than tall zone (92% vs. 29%, respectively; P ! :001; fig. 3F). In addition, transplants produced more flowering stems in the short than tall zone (18% vs. 2%; P p :007; fig. 3E; app. 2). Similarly, both stem den-sity and total biomass were significantly greater in the short than tall zone (density and total biomass: P ! :001; fig. 3C, 3D). Together, these results (figs. 2, 3) highlight that Spartina seedlings had lower fitness in the tall-zone environment compared with the short-zone environment.

Population Genomics
Our analysis of 2,735 SNPs from 11 marshes from the Atlantic and Gulf coasts of the United States identified strong biogeographic differentiation that reflects historical separation of Spartina across this range. For example, PCA de- ). An admixture analysis indicated that five genetic clusters (k p 5) best explained variation among genotypes; these clusters broke across similar phylogeographic breaks (FL, SC, NC, RI, MA1NH; fig. A8). These results mirror previous studies using chloroplast and microsatellite markers (Blum et al. 2007).
Expected heterozygosity (H E ) within southern marshes (FL, SC, and NC) was greater than within northern marshes (RI, MA, and NH) when all individuals were considered (average H E p 0:348 and 0.328, respectively; F 1, 15 p 16:5,  Note: F-statistics are for responses analyzed with linear models and x 2 for logistic regression. Transplant zone is location (tall vs. short) where planted; origin zone is location (tall vs. short) where maternal plant was from. Numerator degrees of freedom is 1 for each zone factor and 2 for each site factor. Boldface denotes significance at P ≤ .05.  P ! :002). This pattern also held when short-zone samples were excluded (F 1, 9 p 6:4, P p :032) and tall-zone samples were excluded (F 1, 9 p 10:6, P ! :010). There was no significant difference between southern and northern marshes in observed heterozygosity (H O p 0:313 and 0.302, respectively; F 1, 15 p 1:34, P p:265). At the six marshes in which tall and short zones were collected, paired t-tests detected no significant difference in either H E (t p 20:206; df p 5; P p :845) orH O (t p 20:023;df p 5;P p :983).
We detected repeated genetic differentiation between talland short-form phenotypes that is consistent with independent evolution in multiple analyses. First, marsh-specific PCAs detected significant differentiation (ANOSIM P ! :001) between tall-and short-form plants at all six marshes for which we had samples from both zones ( fig. 4C). Second, marsh-specific admixture analysis (k p 2; fig. 4C) also detected strong differentiation between tall-and shortform Spartina at all but one marsh, Saint Teresa in Florida, which is a microtidal system with less distinct zonation patterns in Spartina height (R. A. Zerebecki and A. R. Hughes, personal observation). These patterns are consistent whether including all loci ( fig. 4) or excluding a subset of loci (i.e., those with an F ST 1 90% of all loci; fig. A9; ANOSIM P ! :001). The PCA results are also robust when using a relatively small number of called genotypes (i.e., using STACKSgenerated genotypes; fig. A10) instead of genotype likelihoods (ANOSIM P ! :001 for all comparisons).
Third, an AMOVA confirmed significant hierarchical genetic structure among the six marshes from which we had both short-and tall-form samples, as well as between zones within these marshes (table 2). When all six marshes were included, 28% and 13% of the genetic variance was partitioned among marshes and between zones, respec-tively; both hierarchical levels were statistically significant. In contrast, an AMOVA with only the two South Carolina marshes (separated by !2 km) indicated no significant genetic variance between marshes (0.3%; P p:332) but strong genetic structure among zones within marshes (10.3%; P ! :001). An identical pattern was found when analyzing the two Massachusetts marshes separated by 1 km: betweenmarsh variance was negligible (1.0%; P p :331), but amongzone variation was significant (21.7%; P ! :001). The level of differentiation (as measured by genome-wide mean F ST ) between short-and tall-form plants separated by less than 200 m (i.e., within marshes; mean F ST 5 SE p 0:204 5 0:029) was significantly greater than mean F ST between marshes separated by between 650 m and 4.5 km (mean F ST 5SE p 0:12850:019) and comparable to differentiation among marshes separated by 5-15 km (mean F ST 5 SE p 0:209 5 0:006; F 2, 16 p 14:939, P ! :001; fig. A11).
Finally, the admixture analyses of all 11 marshes ( fig. A8) and the six marshes sampled at tall and short zones ( fig. A12) indicated that genetic clusters become increasingly differentiated by tall-versus short-form samples as the a priori number of genetic clusters included increased (i.e., moving from k p 5 to k p 16). These splits were especially evident at five marshes: three in the northeastern United States (MAWS vs. MAWT; MASS vs. MAST; RIT vs. RIS) and the two in the southeastern United States where the field experiment was conducted (SCBT vs. SCBS; SCFT vs. SCFS). If these ecotypes had evolved once and subsequently expanded into all other marshes, we would expect genetic clusters that were specific to only tall-or short-form plants to occur at multiple marshes. The lack of phenotypicspecific clusters across marshes suggests that this differentiation evolved independently. Another view of the marsh-specific genomic response to ecotypic divergence comes from an RDA, which detected no loci that consistently differed between ecotypes across three or more marshes. There were six loci that were shared at two marshes, but the identity of those two marshes differed by locus.

Discussion
Examples of parallel adaptive phenotypic divergence in independent populations across similar environments are particularly fascinating because they imply a repeatability in and predictability to evolution. The extent of parallelism can vary widely among replicated populations (Stuart et al. 2017) along a continuum from parallel (i.e., similar phenotypes) to divergent (i.e., discrete phenotypes; Bolnick et al. 2018;Thompson et al. 2019). However, this too may be predictable, depending on the selective environmental landscape and nonadaptive mechanisms including genetic drift, gene flow, mutation, and standing genetic variation (Bailey et al. 2017;Bolnick et al. 2018;Thompson et al. 2019). Our study provides compelling evidence that differences between talland short-form Spartina in part reflect parallel and independent microgeographic (!200 m) genetic divergence among habitat-associated genotypes and their traits in this wide-spread salt marsh foundation species. We detected a genetic basis for Spartina seedling height in a common greenhouse environment, and we confirmed adaptive divergence in fitness (survival and seed production) in a field reciprocal transplant experiment. We also documented signals of ecotypic divergence in two morphological traits-stem height and biomass allocation-that likely have adaptive value in this system. Shorter stems and reduced allocation to aboveground biomass of short-origin plants in the tall zone ( fig. 5) both result in a higher likelihood of complete submergence when inundated, as well as shading from neighboring plants. In turn, submergence and/or shading can limit biomass production (Smith 1982), an important component of long-term fitness in perennial plants (Younginger et al. 2017). Finally, we detected repeated genetic divergence between zones at multiple independent marshes across the geographic range of Spartina, indicating that our experimental results are likely applicable beyond the specific sites where the experiments were conducted. Our study highlights that species with welldocumented patterns of microgeographic phenotypic divergence that occur repeatedly across environmental gradients, like Spartina, are ideal systems to test for parallel evolution and its underlying ecological and genetic mechanisms.
The most parsimonious explanation for the maintenance of spatial differentiation in plant height and biomass allocation is that these traits, or traits correlated with them, Figure 5: Summary figure to illustrate the results of our reciprocal transplant experiment. Tall-origin Spartina in the tall zone (left) is 60% taller and has 58% greater above-to-belowground biomass than short-origin Spartina, illustrating the genetic constraint on short-origin plants. In contrast, tall-origin and short-origin Spartina in the short zone (right) have similar morphology. Adapted from Mudd and Fagherazzi (2016, fig 12.2) Image credit: M. Freedman. generate local fitness consequences (i.e., are locally adaptive) and that strong, local postsettlement selection removes the phenotype-environmentmismatch ineachgeneration (Schmidt and Rand 2001;Hereford 2009). The observed adaptive phenotypic divergence in Spartina (fig. 2) is likely also reinforced by mechanisms that minimize gene flow and generate the genomic differentiation revealed by our SNP survey. Seed dispersal barriers seem an unlikely mechanism, given the short distances (!200 m) relative to documented seed dispersaldistancesinSpartina(Tayloretal.2004) acrosswhich we detected both genetic divergence in traits and differentiation at SNPs. Barriers to pollen dispersal (e.g., differences in flowering time; Crosby et al. 2012;O'Connell et al. 2020) may act as a reproductive barrier to gene flow between zones, as in other systems (Hall and Willis 2006;Richardson et al. 2014). The relative influence of restricted gene flow and selection on the generation and maintenance of genetic and phenotypic differentiation remains to be tested in this system.
Parallel adaptive phenotypic divergence can arise through either shared genetic changes, alternative genomic routes, or some combination of both (Elmer and Meyer 2011;Soria-Carrasco et al. 2014). The extent of parallelism is typically highest for fitness, followed by phenotypes and then genes and mutations (Bailey et al. 2015;Thompson et al. 2019) due to many-to-one mapping (e.g., many distinct genes or phenotypes yield the same phenotype or function, respectively; Bolnick et al. 2018). Our study demonstrates that while phenotypic differentiation is parallel between marshes (i.e., tall-and short-form Spartina repeatedly evolve), the genomic responses underlying this differentiation do not appear to be parallel, according to multiple lines of evidence. First, differences in Spartina genetic structure and composition across biogeographic zones (figs. 4, A8) reflect thousands of generations of regional differentiation (Wares 2002;Blum et al. 2007) and likely generated regional differences in standing genetic variation on which selection has acted. Gene reuse in parallel phenotypic evolution is more likely in recently diverged populations (Conte et al. 2012), given the impact of historical legacy on standing genetic variation and resulting genetic constraints on responses to similar selection (Langerans and DeWitt 2004;Lucek et al. 2013;Thompson et al. 2019). Second, theory and empirical results predict that changes in polygenic traits (e.g., plant height) are likely mediated by different sets of loci between populations (Bolnick et al. 2018;Walsh and Lynch 2018). Third, a preliminary RDA yielded no anonymous markers that strongly and consistently correlated with ecotypic differentiation. Finally, differences among marshes in the particular abiotic variables underlying environmental differences across tidal elevation may yield marsh-specific genomic responses. For example, the soil salinity gradient across tidal elevation is typically much stronger in southern marshes than northern marshes (Pennings and Bertness 2001), and this spatial variation may alter both environmental selective pressures and adaptive responses across habitats. Cryptic environmental differences between seemingly equivalent habitat types are reflected in deviations from phenotypic parallelism in lake-stream population pairs of sticklebacks (Stuart et al. 2017), and interestingly, slight variation in the selective landscape may underpin the variability in differentiation among tall-and short-form Spartina populations (e.g., Atlantic vs. Gulf of Mexico marshes).
Our finding of microgeographic adaptation to tidal elevation in Spartina is likely conservative for several reasons. First, our experiment was not ideally designed to detect local adaptation in survivorship, as we used 1-year-old seedlings as our field transplants. Thus, we may have missed differences in survivorship that occurred earlier in the life cycle (e.g., seed dormancy, germination, or seedling establishment; Galloway and Fenster 2000;Bischoff et al. 2006). Mortality is known to be highest during early life history stages of Spartina (Metcalfe et al. 1986): 28% of seedlings died even in the optimal environmental conditions of our greenhouse. In addition, the low survivorship of seedlings (!30%) in the tall zone lowered our statistical power to detect origin-zone differences in reproduction in the tall zone, as did the short duration of our experiment relative to the life span of this perennial plant. Second, we did not observe any negative consequences for tall-origin Spartina in the short zone for individual traits (i.e., no antagonistic pleiotropy; fig. 3), consistent with prior results that the magnitude of local adaptation is greater when examining composite fitness than individual traits (Hereford 2009). Perhaps, individual trait trade-offs only become evident under other ecological contexts or environmental extremes (Bischoff et al. 2006;Latta et al. 2012). For instance, short-origin plants allocated a greater proportion of resources to belowground biomass and thus may have been superior competitors to tall-origin plants over longer time frames, particularly in the denser short zone (app. 3).
Our study confirms prior research (e.g., Shea et al. 1975;Valiela et al. 1978;Mendelssohn and McKee 1988;Morris et al. 2002) that local environmental conditions contribute to the observed trait variation in Spartina across tidal elevation. In particular, there was plasticity in height of talland short-origin seedlings when transplanted into the two zones ( fig. 3A). The magnitude of difference in plant height in our common garden experiment (~5 cm) was much smaller than that in the field (~60 cm; fig. 1B, 1C; app. 3), because tall-zone seedlings never grew as tall in the greenhouse as we recorded them growing among natural populations ( fig. 1B) or the field experiment ( fig. 3A). This pattern highlights the importance of environmental context to the tall phenotype. Our results reinforce the notion that genetic differentiation and phenotypic plasticity are not mutually exclusive (Richards et al. 2006;Ghalambor et al. 2007), consistent with variation in plasticity among plant populations to light availability (Winn and Evans 1991;Schmitt 1993), water availability (Heschel et al. 2004;Pratt and Mooney 2013), and habitat type (Gray 1985;Thompson et al. 1991).
Indeed, heritable variation in plasticity has been observed within Spartina (e.g., Elsey-Qurick et al. 2011;Bernik et al. 2018). So it should not be surprising that tall-and shortform Spartina can vary in their plastic response, and such heritable variation in plasticity may help resolve conflicting prior results. Nitrogen addition experiments have shown an increase in height in short-form Spartina; however, this increase in productivity was not equivalent to natural tall-form stands (Valiela et al. 1978;Howes et al. 1986). Although this difference may be due to the combination of other abiotic factors typical of the tall-zone environment (soil waterlogging, sediment oxidation, salinity, etc.; Linthurst and Seneca 1980;Howes et al. 1986;Mendelssohn and Morris 2002) not concurrently manipulated, an alternative hypothesis is that tall-origin seedlings may have either greater recognition of the environmental cues for the tall zone or an underlying genetic architecture that more strongly prompts this plastic response (DeWitt et al. 1998;Callaway et al. 2003). Understanding the relative importance of genetic differentiation and plastic responses and their potential interaction in generating and maintaining these two growth forms will be crucial to fully comprehending the ecological consequences of our observed microgeographic genetic divergence.
Foundation species drive ecological functions of entire ecosystems (Ellison 2019); thus, any adaptive phenotypic divergence within them can have far-reaching impacts on community and ecosystem processes (Norberg 2004;Whitham et al. 2006;Hughes et al. 2008;Matthews et al. 2011;Urban et al. 2020). In this system, Spartina structural traits such as stem height and biomass allocation are key drivers of marsh sediment accumulation rates (Mendelssohn and Morris 2002;Morris et al. 2002), and thus the genetic and environmentally mediated microgeographic divergence of these traits is likely to be critical for the response of marshes to sea level rise (Kirwan and Murray 2007). Our results suggest that with increasing inundation, the ability of tall-form plants to migrate landward may be crucial for stabilizing the retreating marsh edge and preserving marsh function ( fig. 5). We also add to the growing evidence of parallel adaptive phenotypic evolution occurring across a wide range of taxa (e.g., Rundle et al. 2000;Nosil et al. 2002;Roda et al. 2013) and highlight the importance of exploring the commonality of this process in foundation species. Greater understanding of the ecological factors driving and evolutionary processes underlying microgeographic adaptation in foundation species is needed, both to advance eco-evolutionary research and to inform the management and conservation of critical habitat-forming species and their associated ecosystem functions and services.