A Pleistocene legacy structures variation in modern seagrass ecosystems.

Distribution of Earth's biomes is structured by the match between climate and plant traits, which in turn shape associated communities and ecosystem processes and services. However, that climate-trait match can be disrupted by historical events, with lasting ecosystem impacts. As Earth's environment changes faster than at any time in human history, critical questions are whether and how organismal traits and ecosystems can adjust to altered conditions. We quantified the relative importance of current environmental forcing versus evolutionary history in shaping the growth form (stature and biomass) and associated community of eelgrass (Zostera marina), a widespread foundation plant of marine ecosystems along Northern Hemisphere coastlines, which experienced major shifts in distribution and genetic composition during the Pleistocene. We found that eelgrass stature and biomass retain a legacy of the Pleistocene colonization of the Atlantic from the ancestral Pacific range and of more recent within-basin bottlenecks and genetic differentiation. This evolutionary legacy in turn influences the biomass of associated algae and invertebrates that fuel coastal food webs, with effects comparable to or stronger than effects of current environmental forcing. Such historical lags in phenotypic acclimatization may constrain ecosystem adjustments to rapid anthropogenic climate change, thus altering predictions about the future functioning of ecosystems.

The distribution and composition of Earth’s life are shaped by the environment across timescales both long and short. Over millennia to millions of years, climate and other environmental drivers produced characteristic vegetation types that underlie the Earth’s major biomes (1, 2). Although the character of biomes is shaped by the prevailing environmental conditions, legacies of past climatic and geological events also can strongly influence the composition and distribution of life (3, 4).

Such historical legacies have especially far reaching consequences when they act on foundation species, which create the habitat structure that defines biomes and underpins ecosystem functioning (5, 6). Foundation species often give their names to the ecosystems they create, such as redwood forests and coral reefs. Because many foundation species span broad geographic ranges, their populations are shaped by a mix of present and past conditions. Importantly, the traits of these key species—their morphology, phenology, and so on—are also geographically variable, and this variation can ripple through ecosystems to affect associated communities that drive fluxes of matter and energy (7).

Traits are shaped over time by environmental selection acting on genetic variation. At the population level, trait responses to environmental change depend on genetic composition and diversity. However, trait responses also depend on their degree of phenotypic plasticity: Some traits are plastic while others are fixed. Major historical events such as ice ages and sea-level changes represent an extreme form of environmental selection on species traits that has only been appreciated more recently (8). Now, with Earth’s climate changing faster than it has in at least 10,000 y (9), a central question is whether organisms and ecosystems can track changing conditions by acclimatizing and adapting, or are instead constrained by the legacies of such past events (10–12). This is particularly concerning for foundation species, which are challenged by rapid Anthropocene change in coastal areas worldwide, threatening the services they provide to nature and people.

Here, we explore how the growth form and associated community of the coastal foundation species eelgrass (Zostera marina) are shaped by current environmental conditions and a complex history of evolution and dispersal. Eelgrass inhabits shallow marine waters across the Northern Hemisphere, from warm temperate regions to the Arctic in both the Atlantic and Pacific oceans (13). Descended from one of the few flowering plant lineages (Alismatales) that colonized the sea, eelgrass is a quintessential foundation species, forming often monospecific stands that create habitat and food for diverse animals such as waterbirds (14), vulnerable megafauna (15), and invertebrate prey of commercially important fishes (16). Dense stands of eelgrass also accumulate organic matter in sediments and biomass, making them key sinks for natural sequestration of both nutrients and “blue carbon” (17).

Eelgrass has experienced an eventful history. Z. marina originated in the Pacific Ocean between 10 and 5 million years ago (mya) (18) and colonized the Atlantic via the trans-Arctic exchange sometime after ∼3.5 mya (19). In the Atlantic, its distribution fluctuated strongly with the many Pleistocene glacial–interglacial cycles of ice scour, varying temperatures, and changing sea level. The expansive range of eelgrass and its broad environmental tolerances make it an ideal candidate for addressing a general and timely question: How important is the long-term legacy of historical events such as ice ages in shaping ecosystems compared with current environmental conditions? To answer this question, we sampled living eelgrass and associated organisms at 50 sites across its range () and compared the contributions of the current environment versus the legacy of evolutionary history—recorded in neutral genetic markers—in explaining the growth form of eelgrass and the biomass of organisms that associate with it.

Results

Geography of Eelgrass Form and Genetics.

Our global analysis shows that eelgrass growth form and biomass vary markedly across its range and align with strong genetic divergence between the Atlantic and Pacific oceans (19), as well as within oceans, especially in the Pacific (Figs. 1 and 2 and ). Eelgrass form in turn predicts geographic variation in biomass of associated invertebrates and algae (periphyton; Fig. 3). Indeed, we find that global variation in the structure of this key coastal habitat owes as much to the legacy of the Pleistocene colonization of the Atlantic (4) and subsequent selection as it does to current environmental forcing across its expansive range (Fig. 4).

FIG. 1.

Variation in eelgrass growth form and genetic structure across the Northern Hemisphere. (A) Map of the 50 ZEN sites, with symbol color corresponding to ocean (Pacific, green; Atlantic, blue) and light green showing geographic distribution of eelgrass. (B) Variation in eelgrass growth form and biomass by site ( for site names). Numbers list mean canopy height (H), shoot density (D), and belowground biomass (B) for eelgrass populations at the extremes of the distributions. (C) Variation among sites in neutral genetic markers of evolutionary history, indexed by genetic axes FCA1 and FCA2 summarizing variation across 24 microsatellite loci. Eelgrass growth form and biomass represent the first two axes from a PCA of six eelgrass growth and morphological characteristics (). (C, Inset) A neighbor-joining tree of pairwise FST distances among all populations, with the size of the symbol proportional to the inverse value of PCz1; larger symbols denote longer, more forest-like canopies. In both B and C the percentage of variation explained by each axis is shown.

FIG. 2.

Genetic predictors of eelgrass growth form and biomass across Atlantic (blue) (A and B) and Pacific (green) sites (C and D). Plots of partial effects of genetic structure (FCA2) on eelgrass growth form (PCz1) and biomass (PCz2) show residual variation attributable to genetic FCA2 after controlling for all other influences (denoted by | Z) in the best model chosen by AICc (). Values of PCz2 are inverted such that biomass rises along the y axis. Regression lines and 95% CIs are shown for those predictors with P < 0.05.

FIG. 3.

Effects of eelgrass form and environmental predictors on periphyton (A, B, E, F) and invertebrate biomass (C, D, G, H) across Atlantic (blue) and Pacific (green) sites. Plots show partial effects of eelgrass form (PCz1) and biomass (PCz2) on periphyton and invertebrate biomass per bottom (core) area (g dry mass core−1) after controlling for all other influences (denoted by | Z) in the best model chosen by AICc (). Conventions are as in Fig. 2. See for full model results.

FIG. 4.

Influences of environment and evolutionary history on variation in growth form (A, B) and biomass (C, D) of eelgrass and associated organisms (E-H) across the north Atlantic and Pacific oceans. Bars show direct and indirect effects of environmental drivers and eelgrass genetic composition (history) as a weighted sum of the effect sizes (standardized partial regression coefficients) contributed by abiotic environment (PCe1, PCe2, PCe3) and evolutionary history (FCA1, FCA2); see for calculations. Coefficients are from the best model chosen by AICc for each (log10) response variable ().

The most pronounced distinction in eelgrass growth form is between the Atlantic and Pacific oceans and is captured by eelgrass form axis PCz1 (), which distinguishes tall sparse “forests” throughout much of the Pacific range versus short, dense “meadows” most common in the Atlantic (Fig. 1). While eelgrass populations commonly exceed 1 m in height along both the east and west Pacific coasts, this tall forest form was rare in the Atlantic (and its marginal seas) at the shallow depths we sampled (Fig. 1; only 1 of our 30 Atlantic sites was close to 1 m in height). This finding extends a trend found previously at a subset of our sites (20) to a global scale. The differentiation of populations spans a 24-fold range in canopy height and over two orders of magnitude in belowground biomass (Fig. 1). This divergence in growth form between the Atlantic and Pacific aligns with genetic differentiation based on 24 microsatellite loci (), with distinct population clusters in the two oceans defined by multilocus genotypes (genetic axis FCA1; Fig. 1), genetic distances (Fig. 1, Inset), and a bimodal genetic diversity spectrum that reveals strong barriers to gene flow and divergent phylogeographic signals between the oceans ().

Eelgrass form, biomass, and genetic composition also vary within oceans. Genetic variation in FCA2 is much greater and more spatially structured in the Pacific compared with the compact Atlantic cluster (Fig. 1). This reflects the longer history of eelgrass in the Pacific and the genetic bottlenecks associated with subsequent colonization of the Atlantic and Pleistocene glaciations that drove repeated range shifts, local extinctions, and recolonizations, homogenizing eelgrass genetic composition in the Atlantic (19). The independent histories of Pacific and Atlantic eelgrass are also supported by evidence that genetic structure in the Atlantic primarily reflects gene flow, caused by repeated mixing during glacial advances and retreats, whereas that in the Pacific mainly reflects accumulation of mutations in a more stable region (), with isolation by distance saturating beyond about a 3,000-km distance in the Pacific (). Genetic structure of eelgrass across its range thus retains strong imprints of events from many thousands of years ago.

Influences of History and Environment.

To partition this influence of history from that of present-day environmental drivers, we constructed and compared a set of candidate linear models that quantified their association with eelgrass form and its associated community ( and ). These analyses confirmed the divergence of eelgrass growth form between oceans and the strong association of eelgrass growth form and biomass with genetic composition (FCA2) (). Indeed, our estimates of genetic influence on eelgrass form are likely conservative because genetic axis FCA1 splits the sites distinctly between the Atlantic and Pacific (Fig. 1) so that its effects cannot be isolated from other unmeasured factors that differ between the oceans. Illustrating this, an alternative model that excluded ocean basin as a predictor yielded a highly significant association of eelgrass form (PCz1) with genetic FCA1 ().

Given their different histories, we proceeded with separate analyses of the Atlantic and Pacific sites, which provide largely independent estimates of environmental and genetic effects, and thus tested robustness of results from the global analysis. These analyses confirmed that variation in eelgrass form and biomass remain strongly associated with genetic population history within both oceans (Fig. 2 and ), despite the comparatively low genetic differentiation within the Atlantic (Fig. 1). In contrast to the strong influence of genetic structure in both oceans, environmental influences were evident only in the Atlantic, where eelgrass tended to be shorter (higher PCz1) at cool, high-latitude sites, mostly driven by short-statured stands in northern Europe, and biomass was greatest in productive estuaries (environment axis PCe3; ). In sum, our analysis indicates that the well-documented plasticity in eelgrass growth form (21, 22) is bounded and that the bounds in growth form are strongly influenced by genetics and differ between the oceans.

Historical Effects on the Ecosystem.

Eelgrass growth form and biomass, which we found are strongly shaped by evolutionary history (Figs. 1 and 2), in turn influence ecosystem processes and services (16), including via the organisms that live within its canopy. Microalgae and detritus (periphyton) on eelgrass leaves are the principal food for the abundant invertebrates that shelter among eelgrass leaves, which in turn are key forage for fishes (23). Periphyton was denser in the long, forest form of eelgrass (low PCz1) in both oceans, and was denser at sites with higher eelgrass biomass (low PCz2) in the Atlantic (Fig. 3 and ). Invertebrate biomass was higher in the meadow form of eelgrass (high PCz1) in the Atlantic, but was lower in the meadow form in the Pacific (Fig. 3 and ). Invertebrate biomass tended to increase with eelgrass biomass in both oceans, although this trend was only significant in the Atlantic (Fig. 3). Thus, the strong genetic effects of evolutionary history on eelgrass form appear to flow up to influence key components of the food web in this widespread habitat.

To compare the relative influence of environmental forcing vs. evolutionary history on the eelgrass ecosystem, we summed all effects (using standardized regression coefficients from the best model) of each of these types. The sum includes both direct effects, for example, of climate (PCe1) on invertebrate biomass, and indirect effects, for example, climate effects on eelgrass form that in turn affect invertebrate biomass ( and ). Evolutionary history and current environment had similarly strong influences on eelgrass ecosystems overall, but patterns differed between oceans (Fig. 4). For eelgrass form (PCz1), describing the continuum from forests to meadows, current environmental forcing and evolutionary history were comparably important predictors in both oceans (Fig. 4). Variation in eelgrass biomass (PCz2) was influenced more by evolutionary history than by environment in the Pacific, whereas these predictors had roughly equal effects in the Atlantic (Fig. 4). Periphyton mass also was more strongly associated with eelgrass evolutionary history in the Pacific than in the Atlantic (Fig. 4). Invertebrate biomass was equally associated with environment and genetic structure in the Atlantic, whereas environmental effects were greater in the Pacific (Fig. 4). The estimated “direct” effects of evolutionary history on periphyton and invertebrates are, by definition, unrelated to plant form (); these may reflect the parallel evolutionary legacies of Pleistocene events in the species pool of algae and invertebrates that associate with eelgrass. In summary, eelgrass evolutionary history is at least as strong a predictor of eelgrass form and community biomass as current environmental drivers.

Discussion

Our results reveal a lasting legacy of evolutionary history in the form and biomass of a key coastal foundation species, which in turn shapes ecosystem structure and processes via their effects on associated algae and animals. That legacy reflects the longer history of eelgrass in the Pacific, where populations are more genetically diverse and differentiated, than in the Atlantic (), where a colonization bottleneck and Pleistocene glaciations strongly reduced eelgrass genetic structure (19). We speculate that the meadow-type eelgrass dominant throughout the Atlantic is descended from short-statured Pacific genotypes (Fig. 1) that survived the trans-Arctic invasion and subsequent environmental selection through glacial–interglacial cycles. The reduced genetic diversity spectrum () and the broader range of environmental conditions among our Atlantic sites compared with the Pacific () likely explain the stronger environmental influence on eelgrass form and associated community in the Atlantic ().

Our results show that eelgrass evolutionary history and growth form affect important coastal food web components throughout the Northern Hemisphere. The relative contributions of drift and selection to eelgrass form also remain unknown, but the influence of seagrass form on ecosystem processes and services is broadly understood (16, 20, 24). Advances in genome sequence annotation for eelgrass (25, 26), together with expanded studies of local adaptation (27), offer promise for linking genes to traits to ecosystem processes in eelgrass. Similar links have been demonstrated in other marine foundation species (28), and in riparian forest trees, showing that local adaptation can result in genotypic differentiation that influences the entire ecosystem (7, 29).

Much attention has focused on how wild populations adapt evolutionarily to anthropogenic environmental change (30–32), but there has been less attention to the converse question: How might legacies of evolutionary events constrain species responses to new environmental conditions? Our results stress that the influence of evolutionary history on plant growth form can be surprisingly durable, sometimes overriding the influence of present environmental control, and that this legacy affects habitat quality (Fig. 2), faunal community assembly (33), and related ecosystem properties (Fig. 3). Genetic diversity within foundation species often strongly influences associated ecosystems (34), including eelgrass (35), where trait diversity underlies these effects (36). And eelgrass shows remarkable capacity for local adaptation to variation in conditions (21, 37–39), as do other foundation species (7). However, the ecosystem-level consequences of genetic variation have not previously been addressed at the global scale nor linked to biogeographic history as done here. If legacies of past events have indeed constrained phenotypic response to environmental conditions over large regions, as our results suggest, then they may similarly constrain responses of foundation species and their associated ecosystems to the rapid environmental change under way in the Anthropocene.

Methods

Sampling Design.

To explore variation in structure and forcing of eelgrass ecosystems, we sampled 50 sites spanning the geographic range of eelgrass, all chosen to be apparently healthy populations growing at less than a 2-m depth at low tide and in relatively protected waters. At each site, we collected data from 20, 1-m2 plots spaced roughly 2 m apart. The 50 sites were sampled by 25 partner groups in the Zostera Experimental Network (ZEN) () along the east and west coasts of the Atlantic (including its marginal seas) and Pacific oceans (Fig. 1 and ), with 3 to 9 areas along each coast, for a total of 1,000 plots (). In each plot, we sampled eelgrass above- and belowground biomass, shoot density, canopy height (longest leaf length per shoot), and length of the leaf sheath (a nongrowing structure that encloses the basal parts of the growing young leaves) using previously described methods (20, 40), and we collected tissue samples for genotyping 24 microsatellite loci (see next section). We also sampled mass of fouling material (periphyton), and the abundance and biomass of mobile, herbivorous, and detritivorous invertebrates (40). Full methods are available in the ZEN Handbook in .

Genetic Sampling and Analysis.

Sampling and molecular analysis.

We collected eelgrass leaf samples from each of the 20 plots at each of the 50 sites as described above (). Leaf tissue was stored in silica gel within hours after collection. DNA was extracted from ∼20 mg of silica gel–dried tissue in 96-well plates using a silica-based cetyl trimethylammonium bromide (CTAB) protocol (41), except that samples were incubated in CTAB for 1 h at 60 °C. The 24 microsatellite loci, primer sequences, and multiplex combinations are given in . PCRs were performed in 96-well microtiter plates using the Qiagen Type-it Kit in a 6.2-µL reaction volume following the manufacturer’s instructions. The reaction profile consisted of 95 °C for 5 min followed by 30 cycles of 95 °C for 30 s, 56 °C for 1 min 30 s, and 72 °C for 30 s, with a final extension step of 60 °C for 30 min. PCR products were diluted 1:100 (apart from the “4-plex,” which was used undiluted) and fragment analysis was performed on an Applied Biosystems 3730 DNA analyzer with a 350 ROX internal size standard added to each well. Fragments were scored automatically using GeneMapper (Life Technologies) and rechecked by eye for each individual and locus. Samples with ambiguous or rare alleles were reamplified and regenotyped for confirmation. We tested for null alleles with MicroDrop (42) (10,000 permutations and 100 replicates), which does not rely on Hardy–Weinberg equilibrium assumptions to calculate null allele frequencies.

Data properties and basic metrics.

We identified clonal replicates based on multilocus genotypes (MLGs) using the software genclone 2.0 (43) and estimated genotypic richness for each population as described in Dorken and Eckert (44): R = (G − 1)/(N − 1). We excluded replicates of a genotype/clone for all downstream calculations unless described otherwise. Basic population genetic metrics and factorial correspondence analysis (FCA) were calculated in Genetix 4.05 (45) (). FCA is similar to principal-component analysis (PCA) but more suited to multilocus microsatellite data in that it can link clusters of individuals with subsets of variables (microsatellite loci in our case) to which they are significantly related (45).

Genetic diversity spectra.

We chose genetic diversity spectra (GDSs) for the differentiation analyses because it has shown good power to retrace large-scale biogeographical events that have shaped the distribution of polymorphisms (46), which is not possible with other measures such as FST. GDS is an individual-based method that summarizes the scale-dependent structure of genetic diversity and provides inferences about the relative role of processes (e.g., gene flow and mutation) that shape genetic differentiation among and within populations (47) (). We calculated genetic similarity among individuals using weighted Rozenfeld distance (RD) (47) and visualized the resulting GDS (the frequency distribution of all pairwise interindividual genetic distances) with RClone (v. 1.53.3) in R 3.2.2 (48). RD makes use of the information contained in allele length by calculating distance among individuals based on the difference in allele sizes assuming a stepwise mutation model (SMM). We next compared the GDS of the observed dataset against the distribution of two randomized datasets. The first randomized dataset tested the null hypothesis for panmixia (). We used a matrix of the observed unique multilocus genotypes and shuffled alleles randomly with 1,000 iterations using the R package Picante (49). After shuffling and in order to take into account the occurrence of clonality, we added back random genotypes from each simulated population to the “random dataset.” Each random genotype added back into the random dataset corresponded to the observed number of replicate genotypes in the original dataset. In this way, unique and clonal allele frequencies were accounted for. The second randomized dataset tested the null hypothesis for “lack of a phylogeographic footprint” () in the spatial distribution of allelic divergence. We used a matrix with all collected individuals (including replicate genotypes) and recoded each allele size to another randomly chosen size that was detected for that locus in the dataset. In this way, all informative values of allele size were removed, but genetic and genotypic differentiation, as well as allele frequencies and clonality, were maintained (46). The two null hypotheses were tested on the global dataset, the Pacific sites alone, and the Atlantic sites alone. The significance of the deviations between the GDS of the original data and the randomized datasets was tested by Kolmogorov–Smirnov tests in R for the GDSs based on RD.

Simulations of mutation and migration.

To investigate the relative importance of migration vs. mutation (), we used EASYPOP v. 1.7 (50) to model the GDS assuming different rates of migration and mutation, and compared modeled with observed GDS. The following parameters were used to model MLGs after 1,000 generations in EASYPOP: diploidy (but with each individual representing both sexes), and nonrandom mating with a proportion of 20% clonal reproduction (reflecting global results of measured genotypic richness). The proportion of selfing was set to 1%. Estimates for Z. marina outcrossing rates range from 0.61 to 1 (51); our chosen estimate was therefore within the range of actual observations. Since the goal of the analysis was to compare the importance of migration vs. mutation, and we kept clonality and selfing constant for this analysis, the actual values of these parameters are not important in our case, so long as they were kept constant. The simulated dataset consisted of 50 populations with 20 individuals each and 24 loci (representative of our global dataset). Migration and mutation rates were consecutively changed, and four examples of good and bad fit are shown in . A one-dimensional stepping-stone migration model was chosen, as it represents the simplest model of expansion. The 24 loci were allowed to recombine freely and were all assumed to have the same mutation scheme according to the single-step SMM. The number of possible allelic states was set to 99, assuming minimal variability of the initial population. Each simulation was run with 10 replicates.

Isolation by distance with binned GDS.

We visualized the frequency of interindividual RD distances vs. geographic distance bins and then tested for isolation by distance (IBD) using Pearson moment correlation in R. We used the mean of the genetic distance vs. the maximum of each geographic distance class (). Geographic distance was measured as “oceanographic distance” without crossing land in Google maps. Scripts used for analysis of GDSs and IBD are available at https://zenodo.org/record/3660013#.YWF8C0bMLpJ.

Neighbor-joining tree.

A genepop-formatted matrix of multilocus genotypes was imported into R (52) with adegenet (53), and pairwise Weir and Cockerham (55) PhiST values were estimated with hierfstat. A neighbor-joining tree was generated with ape (57).

Data Preparation and Assembly for Integrated Linear Modeling.

We modeled eelgrass form and ecosystem components as a function of seven environmental variables and multilocus genetic variation derived from 24 microsatellite loci (). Environmental variables included sea surface temperature (annual mean and range), photosynthetically active radiation, dissolved phosphate concentration, and sea surface chlorophyll concentration (all as annual means, except for annual temperature range) from the Bio-ORACLE database (52); annual mean cloud cover was from the WorldClim database (54), salinity was measured locally during the experiment, and the nitrogen content of eelgrass leaves was measured at each site as a proxy for local nutrient availability.

Eelgrass form and biomass: PCA.

We summarized eelgrass characteristics in a PCA that incorporated log10-transformed values of above- and belowground biomass, canopy height (longest leaf length), shoot density, and length and width of the leaf sheath (). The first PC axis, hereafter eelgrass form (PCz1), explained 55% of the variation and indexed the continuum from “forests” with tall canopy and sparse shoots (low eelgrass form PCz1) to “meadows” with short canopies and dense shoots (high PCz1) (). The second axis, hereafter eelgrass biomass (PCz2), explained 27% of the variation. High values of PCz2 indicate lower biomass per unit area, especially belowground, and sparser shoots within a given growth form; therefore, we reverse the direction of this axis in the figures to be more intuitive.

Environmental predictors: PCA.

Comparing effects of current conditions and historical evolutionary legacies requires first quantifying them in comparable units. We approached this by condensing effects of environment and genetic composition into separate sets of component axes that can be used as summary indicators, expressed as standardized effect sizes. Ordination approaches such as PCA and factorial correspondence analysis (see below) have two advantages for analyzing our data. First, they reduce the dimensionality of the dataset by reducing the number of predictor variables, releasing the degrees of freedom, and increasing the power of the analysis. Second, they search for the best orthogonal axes explaining variation in the component predictor variables, with the result that they minimize collinearity among the predictors, which can cause problems in model fitting and parameter estimation. We condensed variation in the seven environmental variables via a PCA (). The first three environmental PC axes (PCe1 to 3) captured 78% of the variation in input variables. Generally, positive scores on environmental PCe1 indicate warm, bright conditions characteristic of low latitudes; high values of PCe2 indicate nutrient-rich conditions (primarily dissolved PO4 and eelgrass leaf %N); and PCe3 indicates productive estuarine conditions, with low salinity and high phytoplankton (surface chlorophyll).

Genetic predictors: FCA.

We used FCA to reduce the dimensionality of our 24-locus dataset and to capture differences in genetic distinctiveness among populations as a measure of evolutionary history. FCA is a multivariate ordination technique similar to PCA for use with categorical variables (45), in our case multilocus genotypes. The first two FCA axes (Fig. 1) account for only 24% of the variation, as expected from the high dimensionality of the genetic dataset. The first axis (FCA1) largely separates ocean basins (Atlantic and Pacific). The second axis (FCA2) mainly separates populations within the Pacific. The third FCA axis explained less than 8% of the variation, so we limited the analysis to the first two genetic axes to limit the total number of variables and streamline the models.

Linear Modeling Approach.

We sampled existing, apparently healthy eelgrass populations, so our analysis focuses on the controls on eelgrass growth form and associated organisms where it occurs, namely conditional on eelgrass presence. Our analysis does not speak to the factors determining eelgrass presence versus absence, which is well-known to depend on water quality and to some degree on top-down control (14, 58–60). However, since we sampled nearly the entire geographic and environmental range of the species, including marginal populations in the Baltic, Mediterranean, Adriatic, and White seas, we believe our conclusions are robust for the factors explaining global variation in growth form among eelgrass populations and the biomass of associated organisms.

Our analysis had two related goals: to identify the best set of environmental and genetic predictors for each eelgrass ecosystem response variable, and to obtain the best estimate of effect strength (partial regression coefficient) for each predictor in the final models. We approached both goals using an information-theoretic approach (61) by fitting and comparing a set of general linear models for each response variable, and comparing them using Akaike’s information criterion corrected for small sample size (AICc). Response variables were eelgrass growth form (PCz1), eelgrass biomass (PCz2), periphyton mass (mass core area−1), and invertebrate biomass (mass core area−1). Models of periphyton mass omitted data from one site, Sweden A, where periphyton was not measured. As environmental predictors we used the principal-component axes PCe1 to PCe3, and for evolutionary history (genetic composition) we used FCA1 and FCA2. Prior research shows substantial genetic and environmental differentiation between eelgrass communities in the Atlantic and Pacific oceans (19, 40), so ocean was also included as a categorical predictor with levels Atlantic and Pacific. With the exception of ocean, all other response and predictor variables were continuous. All variables were range-standardized prior to analyses (62) so that effect sizes could be compared in comparable units.

We used site-level mean values as inputs in the models (n = 50 sites in global analyses). The main reason was that data for most environmental predictors (with the exception of leaf nitrogen content measured during the study) were only available at a coarse resolution that corresponded to site scale. Using site-level data for all predictors and response variables ensured that degrees of freedom and spatial resolution were comparable for all variables and models, and also simplified interpretation, as we were able to use standard linear models rather than hierarchical mixed models.

Building a candidate set of models and comparing them can be approached in a variety of ways, and involves both philosophical and practical considerations. We used the following approach. First, a central goal of the analysis was to compare the influences of environmental and genetic predictors, so we opted to retain main effects of the five predictor variables (PCe1 to 3, FCA1, FCA2) in all candidate models. This ensured that we captured estimates of the effects of each variable, whether or not they were considered significant. Second, because we used site means as input, the total number of observations (50 site means) was low relative to the number of predictor variables and possible interactions, so we restricted our set of candidate models to include only main effects and two-way interactions. For each response variable, the candidate model set included main effects alone and main effects plus each of the possible two-way interactions among predictors (). The candidate models therefore differed only in which interaction terms were included. After fitting each candidate model set, we identified those models within 2 AICc units of the best one (61) (lowest AICc score). Where more than one model was within 2.0 units from the top score, we compared them manually to check whether differences in statistical significance and/or coefficient estimates of predictors between the top models changed interpretation of the results (see below).

Comparing Effects of Evolutionary History and Environment.

Finally, we compared the total estimated effects of environment vs. evolutionary history on each eelgrass ecosystem component. To do so, we first traced all possible direct and indirect causal paths from a given driver (e.g., FC1) to an eelgrass ecosystem component (e.g., invertebrate biomass; ). Direct effects were estimated as partial regression coefficients from the best model (). Indirect effects were visualized as chains of causal paths from the predictor (e.g., FC1) through intermediate variables (e.g., eelgrass form, PCz1) to the response variable (e.g., invertebrate biomass), and were estimated by multiplying the partial regression coefficients associated with each path in the chain. We used the absolute values of effect sizes to estimate relative importance, recognizing that effects can be important whether they are negative or positive. Finally, to estimate the total effect of, for example, evolutionary history on invertebrate biomass, we took a weighted sum of the paths from all predictors of that class (for evolutionary history, these are FC1 and FC2); that is, each estimated effect size was divided by its SE before summing. Weighting the components before summing assigned importance to each effect in proportion to the confidence we have in its estimate.

Tests of Robustness.

We conducted three additional analyses to test the robustness of our results and conclusions. First, the environmental data we used are mostly interpolated from the Bio-ORACLE dataset (52) and therefore coarse in geographic scale and often are based primarily on measurements offshore of the shallow estuaries we sampled. As a partial check on these limitations, we fit an alternate version of the main-effects model for each response variable but substituted water temperature, salinity, and eelgrass leaf nitrogen content measured locally, as well as estimated day length, in place of the regional Bio-ORACLE environmental variables (). The models using locally measured environmental predictors retained strong effects of evolutionary history (FCA1 and FCA2) on eelgrass growth form and also showed a positive association of invertebrate biomass with eelgrass leaf nitrogen content, echoing the similar finding of a positive effect of productive estuarine conditions (PCe3) on invertebrates using the Bio-ORACLE data ().

The second test addressed the challenge that the first axis of genetic variation, FCA1, was nonoverlapping between ocean basins (Fig. 1), and therefore partially confounded with the categorical variable ocean. We assessed the effects of this collinearity by fitting and comparing three models of the global dataset: 1) a full model including main effects of all environmental and genetic variables (PCe1 to 3, FCA1, FCA2) plus ocean as predictors; 2) the full model excluding ocean; and 3) the full model excluding FCA1. We then used AICc to compare among the three models. Models excluding ocean versus FCA1 generally had similar AICc values, meaning that the categorical predictor ocean and genetic FCA1 provided roughly equivalent but overlapping predictive information, and switching between them had no appreciable effect on coefficients or P values for any response variable ().

Third, we used random forest analysis to explore whether two local-scale predictor variables, not available in regional-scale environmental databases, might explain variation in eelgrass canopy height across the ZEN sites. Eelgrass height (leaf length) may be influenced by exposure to wave energy and by ambient light levels. In the absence of local wind speed and direction data from all sites, we estimated wave exposure as fetch across the water adjacent to the site by measuring the straight distance (in km) from the center of each of the 50 sampling sites to the nearest land in 8 directions (every 45°, setting the maximum distance to 20 km), and then averaging across the 8 directions. As an indicator of light limitation, we measured the stoichiometric ratio of C:N in leaf tissue (63). Aliquots (1 to 10 mg) of dried leaf tissue were analyzed for total organic carbon and total nitrogen using a Thermo Flash EA Series 1112 NC soil analyzer. Random forest analysis showed that both genetic axes FC2 and FC1 were stronger predictors of canopy height than either fetch or leaf C:N ratio ().