Use of settlement patterns and geochemical tagging to test population connectivity of eastern oysters Crassostrea virginica

Freshwater-dominated estuaries experience large fluctuations in their physical and chemical environments which may influence larval dispersal, settlement, and connectivity of populations with pelagic larval stages. We used settlement patterns and natural tagging along with numerical hydrodynamic model results to assess settlement and connectivity among oysters across the freshwater-dominated Mobile Bay−eastern Mississippi Sound (MB−EMS) system. Specifically, we (1) tested how freshwater inputs and associated environmental attributes in fluenced settlement patterns during high and low discharge conditions in 2014 and 2016, respectively, and (2) analyzed trace element (TE) ratios incorporated into multiple shell types (larval and settled shell of spat and adult shells) to determine if shells collected in situ incorporate temporally stable site-specific signatures. We also assessed if TE ratios compared between adult (TE natal signature proxy) and larval shells could infer connectivity. Larval settlement was 4× higher during low discharge than during high discharge when oyster larvae only settled in higher salinity regions (EMS). Spat and adult shells incorporated site-specific TE ratios that varied from weeks to months. Connectivity results (May−June 2016 only) suggest that EMS is an important larval source to EMS and lower MB. While we were able to infer probable connectivity patterns using adult and larval shells, more study is needed to assess the utility of adult shells as proxies for natallocation TE signatures. Results provide a baseline for measuring future larval connectivity and adult distribution changes in the MB−EMS system. Biological and geochemical data demonstrate the potential to identify environmental attributes that spatiotemporally mediate settlement and connectivity in dynamic systems.


INTRODUCTION
Many marine species have pelagic larval stages, and knowledge of larval dispersal and subsequent population connectivity is imperative to better understand population dynamics and manage populations (Thorrold et al. 2007). Larval dispersal and transport are difficult to quantify due to (1) uncertainty in lar-val sources and (2) difficulty of tracking microscopic larval movements in situ due to advection and high larval mortality (Levin 2006). Many species, including bivalves, also inhabit dynamic habitats often dominated by freshwater inputs (Pollack et al. 2011), which can modify larval transport and delivery of new recruits for settlement (Dong et al. 2012, Kim et al. 2013. Freshwater discharge, along with other environmental factors such as winds, tides, and surface heat exchange, alters the chemical and physical environment (e.g. salinity, water temperature, and circulation patterns) that mediates reproduction, survival, growth, distribution, and settlement processes, which ultimately mediate connectivity (Wilber 1992, Kim et al. 2010. While it is relatively well-established that fresh water inputs affect the distribution and abundance of bivalves, with lower flows into high freshwater input systems resulting in more abundant adult bivalve populations and vice versa (Wilber 1992, Soniat et al. 2012, specifically how changes in freshwater discharge influence dispersal of larvae, settlement of recruits, and ultimately population connectivity has been less studied. Geochemical tracers have been used as natural tags to indirectly track larvae, determine larval sources, define transport, and infer population connectivity (defined as where larvae originate from; Levin 2006) (Thorrold et al. 2002, Becker et al. 2005, Bradbury et al. 2011. Geochemical tracers such as trace elements (TEs) are incorporated into larval shells from the natal environment and retained throughout early development and settlement (Levin et al. 1993, Anastasia et al. 1998, DiBacco & Levin 2000. If the ratio of these elements is location-specific, then larval origins can be identified by comparing larvae of unknown origin to natal reference signatures. Reference TE signatures need to be spatially distinct, encompass the potential dispersal distance (Miller et al. 2013a), and have low temporal variation during the pelagic period (Campa na et al. 2000). Typically, natal reference signatures are determined using larval outplant studies, where hatchery-reared larvae are transplanted to larval homes at specific locations to measure TE ratio incorporation (Becker et al. 2007, Kroll et al. 2018. Outplant methods can be costly and time ly to do, and thus, identifying proxies, such as juvenile or adult shells collected in situ to represent natal reference signatures, can be potentially useful (Sorte et al. 2013). Difficulties arise when using different shell types to classify larval origins. For example, mineralogical differences between adult and larval shell are present; adult bivalves are sedentary and can integrate TE signatures over longer time periods (i.e. weeks to years) than larvae, reflecting TE ratios of multiple larval cohorts or time periods where larvae were not present in the water column. Therefore, site-specific TE ratios in adult bivalve shell can potentially serve as a time-integrated proxy for natalsite reference TE signatures, but more study is warranted prior to the use of adult shells.
To determine the utility of geochemical tracers to track larval origins and movements within a system, the spatial and temporal variation of site-specific TE ratios needs to be defined. Elemental spatiotemporal variability is determined by differences in environmental (physical and chemical) and biological factors (Lorens & Bender 1980, Vander Putten et al. 2000. Estuaries are ideal areas to study variation in TE ratios because they have high spatiotemporal environmental variability due to different geomorphologies, pollution sources, and freshwater inputs (Swearer et al. 2003, Thorrold et al. 2007). Previous estuarinebased studies have found spatial variation in TE ratios in bivalve shells across 1−10s of km (Norrie et al. 2016) and temporal variation on a scale of weeks to months in bivalve shells (Becker et al. 2005) and years in gastropod shells (Zacherl 2005). These findings suggest that spatiotemporal variation within estuaries should be adequate to distinguish larvae among natal sites using TE ratios. The ability to define spatial and temporal variation in TE ratios is of particular concern in freshwater-dominated systems where large-scale mixing may affect larval dispersal processes and reduce the spatial distinction of TE ratios among sites through time.
Oysters are commercially, ecologically, and biogeochemically important, particularly in the north-central Gulf of Mexico (nGoM), which is home to a few of the remaining commercially harvestable oyster populations worldwide (Beck et al. 2011). Eastern oysters are a spawning species that produce planktonic larvae lasting ~2−3 wk depending on environmental conditions (Medcof 1939). In the nGoM, eastern oysters are capable of spawning multiple larval cohorts throughout the year, typically between April and October (Ingle 1951, Butler 1965, Hayes & Menzel 1981, Kim et al. 2010. Spawning is temperaturedependent and is induced when temperatures reach > 25°C or following a rapid increase or decrease in temperature (Hayes & Menzel 1981, Saoud et al. 2000. Timing of spawning can influence connectivity because it regulates if gametes are re leased into optimal conditions for larval development and survival, which ultimately determines the larval pool for connectivity (Prytherch 1929). Larvae can disperse 0.1−100s of km (North et al. 2008, Haase et al. 2012, during which a larval aragonite shell develops (Stenzel 1963(Stenzel , 1964. Freshwater discharge can influence connectivity by affecting the circulation pattern as well as the growth and survival of early life stages, which are salinitydependent, with optimum salinities for larval development and spat settlement between 7.5−22.5 and 10−30, respectively (Calabrese & Davis 1970, Chatry et al. 1983. Following the pelagic larval period, oysters settle under suitable conditions and remain in the settlement location throughout juvenile and adult life stages, incorporating TE ratios of settlement locations into a calcite shell (Stenzel 1963(Stenzel , 1964 representing long-term temporal scales. Thus, oysters are promising model organisms to understand how changes in freshwater discharge and associated environmental attributes mediate larval settlement and connectivity on different spatial and temporal scales. To determine how changes in freshwater discharge influence oyster populations, we conducted a 2 part study to define larval oyster settlement patterns and connectivity using geochemical tags. We used the Mobile Bay−eastern Mississippi Sound (MB−EMS) system, a nGoM estuarine system with the 6 th largest freshwater drainage basin in the USA (Isphording & Flowers 1987), as a model system. Previous studies have found oyster settlement to be highest in the southwest side of MB−EMS, following an increasing salinity gradient from northern to southern MB and westward into EMS, with negligible settlement in the middle, east, and upper portion of MB (Hoese et al. 1972, Lee 1979, Saoud et al. 2000, Kim et al. 2010). In the first part of this study, oyster spat settlement data were collected from areas of known oyster settlement in 2 different years with low and high discharge conditions to examine how the flow patterns and asso ciated environmental attributes (salinity and temperature) determined by freshwater discharge, winds, and tides affect settlement of eastern oysters in the MB−EMS system. In the second part of the study, we (1) used newly settled oysters collected in the first part of the study and adult oysters to determine if larval, settled, and adult shells have sitespecific TE ra tios and (2) conducted an exploratory analysis to as sess connectivity using adult shells as natal reference TE signatures, comparing the larval component of shell from newly settled spat to adult shells to estimate larval origins. Because connectivity studies using TE ratios have not been conducted in freshwater-dominated systems, we verified our predictive me thods by testing if (1) larval and post-settlement components of spat shell were sufficiently different to distinguish between natal origins and settlement locations, (2) larval, settled, and adult shells differed among sites to determine spatial distinctions due to variability in the physiochemical environment, and (3) TE ratios in larval, settled, and adult shells varied temporally on seasonal or annual scales. These data may have important implications for understanding how the spa tial extent of oyster populations shift under changing discharge regimes in freshwater-dominated systems.

Sampling scheme
To determine spat settlement patterns, settlement plates were deployed at settlement sites (Table S1 in the Supplement at www. int-res. com/ articles/ suppl/ m673 p085 _ supp. pdf, Fig. 1) in the freshwater-dominated MB−EMS system throughout a region of historically known oyster habitat (Hoese et al. 1972, Lee 1979, Saoud et al. 2000, with sites EMS1−3 located 87 Fig. 1. Sampling locations included in the Mobile Bay−eastern Mississippi Sound system, including spat settlement plate and adult oyster collection sites. Stippled areas indicate known historical oyster reef locations (data from May 1971, Tatum et al. 1995, Alabama Department of Conservation and Natural Resources 2001 west to east in EMS and MB1, MB4−7 located south to north along a salinity gradient in MB. In 2016, 3 new settlement plate sites were added to include a natural reef site (EMS3) and historically productive oyster areas that have been unproductive in recent years (MB5 and MB7) (Stout et al. 1998). Plates were changed bimonthly from May to mid-September of 2014 and 2016.

Settlement plate moorings
Settlement plate moorings were attached to existing pilings, below the water surface to reduce potential vandalism. Moorings consisted of an 18.9 l bucket with a steel pole (25.4 mm × 2.1 m) cemented in the center with 18 kg of concrete, ratchet strapped to a piling. A PVC sheath holding 4 settlement plates was bolted in place 1 m above the bottom and accessible via snorkeling to allow easy access to plates without removing the mooring. Settlement plates (17 × 17 cm) made of HardieBacker ® Cement Board were deployed at a 45° angle (Carleton & Sammarco 1987, Dayton et al. 1989). Plates were pre-soaked in 0.7 μm filtered seawater for ~24 h to promote biofilm development in situ to increase settlement (Quayle & Newkirk 1989, Tamburri et al. 2008). Prefiltering of seawater was done to remove organisms that could potentially affect settlement. To reduce predation while allowing adequate water flow, a half-caging method was used by placing 3.5 mm plastic mesh with 2 sides open with the top and bottom closed (Kim et al. 2010).

Spat abundance
Two settlement plates were analyzed per site; only the inner 16 × 16 cm of the plate was analyzed (outer 1 cm excluded) to eliminate handling effects. Spat on settlement plates were examined with a Fisher Stereomaster zoom dissecting microscope (10−60× magnification). All spat were counted and distinguished as live (non-gaping and/or tissue present) or dead (spat scars or gaping oysters with no visible tissue).

Environmental attributes
To determine freshwater discharge, daily discharge data in 2014 and 2016 were downloaded from the US Geological Survey gauging stations 02428400 (Alabama River at Claiborne near Monroe) and 02469761 (Tombigbee River at Coffeeville); their sum was used as the total discharge into the MB−EMS system (Park et al. 2007). To determine environmental attributes at the time of sampling, temperature, salinity, and dissolved oxygen were measured at the depth of plates when changing plates bimonthly using a YSI pro 2030 handheld data sonde. To determine wind (direction and velocity) and water temperature, hourly data were downloaded from the Dauphin Island station DPIA1 (NOAA Tides and Currents). A wellvalidated 3-dimensional hydrodynamic model was used to simulate salinity at finer spatial and temporal resolution throughout the study area in 2014 and 2016. Detailed information of the model configuration and its application and validation can be found in Kim & Park (2012) and Park et al. (2014).

Settlement sampling statistical analyses
To determine spatial variability in spat settlement, regression analysis was used to compare average numbers of spat settled to distance from freshwater input (i.e. head of Mobile Bay), to compare to previous studies done in the MB−EMS system. To compare settlement patterns found in this study to those from previous studies conducted in the MB−EMS system, literature data were collected from Hoese et al. (1972), Lee (1979), Saoud et al. (2000), and Kim et al. (2010), and their average settlement trends versus distance from freshwater input were compared to the results from this study. Data collected from Lee (1979) and Saoud et al. (2000) were not included in regression analyses due to having too few sites where spat settlement was measured. To determine spatial and temporal variability in spat settlement, general linear models (GLMs) were used with settlement site and week as factors. Negative binomial GLMs with log-links were run for 2014 spat settlement. To account for over-dispersed count data with zero inflation, 2016 spat settlement data were analyzed using a zero-altered negative binomial GLM with a log-link using the 'pscl' package in R (Jackman 2008).
To determine if discharge differed between years, a generalized additive model was run with a reduced maximum likelihood smoothing function using the 'mgcv' package in R (Wood 2012). To determine spatial and temporal variability in salinity and temperature, GLMs were run with settlement site and week as factors for sites and time periods that had appreciable settlement. For comparisons between years, GLMs with year, week, and settlement site as factors were run with salinity data for the same sites and time periods, and a t-test was performed for temperature comparisons. Because environmental variables measured at the time of plate collection represented a single time point that is not necessarily representative of conditions in the MB−EMS system during settlement in the days prior to collection, settlement data were not directly compared to environmental data. Rather, temperature and salinity patterns throughout the settlement season were analyzed, and time points at the beginning of consistent settlement were considered 'settlement events.' To determine environmental conditions potentially indicative of spawn ing or more suitable conditions for settlement, temperature and salinity changes 2 and 4 wk prior to the settlement event were assessed. These times were chosen based on an average 2−3 wk pelagic larval duration where spawning and subsequent settlement could potentially be explained by environmental conditions during these time periods. For example, previous studies have indicated that settlement events commence about 3 wk after a decline in temperature triggers spawning (Saoud et al. 2000). Two weeks prior to the settlement events were termed 'spawning events' and 4 wk before were termed 'prior to spawning'.
Data were statistically modeled independently for 2014 and 2016 due to large differences in settlement between years and the addition of settlement plate sites in 2016. Due to overfitting of models, settlement sites and sampling periods with few (n ≤ 3) spat settled were removed from models, such that data from MB4 and MB6 were removed from 2014 temporal settlement models, data from MB5 were removed from 2016 spatial and temporal settlement models, and data from MB7 were removed from 2016 temporal settlement models. For comparison between years, spatial and temporal settlement models were run from mid-May through mid-August to cover periods of common spat settlement and data acquisition.
All models were checked for the assumptions of linearity, normality, and heteroscedasticity of residuals. If assumptions were not met, appropriate model relationships, i.e. polynomial relationships, were used or data were log transformed. GLMs were re-run without significant factors and/or interactions in cases where factors and/or interactions were not significant. Final models were chosen using partial likelihood ratio tests. Week as a unit of time was used as a continuous variable in all models. An α of 0.05 was used for all tests, and error is presented as ± SE. GLMs and the negative binomial GLM were run using the MASS package in R (Venables & Ripley 2002). All analyses were conducted in RStudio v.1.1.453 and R v.3.5.0 (R Core Team 2017).

Shell collection
Spat shells (1.3 ± 0.1 mm) collected from 2016 settlement plates were used for TE analysis. Spat from 2014 could not be used due to low and inconsistent spat settlement and shell fragility. Spat shells were collected at all settlement plate sites in Fig. 1 except MB7 due to inconsistent plate recovery. Adult oysters (5.0 ± 0.6 cm) were collected from The Nature Conservancy (TNC) Mobile Bay Restoration Project intertidal oyster reefs (sites EMS4, EMS5, MB2, and MB3) and opportunistically during field sampling from sub tidal settlement plate sites (shell symbols in Fig. 1). Sample size (n = 3−5) at each site was limited by the number of available adult oysters, which are largely only present in southeastern Mobile Bay (indicated by red historical area in Fig. 1). Only adult oysters collected after September 2016, which were in the water for the duration of spat settlement sampling, were used for analysis to ensure that the time period of TE ratio acquisition overlapped between adult and spat shell.

Shell preparation
Spat shells were collected off settlement plates using a tungsten probe. To remove excess organic matter (i.e. mud and bryozoans), spat shells underwent 3 cycles of rinsing with ultrapure water and drying at 60°C. Additional chemical treatments were not used to avoid altering the elemental composition of biogenic carbonate (Krause-Nehring et al. 2011) or degrading the fragile larval shells (Kroll et al. 2016). Adult oyster shells were radially sectioned (hinge to outer edge) using a Buehler Isomet ® 1000 Precision Saw using Isomet™ Diamond Wafering Blades (152.4 × 5.1 × 12.7 mm) to create a 2 mm thick section. For analysis, all shells (whole spat and adult sections) were mounted to slides using Scotch™ double-sided tape.

Laser ablation inductively coupled plasma mass spectrometry
Shells were ablated with an Nd: YAG NWR213 (Electro Scientific Industries) laser ablation system coupled with Agilent Technologies 7700× inductively coupled plasma mass spectrometry (ICP-MS). Spat were ablated across the whole shell (umbo to outer edge; Fig. 2) with a 10 μm spot size, scan speed of 10 μm s −1 , 25% intensity, and 10 Hz following pre-ablation at 20 μm spot size, 50 μm s −1 scan speed, and 20% intensity. Laser intensities and scan speeds were chosen to ensure that the laser did not burn through the entire depth of spat shell, and all spat were checked under the microscope post-ablation. Adult shells were ablated on a transect from inner to outer surface ( Fig. 2), perpendicular to internal lines of growth at a 25 μm spot size, scan speed of 15 μm s −1 , 45% intensity, and 10 Hz following pre-ablation at 40 μm spot size, 50 μm s −1 scan speed, and 40% intensity. Three transects were run per shell, with a laser warm-up and washout delay of 30 s between transects. NIST 612 glass (National Institute of Standards and Technology) and MACS-3 calcium carbonate (US Geological Survey) reference standards (2 transects with identical parameters to shells) were run every hour and at the beginning and end of sampling.
Shells were sampled for 24 Mg, 88 Sr, 137 Ba, 208 Pb, 57 Fe, 111 Cd, 63 Cu, 55 Mn, 66 Zn, 75 As, 60 Ni, 51 V, 52 Cr, and 59 Co with 43 Ca as the internal standard. Cd and As were not detected in spat shells, while Cd, As, and V were not detected in adult shells. Elements were chosen based on previous geochemical tracking studies conducted in estuarine environments (Becker et al. 2005, Kroll et al. 2018. Data reduction and limits of detection were calculated in Iolite software (version 3.63) using MACS-3 as the reference material and Trace_Elements as the data reduction scheme. Elements are reported as metal (Me):Ca in mmol mol −1 . Spat transects (~1300 μm) were averaged along the first 150 μm from the umbo and the last 150 μm near the shell margin, representing larval and settled shell, respectively. Concentrations in adult shells were averaged for the whole transect (multi-year average) and 300 μm of most recent growth (approximately a single year). Distance stratification along the shell was based on known growth relationships for adult and larval oysters in the region (Gallager et al. 1986, Kirby et al. 1998, Rikard & Walton 2010.

TE statistical analyses
To determine if spawning and settlement sites had different elemental signatures, TE ratios of the larval versus settled shell were directly compared using 2-way MANOVAs (multi-elemental comparisons; Pillai's trace test statistic) with site and shell type as factors, followed by univariate ANOVAs (individual elemental comparisons). To determine if larval (i.e. natal site signature), settled (i.e. settlement location signature), and adult shells had sitespecific TE signatures, linear discriminant ana lyses (LDAs) were run for larval, settled, and adult shells separately.
To determine if the elemental composition in larval and settled shell grown during the study was similar Adult shell transects were run perpendicular to growth lines from inner to outer surface. Three transect lines were run for each shell within the spawning season, separate LDAs were run for larval and settled shells sampled during 3 spat collection time periods (May−June, July−August, and late-August−September 2016) and compared to determine variability in site separation. To determine if the elemental composition in adult shell grown during the study (most recent 300 μm) was similar to all previous years of shell growth, separate LDAs were run for the recent shell (single year) and whole shell (multiple years) and compared to determine variability in site separation.
To determine if adult shells could be used as a proxy of natal location TE signatures and provide probable connectivity patterns in the MB−EMS system, we determined larval origins by classifying larval shells of settled spat to site-specific TE ratios in the inner margin of adult shell. To predict larval origins, larvae were classified to sites using the most recent adult shell TE ratio LDAs that were run with (1) TE ratios chosen from larval shell LDA results (except V because it was not present in adult shell) and (2) without the northernmost MB site (MB8) from where larvae were unlikely to originate and could confound results (Gomes et al. 2016). Larval predictions were also attempted with whole-shell TE ratios, but differences among sites were greater for the recent shell LDA, and thus the whole-shell LDA was not used. As a validation for connectivity results and to ensure that differences between aragonite and calcite shell did not influence results, we re-ran larval predication analyses without the major 2 + cations (Mg, Sr, and Ba) that are readily substituted for Ca in the shell, and connectivity results were largely similar between the 2 results.
For all statistical tests, TE ratios were Box-Cox transformed. Multivariate outliers were identified by plotting robust squared Mahalanobis distances of the residuals against the corresponding quantiles (Q−Q plot) of chi-squared distribution. Multivariate normality was tested using a multivariate Shapiro-Wilks test. We had an equal number of groups and small sample sizes, and thus homogeneity of variancecovariance (Box M test) was not tested and LDAs (MASS R package; Venables & Ripley 2002) were chosen over quadratic discriminant function analysis. Adult shell LDAs were run in backward stepwise fashion using the stepclass function in 'klaR' package in R (Roever et al. 2018) to determine the suite of elements to include in the final models. Larval and settled shell of spat LDAs were run in a forward stepwise fashion using the Wilks lambda statistic to determine the order of variable entry, and F-statistic probabilities were used to evaluate model improve-ment using the 'klaR' package in R (Roever et al. 2018). Backward LDAs did not perform well in spat analyses because of high collinearity among most TE ratios. Prior to running final LDAs, MANOVAs were run using Pillai's trace test statistic to ensure site separation in multi-elemental TE ratios to validate the use of an LDA. Prior probabilities were computed with equal group sample sizes for all LDAs. Jackknife reclassification success was used to determine classification accuracy (i.e. the success rate of the LDA to assign shells to a site). Standardized coefficients were used to assess the relative contribution of each TE ratio to determining site separation. Because the LDA model was forced to predict all larvae to a potential sampled spawning site while not all potential spawning sites in the system were sampled, when we assigned larvae to a natal source (i.e. adult shell proxy source), any larvae with a group probability > 0.9 had a strong likelihood of originating from the predicted site (Gomes et al. 2016). All analyses were conducted in RStudio v.1.1.453 and R v.3.5.0 (R Core Team 2017).

Spat settlement
3.1.1. Spatial and temporal variation Settlement was higher in 2016 (mean and maximum settlement: 70 ± 19 and 1137 spat plate −1 , respectively) than 2014 (18 ± 6 and 174 spat plate −1 , respectively). In both years, settlement was higher in EMS than MB. MB sites had up to 1 ± 1 spat plate −1 in 2014 and up to 25 ± 11 spat plate −1 in 2016. In 2014, spat settlement exponentially increased westward from MB6 to EMS1, with the highest settlement at the westernmost site, EMS1 (34 ± 14 spat plate −1 ; Fig. 3A). In 2016, spat settlement did not increase continuously westward because settlement at MB1 (56 ± 23 spat plate −1 ) was as high as that at EMS1 (54 ± 23 spat plate −1 ) (Fig. 3A). EMS3, sampled only in 2016, had the highest settlement (200 ± 91 spat plate −1 ). Spat settlement increased exponentially westward from MB7 to EMS3, shifting peak settlement from west to east EMS in 2016 compared to 2014. Previous studies conducted in the MB−EMS system also showed patterns of increasing settlement westward but reported higher maximum settlement than found in this study (Fig. 3B).
An exponential increase in spat settlement began in mid-July 2014 (Day 192) and 2016 (Day 234) (dashed lines in Fig. 4). In 2014, site and week were the most important drivers of the in crease in settlement (negative binomial GLM: deviance = 49%; dispersion parameter = 0.58). Similarly, site and week in 2016 were the most important drivers of the increase in settlement (zero-altered negative binomial GLM: McFadden's pseudo R 2 = 0.13; dispersion para meter = 0.62). In both years, oyster settlement increased exponentially with time at similar rates (slope = 0.38 in 2014 and 0.34 in 2016) at all sites, but intercepts differed among sites (Table S2, Fig. 4).

Environmental attributes
Freshwater discharge differed between years (F 1,303 = 93.03, effective df = 8.29, p < 0.0001, Gaussian process smoother). In both years, discharge was highest at the beginning of the sampling period, with discharges of 3000−5000 m 3 s −1 in May− June of 2014 and ~2500 m 3 s −1 in early May of 2016 and discharges for the remaining sampling periods <1500 and <1000 m 3 s −1 , respectively (Fig. 5A). In both years, winds were primarily SW and SE with speeds of 3−6 m s −1 , and leading up to the settlement events, winds were primarily SW (Fig. S1).
Salinity (2014) (Table S3). Among the sites tested when comparing 2014 and 2016 salinity data (EMS1, EMS2, MB1), there was a significant interaction between year and week (F 1, 32 = 25.17, p < 0.0001). For every 2 wk period, there was a 1.4× greater increase in salinity in 2014 compared to 2016, driven by the higher discharge and resulting lower salinity in May− June of 2014 (Figs. 5 & S2, Table S4) and the sustained low discharge in July−August that in creased salinity to comparable values in August of both years. Additionally, modeled salinities were lower in the MB−EMS system in 2014, when salinities within MB were rarely >10 in May−June (Fig. 6). Temperature differed between years (t = 2.36, p = 0.01) and was lower at the beginning of 2014 (26°C) than in 2016 (29°C) (Fig. 5B).
Salinity generally increased prior to potential spawn ing and settlement events, while temperature changed little before either event. During the 2 wk period in mid to late June (Days 164−178), prior to the estimated spawning event in 2014, salinity and  Saoud et al. 2000Saoud et al. (1999 is not shown for scaling clarity temperature increased 9.5 units and 1.8°C, respectively (Fig. 5). Similarly, during the 2 wk period in late June to mid-July (Days 178−192) prior to beginning of the subsequent 2014 exponential increase in settlement (settlement event), salinity increased 4.2 units and temperature dropped by 1°C. The modeled results also indicated salinity increased in 2014 prior to the settlement event, with salinities throughout MB and EMS <10 and <15, respectively, before the event but ~20 in EMS during the settlement period (Fig. 6A). In 2016, 2 wk prior to the estimated spawn-ing event (June to early July, Days 206−220), salinity increased by 3.7 units among all settlement sites and temperature was stable within 0.1°C. Conversely, during the 2 wk before the 2016 settlement event (mid to late July, Days 220−234), salinity decreased by 2.9 units at all sites except for one of the westernmost sites in EMS (EMS2; Fig. 5), and temperature was stable within 0.2°C. In contrast to YSI measurements, the modeled salinities in 2016 rose leading up to the settlement event on Day 234, with salinities ranging from 10−25 in lower MB and 25−30 in EMS during this time (Fig. 6B). While YSI measurements did not detect large temperature variations, finer resolution hourly buoy temperature data located at a site representative of temperature throughout the system showed a 3.5 and 1°C decrease in 2014 and 2016, respectively, during the 2 wk prior to potential spawning events, and a 3 and 4°C temperature in crease in 2014 and 2016, respectively, during the 2 wk before the settlement events (Fig. S3).

TE ratios among shell types
TE ratios differed among shell types (larval, settled, recent, whole; Fig. 7). When comparing larval versus settled shell TE ratios, differences were seen between larval and settled shells of spat and among sites for all 3 time periods tested (all MANOVAs: p < 0.001; Table S5), with Sr being the only element to differ among sites for all time periods tested (Table S6). Accordingly, TE ratios differed between larval and settled shells for all elements except Cu in May−June, Mg and Sr in July−August, and all elements except Mn, Cu, and Sr in August−September. TE ratios differ ed among sites for Mg, Mn, Co, and Sr in May−June, all elements in July−August, and Sr in August−September. Sites that drove the spatial elemental differences were MB1 (at historical reef sites) in May−June and July−August, EMS3 in July− August, and EMS2 in August−September.  (Table S2)

Site-specific TE signatures
Larval and settled shell LDAs showed site-specific TE signatures during most of the time periods tested (Fig. 8A,B: only time periods with significant MAN -OVAs are shown). TE ratios in larval shells differed among sites (i.e. distinct natal sites) in May−June (MANOVA: p < 0.001) and August−September (MAN OVA: p = 0.02; Table S7). TE ratios in settled shells differed among sites (i.e. distinct settlement locations) in July−August and August−September (MAN OVAs: p < 0.001). Larval shell site differences in May−June and August−September were driven by differences in Sr and Ni (LD1 explaining 70% of site variation) and Mn and Sr (LD1 explaining 62% of site variation), respectively, while settled shell site differences in July−August and August−September were driven by differences in Sr and Cu (LD1 explaining 60% of site variation) and Sr and V (LD1 explaining 75% of site variation), respectively (Fig. 8A,B; Table S8).
Adult shell (recent and whole) LDAs showed sitespecific TE signatures, with better site distinction for recent shells, representing incorporation of TE ratios throughout the spawning season (Fig. 8C). According ly, TE ratios in shells of adult oysters differed among sites (MANOVAs: p < 0.001; Table S9, Fig. 8C) with recent shell site differences driven by differences in Sr and Fe (LD1 explaining 69% of site variation), and whole-shell site differences driven by Sr and Mg (LD1 explaining 79% of site variation; Table S10).

Temporal variation in site-specific TE signatures
Larval and settled shells had different site separation patterns (Fig. 8A,B) and classification successes (Table S11) during each of the 3 time periods tested, indicating that TE signatures varied over multiple months. For example, larval (May−June) and settled (August−September) shell differed in site separation (Fig. 8A,B) and were temporally variable (classification [100 and 89%, respectively] and jack-knifed reclassification [73 and 67%, respectively] success; Table S11).  Table S3. Arrows indicate estimated spawning events, and dashed lines indicate the subsequent start of an exponential increase in spat settlement Recent and adult shells also had different site separation patterns (Fig. 8C) and classification successes (Table S12), indicating that TE ratios varied from yearly (recent shell) to multi-year (whole shell) scales (classification [100 and 96%, respectively] and jackknifed reclassification [74 and 63%, respectively] success; Table S12).

Larval predictions
Larval origins were only predicted for the May− June time period because TE ratios had the most differentiation among sites for this period. In recent adult shells, TE ratios used for larval predictions differed among sites (MANOVA: p < 0.001; Table S13) Table S14, Fig. S4), indicating that adult shells could potentially be used as a proxy for natalsite TE ratios. Recent shell site differences were driven by differences in Sr (LD1 explaining 71% of site variation; Table S15). All larvae that settled at site EMS1 were predicted to originate from site EMS1 (probability > 0.9; Table 1). Additionally, some larvae from sites EMS3, MB4, and MB5 were predicted to originate from site EMS1 (probability > 0.8).

DISCUSSION
Exploration of settlement and connectivity can provide insight into the persistence of larval spe-cies under differing freshwater regimes. Settlement through out the MB−EMS system differed between the 2 years (i.e. high and low discharge years). High discharge during 2014 resulted in lower overall settlement, possibly due to lower survival and flushing of larvae out of the system, while oyster settlement was higher and more widespread during the low discharge year (2016). While we would have liked to test connectivity in both years, high discharge in 2014 resulted in larval shells that were too fragile for TE analyses. Additionally, while our results for larval predictions should be taken with caution due to the use of classifying larvae to adult shell, predictions did indicate EMS as an important larval source, which is in line with settlement results (i.e. settlement was highest in EMS for both years) and previous studies, showing the promise of this method. Furthermore, if our connectivity results are to be trusted, we can postulate that connectivity was likely reduced in the high discharge year compared to low discharge year, where it is likely that only a few reefs in the system were spawning, as evidenced by lower and more spatially confined settlement in 2014. These results have implications for the long-term settlement and larval connectivity of this system, especially as climate change introduces large-scale changes in precipitation and runoff patterns to coastal systems (Vörösmarty et al. 2000).

Spat settlement
Larval transport processes are inherently stochastic (Pineda et al. 2007) and result in variable spatial and temporal larval settlement patterns (Siegel et al. 2008). Accordingly, oyster settlement in the MB− EMS system was variable among sites, through sampling seasons, and between years, showing overall patterns similar to previous observations (Hoese et al. 1972, Lee 1979, Saoud et al. 2000, Kim et al. 2010). Collectively, this study  Table 1. Predictions of larval origins (May−June 2016) from the larval origin prediction linear discriminant function analysis using trace element (TE) ratios in recent adult shell (i.e. proxy for natal site TE ratios). Site indicates spat collection site and predicted site indicates adult shell collection site where larvae were predicted to originate from. Probability of group classification is the probability of larvae correctly originating from a site. Bold indicates the highest probability of site origination and thus indicates the predicted site. Each row represents an individual larva shell and previous observations indicate a persistent gradient of in creasing spat settlement westward from MB into EMS during the past ~40 yr (Fig. 3B). In this study, however, peak settlement was lower overall than in previous studies, with the highest settlement within EMS in a known productive commercial harvesting area in Alabama (EMS3; VanderKooy 2012). Kim et al. (2010) also ob served highest settlement in this commercial harvesting area during one yearly survey, while the remaining surveys had highest settlement in western EMS. Similarly, Saoud et al. (2000) found 7× greater spat settlement in this harvesting area in 1999 compared to 1998 but did not measure settlement in EMS, limiting conclusions about EMS settlement during that time (Fig. 3B). Settlement timing (June−October) also occurred within the range rep orted in other studies in the MB−EMS system (Hoese et al. 1972, Lee 1979, Saoud et al. 2000, Kim et al. 2010. Observed peak settlement coincided with published summer peaks (late July− early August) in settlement (Hoese et al. 1972, Lee 1979, Saoud et al. 2000. Other published studies have reported summer and fall peaks in settlement. While this study did not assess fall settlement, we did find peak settlement occurring in late July and late September in 2014 and 2016, respectively. In this study, spat settlement did not taper off during the sampling periods in either year, even through the end of September during 2016, suggesting that additional settlement could have occurred later in the season. While variability in the timing and location of spawning and number of larvae necessarily influence interannual variation in settlement, this and previous studies suggest that spatial and temporal settlement patterns within the MB−EMS system have been largely consistent across decades. Yet lower spat settlement observed in recent compared to previous years may lead to a loss of oysters in the future if settlement continues to decline. Changes in salinity in response to varying freshwater discharge conditions in 2014 and 2016, and to a lesser extent in temperature, likely contributed to differences in settlement and potentially spawning patterns during this study. Changes in salinity and temperature preceded spawning and peak settlement by 2−4 wk, consistent with findings in previous studies (Hoese et al. 1972, Hayes & Menzel 1981. Spe cifically, a salinity increase of 4−10 in both years and a temperature increase of 2°C in 2014 prior to the potential spawning event likely triggered spawning in the MB−EMS system. Furthermore, settlement in the high discharge year (2014) was only seen in EMS potentially because high discharge transported lar-vae westward (Kim et al. 2013) and/or spawning and survival of larvae were limited in lower salinity areas of the MB−EMS system (Gancel et al. 2019). In contrast, during the low discharge year (2016) settlement was higher and more spatially widespread. Thus, the spatial and temporal scale of settlement changed between low and high discharge years.
It is important to note that the discharge rates, salinities, and corresponding temperatures measured in 2014 and 2016 do not represent historical extremes. The observed environmental conditions were within the ranges reported in previous settlement studies in the MB−EMS system (Hoese et al. 1972, Lee 1979, Saoud et al. 2000, Kim et al. 2010. Lower spat settlement under typical environmental conditions could indicate a loss of broodstock or a lack of suitable substrate for settlement. Future studies also could consider other environmental factors, such as pH and dissolved oxygen, which may be linked to freshwater discharge or interact with salinity (Patterson & Carmichael 2018) to affect survival and are of growing importance to larval demographics in global oceans (Peguero-Icaza et al. 2011, Gerber et al. 2014).

TEs
TE ratios recorded in oyster shells showed distinct site signatures, meaning that shells collected from certain sites could be distinguished from one another. Settled shell TE ratios could distinguish sites 11 km apart (e.g. MB1 and MB4), while adult shell TE ratios could distinguish sites 2.5 km apart (e.g. EMS1 and EMS2) (Table 2, Fig. 8). While larval shells also had distinct site TE ratios, we cannot conclude that this result is representative of larval collection sites because larvae were likely transported during their larval stage, homogenizing TE signatures from various locations. Thus, without larval in cubation experiments to ground-truth these results, caution should be given to our larval shell spatial results. Furthermore, larval and settled shell sites were harder to distinguish from one another using TE ratios, as evidenced by the similarity of TE ratios among multiple sites (Fig. 8A,B). This could be ex plained by collinearity of TE ratios in larval and settled shells, possibly due to high freshwater input into the system, which is known to increase variation in TE ratios due to mixing of environmental gradients (Miller et al. 2013a, Kroll et al. 2016. Spat shells likely were more affected by changes in freshwater input compared to adult shells due to (1) a shorter time spent (~2−3 wk larval period) in the environment resulting in less reliable incorporation of elements into the shell and (2) a larval stage that ex periences multiple water masses during transport increasing the likelihood of encountering freshwater and potentially homogenizing TE signatures (Miller et al. 2013b). Other factors that may affect TE resolution in spat shells by affecting intake, assimilation, and retention of elements among sites include (1) growth rate and physiological state, such that higher elemental incorporation occurs at faster growth rates (Carré et al. 2006) and no elemental uptake occurs during growth cessation (Schöne 2008); (2) food composition that influences the availability of elements for incorporation into shell (Thé bault et al. 2009); (3) shell matrix composition, because elements are in corporated differently between calcite and aragonite (Lorens & Bender 1980, Weiss et al. 2002 and (4) ontogenetic shifts (Strasser et al. 2008, Génio et al. 2015. Sr concentration, which has been documented to vary with salinity in taxa in many systems, consistently contributed to site separation in this study. Sr was the most important TE ratio in larval and settled shells of spat to discriminate among sites in this study, and Sr was not correlated with other TE ratios. No col linearity between Sr and other elements suggests Sr was likely incorporated differently than other TEs (Lazareth et al. 2003). Sr is a commonly used salinity indicator (Dodd & Crisp 1982) despite some variation among stu dy systems and species (e.g. Secor et al. 1995, Kroll et al. 2016. Here, we found that Sr in settled shells had a weak positive relationship to salinity, but larval shells did not (re sults not shown: R 2 = 0.10, p = 0.05), which is not surprising given the potential move-  Lazareth et al. 2003). Localized sources of TE ratios, such as industry, agriculture, groundwater, and river discharge, may further contribute to variation among sites (Charette & Sholkovitz 2002, Cor di et al. 2003, Carson et al. 2013. Overall, the potential for relatively high variation in TE ra tios in larval and settled oyster shells within sites on short timescales (weeks) suggests the need to determine site-and time-specific elemental signatures for each study. This need may be particularly great in freshwaterdominated and urbanized estuaries like the MB− EMS system, which receive elements from multiple and poorly defined upstream sources. Recent versus whole adult shell also showed different TE ratios, indicating that different TE ratios are incorporated during different time periods (Fig. 8C). Temporal variation in TE ratios makes it difficult to use adult shells as a proxy for natal site signatures because TE ratios change through time and the time period of TE ratio acquisition needs to match up with the larval shell acquisition period. In this study, recent shell, representing roughly 1 yr of growth, was able to provide better classification success and site distinctions than whole shell, representing multiple years of growth. These results are promising because the recent shell represented the current spawning season time period and thus was more representative of TE ratios present in the system during larval transport and settlement than the whole shell. While not perfect due to recent and larval shell integrating different time frames (months versus weeks, respectively), adult shells show some promise for being proxies for natal reference signatures because there was adequate site distinction using recent shells. We recommend that future studies perform larval outplant studies to better define those natal signatures and compare results with adult shells to better assess the utility of adult shells as natal proxies.

Larval connectivity
While our connectivity results should be taken with caution due to mineralogical differences and different integrated time periods between shell types, the use of adult shell as a proxy for natal site reference signatures did provide evidence of probable connectivity patterns between adults and larvae in the MB− EMS system. There was high self-recruitment and recruitment to nearby locations among larvae collected from EMS, where oysters are historically harvested (May 1971) and larval retention is known to be high (Kim et al. 2010, this study). Although self-and nearbyrecruitment was seen among EMS sites, probabilities were highest (> 0.9) for the site farthest west in EMS; other predicted sites had probabilities of > 0.5, indicating potentially high larval mixing in EMS or larvae originating from natal locations outside the study area (Table 1, Fig. 9). Interestingly, larval shells collected from lower and lower-mid MB near the MB ship channel had some larvae predicted to originate from the site farthest west in EMS, suggesting that oysters in EMS are important larval sources to this system and may be the primary source of larval oysters to some of the most productive harvest areas in the Gulf of Mexico (Zu Ermgassen et al. 2012).
Other studies have concluded that larvae are unlikely to pass from west to east MB due to flushing through the mouth. Under certain physical conditions, i.e. low river discharge (< 500 m 3 s −1 ), S to SW winds, and tropic tides, however, small numbers of larvae may be transported to the lower-mid MB re gion against the expected dominant flow pathway (Kim et al. 2013). Conditions during May−June 2016 were conducive for larval transport from EMS to lower-mid MB (discharge: < 298 m 3 s −1 ; winds: SW 4−8 m s −1 ). Additionally, when discharge conditions are <1715 m 3 s −1 , tidal transport may be important for the lower MB−EMS (Webb & Marr 2016), and larvae may use selective tidal stream (i.e. biological) transport to move from EMS to MB across multiple tidal cycles (Wood & Hargis 1971, Newell et al. 2005. Previous studies concluded that although biological transport is negligible to the overall pattern of larval transport in MB−EMS it does contribute to larval re tention in EMS (Kim et al. 2010), which is in line with the results of this study. Thus, a combination of low river discharge, winddriven circulation, and tidal transport could transport larvae from EMS to lower-mid MB according to TEbased model predictions. It is also possible that similarities in the physiochemical environment, such as salinity and temperature, among sites led to difficulty in discriminating sites using the LDA model (54% jack-knife classification success) or adult shell was not an adequate TE reference signature. Overall, preliminary larval predictions showed self-recruitment and major connectivity to the EMS region.

CONCLUSIONS
Settlement data in conjunction with geochemical tagging data can provide information on larval origins and population connectivity that are critical to define priority areas for settlement and recruitment in freshwater-dominated estuarine systems. Environ-mental attributes, particularly salinity and temperature, mediated settlement patterns by affecting the magnitude, timing, and location of potential spawning events. In the MB−EMS system, these key environmental attributes were largely driven by freshwater discharge, winds, and tides. Lower overall settlement in recent years, however, suggests a po tential for loss of oyster populations through de creased spat settlement and loss of adult oyster habitat, if this trend continues. Higher spat settlement in conjunction with larvae predicted to originate from EMS indicate that the EMS could be an important source of larvae to the region and may be of special interest for future restoration or propagation efforts. This work highlights that a combined biological and geochemical approach can help identify and predict how subsequent adult distributions may change through time in highly dynamic environments. This study provides a baseline and approach for monitoring this type of future change in the MB−EMS system that can be applied to other regions of the world where freshwater inputs are high or increasing due to habitat alteration or climate change.
Due to global changes in precipitation patterns, melting of ice due to warming temperatures, and anthropogenic water diversions, freshwater inputs to coastal systems are expected to change substantially in coming years (Milliman et al. 2008, Dai et al. 2009). Large-scale increases in freshwater discharge into already freshwater-dominated systems may push suitable salinity habitat farther down-estuary, reducing or altering recruitment and connectivity patterns of many species with salinity-dependent pelagic larval stages (Powell et al. 2003, Shoji et al. 2006. Con versely, decreases in freshwater discharge to estuarine environments can cause saltwater intrusion and also change the salinity regime and habitat for salinity-dependent pelagic larval stages (e.g. Apalachi cola Bay, Wilber 1992; San Francisco Bay, Cloern & Jassby 2012). This element ratios were distinct among sites. In (A), larger circles indicate higher average spat settlement, and thicker hatched arrows indicate higher freshwater discharge. In (B), thicker arrows indicate higher probability of correctly predicting larval origins from a natal site study demonstrates how changes in freshwater discharge between years can result in different settlement patterns. Studies looking at the persistence of commercially important species that have pelagic larval stages need to consider how future long-term changes to freshwater inputs will affect long-term settlement and connectivity.