Populations within species often exhibit variation in traits that reflect local adaptation and further shape existing adaptive potential for species to respond to climate change. However, our mechanistic understanding of how the environment shapes trait variation remains poor. Here, we used common garden experiments to quantify thermal performance in eight populations of the marine snail Urosalpinx cinerea across thermal gradients on the Atlantic and the Pacific coasts of North America. We then evaluated the relationship between thermal performance and environmental metrics derived from time-series data. Our results reveal a novel pattern of “mixed” trait performance adaptation, where thermal optima was positively correlated with spawning temperature (cogradient variation), while maximum trait performance was negatively correlated with season length (countergradient variation). This counterintuitive pattern likely arises because of phenological shifts in the spawning season, whereby “cold” populations delay spawning until later in the year when temperatures are warmer compared to “warm” populations that spawn earlier in the year when temperatures are cooler. Our results show that variation in thermal performance can be shaped by multiple facets of the environment and are linked to organismal phenology and natural history. Understanding the impacts of climate change on organisms therefore requires the knowledge of how climate change will alter different aspects of the thermal environment.

(a) Natural History and Environmental Context

Urosalpinx is a predatory snail that is native from Florida to Nova Scotia and was introduced to the Pacific coast of North America in the late 1800s via American oyster (Crassostrea virginica) culture [34,35]. We quantified patterns of thermal performance from populations sampled across both the native and introduced coasts because they experience radically different thermal regimes. For example, while mean temperature and growing season length both decrease strongly from south to north along the Atlantic coast of North America, the gradient is much weaker and cooler on the Pacific coast [26,36]. Non-native species, such as Urosalpinx, provide an opportunity to compare intraspecific trait performance across different environmental gradients [37]. Although demographic history and founder effects have the potential to alter population responses to environmental regime and can confound interpretation of physiological trait data [38,39], the introduction of Urosalpinx to the west coast largely ended by the 1900s when transcontinental oyster imports ceased [40,41], allowing for 120 years of potential adaptation to the non-native climate regime. Further, in San Francisco Bay alone, 1.7 million kg of eastern oysters (Crassostrea virginica) were transplanted. Due to unregulated collection and transport processes at the time, large amounts of Urosalpinx were probably also introduced [41], which would greatly decrease the likelihood of a strong genetic bottleneck with drastic reduction of genetic diversity compared to source populations. Further, because this species undergoes direct development (i.e. there is no planktonic larval stage), dispersal and gene flow are likely limited among populations, suggesting a high potential for local adaptation [13,42].

(b) Broodstock Field Collection and Common Garden Experiment

We examined physiological performance of lab-reared offspring from broodstock mothers sourced from eight populations of Urosalpinx to evaluate the effects of environmental drivers on local adaptation. Experiments were conducted on juveniles that experienced controlled environmental conditions for their entire embryonic and juvenile life until cessation of experiments described below. To produce juvenile Urosalpinx, we collected broodstock adult Urosalpinx from eight sites, six from the Atlantic and two from the Pacific from March 15-June 9, 2019 (figure 2, table S1, [42]). 

At each site, we hand-collected at least 30 adults in the low intertidal to shallow sub-tidal zone from oyster reefs, pier pilings, and boulders across a sampling area of 300 m². We transported Atlantic-collected snails in aerated coolers of seawater from collection sites to the University of Massachusetts Amherst Gloucester Marine Station. Samples from Humboldt Bay and Willapa Bay were collected in a similar fashion except that they were overnight shipped in plastic bags with saltwater-moistened paper towels and immediately placed in a holding tank upon arrival. No mortalities occurred as a result of collection or shipping.

We maintained broodstock Urosalpinx in a recirculating seawater system at the Gloucester Marine Station at 12 °C (salinity 30 PSU) until needed for experimentation and as other populations were collected. Once all populations were established in the lab, we raised the water temperature to 20°C over the course of a week (1 °C/day) based on previous work indicating this was a suitable spawning temperature [2]. Populations were kept separate in plastic aquaria with aeration and fed blue mussels (Mytilus edulis), acorn barnacles (Semibalanus balanoides), and eastern oysters (Crassostrea virginica) ad libitum. Urosalpinx mothers began laying egg capsules on June 6^th 2019, with eggs being laid by 4-8 mothers per population. We collected egg capsules laid between June 6^th and July 4^th, 2019 and kept the egg capsules separated by population in labelled tea strainers (Tops Permabrew, Darien, CT). We maintained eggs at 20 °C and 30 PSU and checked strainers daily for hatchling emergence. We collected hatchlings for use in the common garden experiment within two days of hatching.

To maintain experimental temperature, hatchlings were enclosed within individual, labeled tea strainer halves that were floated within 30 L bins of 20 L aerated seawater (salinity = 30 PSU), which were in turn immersed in temperature-controlled seawater tables. These bins served as our temperature sub-replicates; we used three bins per experimental temperature. We placed nine snails from each of the eight tested populations in the six temperature treatments, distributing snails equally across the three replicate bins per temperature treatment (see figure S1). Each hatchling had a unique identification number to track date of entry into experiment and initial shell size. Hatchlings snails were supplied ad libitum with 3-5 mm oyster spat (Crassostrea virginica) per hatchling (Muscongus Bay Aquaculture, Bremen, Maine) for the duration of the common garden experiment. We checked tea strainers every 3-5 days (depending on temperature treatment) and all oysters were inspected using a stereomicroscope for signs of consumption, indicated by the presence of a drill hole (figure S2). During these checks, we replenished prey oysters with at least twice the number of observed consumed oysters to ensure ad libitum conditions were met.

For hatchlings, we digitally photographed each snail aperture-down using a stereomicroscope (Leica s9i Leica Microsystems GmBH, Wetzlar, Germany) and measured shell height from spire tip to distal siphonal canal using ImageJ (National Institutes of Health, Bethesda, Maryland) [3]. Digital measurement of hatchlings was used to assess initial height as opposed to using vernier calipers due to the small size of hatchlings (1.55 mm ± 0.202, mean ± SD) and the potential for calipers to damage shells. We recorded final snail measurements after the duration of the common garden experiment using vernier calipers, as snails were now large enough to be handled.

Details of Environmental Data Sources

We sourced environmental data from NOAA National Data Buoy Center (NBCC), the National Estuarine Research Reserve System (NERRS), the University of Virginia, the Pacific Shellfish Institute, and the Wiyot Tribe (table S1), based on 1) proximity to Urosalpinx collection site, 2) hydrology similarity between the data source and the collection site, and 3) completeness of annual data. To account for potential outlier years, we identified four years of complete annual water temperature data from each site between 2012 – 2019 to maintain data quality and completeness (table S1). We sourced data from all but one site (Skidaway, GA) within 15 km of the collection site. Skidaway’s closest environmental data source was only complete for one year (2018), and thus we used the next closest and most complete dataset 70 km away. The water temperature between the incomplete, proximal Georgia data source was highly correlated with the complete, distant Georgia data source (slope = 1.02, intercept = -0.501, R2 = 0.981).

To construct Urosalpinx thermal performance curves, we conducted a common garden experiment. We exposed hatchlings from the 8 populations to 6 chronic experimental temperatures (16, 20, 24, 26, 28, and 30°C) chosen to capture T_opt based on past experiments [43]. These temperatures are also realistic when compared to habitat temperature across populations (min-max: 2–37.5°C). Growth rate was measured using snail shell height, which is correlated with body mass [42]. First, juvenile snails were measured for an initial shell height and then randomly assigned into common garden temperature treatments. Snails were less than 24 hours of age (post-hatch) when they entered the common garden experiment that lasted for 24 days. On the last day, we measured shell height and calculated growth rate as the difference between initial and final size. We counted snails that died over the duration of the experiment to quantify survivorship in the common garden experiment, but these data points were excluded from growth analyses. See Supplementary Material S1 for further details about common garden experimental design and measurement of Urosalpinx growth rate.

(c) Environmental Predictors

In order to quantify environmental drivers of variation in growth rates in Urosalpinx, we derived nine metrics from four years of temperature data sourced from environmental loggers co-located within 14 km of the collection site for each population. One exception was the Georgia data logger, which was located 70 km away from the collection location but was highly correlated with local temperature logger that had a shorter record (Supplementary Material S1). We selected four years of temperature data from 2012-2019 based on the completeness of the record and to maximize temporal overlap among sites (table S1 , Supplementary Material S1). From this data, we calculated: 1) latitude, 2) summer mean temperature (June 1 – Sept 30), 3) upper 25^th temperature percentile of the summer period, 4) upper 10^th temperature percentile of the summer period, 5) maximum recorded temperature, 6) season length (number of days) where daily mean exceeded 10°C, 7) season length where daily mean exceeded 12.5°C, 8) the mean temperature for the first month of spawning, and finally 9) the mean temperature for the maximum period of spawning (table S2). We included length of season as a predictor because theory predicts organisms exposed to shorter growing seasons (i.e. high latitudes) are selected for faster growth [5,6,27]. We selected two likely lower temperature limits to calculate season length for Urosalpinx, 10°C and 12.5°C, based on reported absolute lower limit for feeding [44,45] and a breakpoint in oxygen consumption rates [46], respectively. We included mean temperature during spawning, as one of our hypotheses of T_opt behavior with environment is that high latitude populations experience warmer spawning periods than do low latitude populations [18]. We determined initial and maximum spawning periods as reported by [34] from the Atlantic and observations from the Pacific [47]; where no records of spawning periods could be found for a site, we used the closest neighbor site (table S4). We selected this broad range of variables as we did not necessarily know a priori which were relevant, and because the grouping of correlated predictors in AIC tables indicate what aspects of the environment best explain trait adaptation (i.e. season length generally, as opposed to season length above a specific threshold).

(d) Statistical Analysis

We used a two-step analysis framework to determine the environmental mechanisms driving growth rates in Urosalpinx populations in R [48]. First, we constructed and fit  nonlinear regression models to thermal performance curves with initial snail size as a random effect (contributing 2.4% variance, means ranging from 1.44 ± 0.180 to 1.69 ± 0.204 mm among populations) using the rTPC and nls.multstart packages for each population to quantify thermal performance curve attributes (T_opt and MTP, temperature at which maximum growth occurs and the maximum growth rate, respectively) for each population across the six common garden temperatures using the Rezende equation [49,50].  For each of the eight populations we fit three models based on the three replicate bins across the six common garden temperatures where populations were randomly assigned (figure S1). To produce 95% confidence intervals about each model prediction, we used non-parametric case resample bootstrapping on each population-bin model using rTPC and  (table S5, [49]). Once  models were fit to the data, we extracted the thermal optima (T_opt) and maximum trait performance (MTP) of each thermal performance curve (table S5, figure S3). We then modeled the T_opt and MTP for each population against a suite of environmental metric predictors (table S2) in a model-selection framework using generalized linear mixed models with gaussian error distribution and with population as a random effect using glmmTMB [51]. Each environmental predictor was used only once per model to identify which predictor best describes patterns in trait performance and to avoid  the issue of multicollinearity in models with multiple correlated predictors (i.e. where VIF > 4; table S3). We used Akaike’s Information Criterion (AICc) to select the greatest supported model, with a cutoff of ΔAICc < 2 [52]. For MTP, multiple predictors fell within the model selection criterion, and so we performed model averaging. To evaluate the possibility that analyses were influenced by invasive populations from the Pacific, we conducted a sensitivity analysis by constructing identical models with Pacific populations excluded. Our analyses were not sensitive to the removal of Pacific sites from analysis; both the best-supported environmental parameters and the significance level (p < 0.05) were maintained for the MTP and T_opt analyses. We therefore present the full analysis of Atlantic and Pacific sites in the results. For survival, we used type II ANOVA from the car package [53] on generalized linear models with a binomial error distribution and logit link to analyze if population, common garden temperature, or the interactive effects of population and temperature affected Urosalpinx survival in the common garden experiment.

The analyses for this manuscript were performed in R. To view our scripts and reproduce our analysis, you must download the latest version of R and RStudio.

Several files are included in this repository, described below.

1) README.md - This ReadMe document must be opened in R. Information and metadata from each script and data file are outlined here. 

2) Uro_ch_2.Rproj - R project file 

3) data.zip - Zip file of all data used in this analysis

4) scripts.zip - Zip file of all R scripts used in this analysis