Reproductive Phenology of the Eastern Oyster, Crassostrea virginica (Gmelin, 1791), Along a Temperate Estuarine Salinity Gradient

Low salinity can negatively affect reproduction in estuarine bivalves. The spatial and temporal extents of these effects are important to inform models of population dynamics, environmental risk assessments, restoration efforts, and predictions of climate change effects. A hypothesis of delayed gametogenesis for oysters at low salinity sites was tested relative to their higher salinity counterparts in downstream experimental cages. In 2018, the timing of gametogenesis and spawning was observed June–August for 2-year-old oysters from three distinct ancestries (native, hatchery, aquaculture), outplanted at age 1 month along the salinity gradient (3–30 psu) of a temperate estuary. A second season of data was collected in 2019 from a 3-year-old aquaculture line and mixed-age native adult oysters dredged and transplanted 1 year prior. Dermo was tested in 2019 and prevalence was 1.3% (n = 240). Gametogenesis and spawning were retarded for all ancestries at low salinity relative to higher salinity sites during July and August. The reverse pattern was found in June, with low salinity sites having more advanced gonad index than at a high salinity site. This difference in average gonad index was 2.65 vs 1.46, respectively, for the native line and 2.62 vs 2.08 for aquaculture. Low salinity seemed to not only induce earlier gametogenesis in June, but also extended the reproductive season relative to higher salinity sites. Among oyster ancestries, the aquaculture line stood out as having 30–48% lower gametogenic synchrony within sites, but only in 2018. Because the native oysters used in this study have been restricted to low salinity conditions for many generations, demonstration of their reproductive plasticity across salinities is notable and broadens the range of potential future restoration strategies.


Introduction
Estuaries are ecologically vital ecotones characterized by a salinity gradient (Telesh and Khlebovich 2010). Average salinity levels fluctuate with tide, weather, river discharge, and season, with a general rise in salinity as proximity to the ocean increases (Warner et al. 2005). Many species are adapted to the variable salinity conditions of estuaries, whereas others are interlopers that take advantage of estuarine habitats during one or more life stages (Able 2005;Lellis-Dibble et al. 2008). Among estuarine adapted species, physical niche is often described using a habitat suitability index constructed from viability or growth rate across different combinations of estuarine environments (Barnes et al. 2007;Swannack et al. 2014;Linhoss et al. 2016). These indices help guide habitat protection priorities and restoration planning by modeling how species would fare in specified areas. There are two aspects of population variation that are rarely accounted for with these indices: (1) phenotypic plasticity and (2) genetic differences among local or supplemented populations. Phenotypic plasticity can dramatically broaden the realized niche of a population, so experiments that inform performance indices should ideally include a full range of acclimation effects to accurately predict performance at the margin. Genetic differences among local populations, originating from within-generation selection rather than multigenerational adaptation (Marshall et al. Communicated by Eric N. Powell. 1 3 2010; Sanford and Kelly 2011), have potential to affect population resilience through variation in life history traits, including the degree of phenotypic plasticity (Eierman and Hare 2015).
With metapopulation structure as our within-estuary model, and given the added productivity and resilience expected from portfolio effects when component populations have somewhat independent life history or population dynamics (Lipcius et al. 2008;Schindler et al. 2010), it is important to document performance variation among the diverse habitats within an estuary. Here, we focus on a trait that is poorly known for many estuarine organisms-gametogenic and spawning phenology relative to the salinity gradient. Spawning phenology differences across an estuarine metapopulation could have a large impact on recruitment and population dynamics because the timing of larval production, relative to hydrodynamic, microalgal, and community ecology dynamics, can be expected to constrain or promote larval survivorship (Starr et al. 1990;Morgan 1995).
Our study organism is the eastern oyster (Crassostrea virginica, Gmelin, 1791) because of its keystone role, contributing valuable ecosystem services to estuarine communities such as water filtration, benthic pelagic coupling, and a structurally complex benthic reef habitat used by many commercial fish species (Coen et al. 2007;Beck et al. 2011;Bricker et al. 2020;Rose et al. 2021). The synergistic mix of coastal degradation, overharvesting, disease, eutrophication, and climate change threatens to degrade oyster populations to functional extinction in some regions, if it has not done so already (Beck et al. 2011). While salinity is considered to be a weak or negligible environmental factor for triggering spawning in eastern oysters relative to temperature and chemical cues (Thompson et al. 1996), its effect on prespawning gametogenic development and overall reproductive timing is much less studied. More fully understanding the effect of salinity variation on reproductive phenology has the potential to inform eastern oyster population dynamic models, environmental risk assessments, and restoration efforts and improve predictions of the effects of climate change (Barnes et al. 2007;Marshall et al. 2010;Levinton et al. 2011). For these reasons, this study primarily examines the effect of natural estuarine salinity variation on the reproductive phenology of the eastern oyster.
Unfortunately, many estuaries no longer have extant populations of eastern oysters along the entire salinity gradient. In fact, the desire to restore full metapopulation structure within estuaries provides an important motivation for understanding how salinity variation affects phenology. The effects of salinity variation on overall reproductive timing have received little attention outside of aquaculture settings. It is well-established that the timing of spawning for eastern oysters is closely correlated with rising temperatures, particularly in temperate estuaries where reproduction has a more pronounced seasonality (Loosanof and Davis 1952;Davis and Calabrese 1964;Bourlès et al. 2009). The interactive effects of salinity and temperature influence multiple aspects of oyster life history traits such as mortality, growth, and survival at all life stages (Davis and Calabrese 1964;Rybovich et al. 2016;McFarland et al. 2022). Multiple studies have shown that low salinity environments (< 6 psu) can inhibit or depress gametogenesis in oysters (Butler 1949;Loosanof 1953;Shafee and Daoudi 1991;Shumway 1996;Honig et al. 2014;Volety et al. 2017), but the effect of salinity variation specifically on reproductive phenology has been much less studied. Butler (1949) addressed how natural salinity variation affects gametogenic progression using histology, yet that case involved unusually protracted and extreme low salinity conditions in the upper Chesapeake Bay rather than the average salinity gradient experienced by an estuarine oyster metapopulation. The strongest data relating salinity to reproductive phenology that we are aware of come from a 15-year time series of monthly oyster collections at five sites along the Caloosahatchee River Estuary, Florida, ranging from 12 to 28 average salinity in September (the annual salinity low point). Based on gonad index averages over 15 years for wet, moderate, and dry subsets, McFarland et al. (2022) found that the upper-most site had the earliest gametogenic development and the longest overall reproductive season. Dermo disease was negatively correlated with salinity and average gonad index, making it difficult to infer the separate phenology effects of these two factors (McFarland et al. 2022).
Dermo disease is caused by the pathogen Perkinsus marinus. Moderate to heavy infection reduces an oyster's investment in reproduction and shell growth, with the potential to cause mortality and devastate oyster populations (Chu and La Peyre 1993;Powell et al. 1996). This disease poses a pervasive risk to populations across the eastern seaboard and is thus of great concern to management and restoration planning (Andrews 1988;Dittman et al. 2001;Mackenzie 2007). Non-lethal infection intensities can have significant negative effects on gonad size and gametogenic development (Dittman et al. 2001). Dermo disease risk is largely limited to high and moderate salinity environments (≥ 12 psu) because of the pathogen's environmental limitations (Powell et al. 1996;Levinton et al. 2011). Thus, while oyster gametogenesis is inhibited in extreme low salinity environments, such areas can also serve as a refuge from disease (La Peyre et al. 2003. Documenting Dermo prevalence can not only help to interpret potential mechanisms limiting gametogenesis between high and low salinity sites, but it also provides a critical baseline metric to inform restoration planning and monitoring the health of restored populations. The Hudson River Estuary (HRE) provides a particularly interesting case study. This estuary includes multiple named water bodies in the brackish zone above and below New York City (Fig. 1). The HRE was once home to a highly abundant, economically valuable population of eastern oysters (Franz 1982;Kurlansky 2006). Overharvesting, water pollution, habitat destruction, and urbanization led to a precipitous decline of the population in the early twentieth century (Kurlansky 2006;Mackenzie 2007). Today, a self-sustaining eastern oyster population is largely absent from the HRE except for a small remnant population in the Tappan Zee-Haverstraw Bay (TZHB) portion of the estuary near Irvington, NY (Fig. 1), where average salinity is low (3-12 psu) (Medley 2010;Starke et al. 2011). Recent unpublished data (Hare) suggest that the TZHB population, while reproducing annually and self-recruiting consistently McFarland and Hare 2018;AKRF Inc. et al. 2021), is barely contributing any spat recruitment downstream that could help recover a self-sustaining population in lower parts of the estuary where water quality has improved and restoration efforts abound (Stinnette et al. 2018;McCann 2019). Based on published habitat suitability indices for C. virginica, the persistence and reproductive output of the TZHB population at very low salinities are enigmatic but provide hope that when restoration broodstock is needed, local native oysters will be considered a viable option. Our goal was to test the spawning capacity of TZHB oysters when transplanted to other parts of the estuary, and measure salinity effects on gametogenic timing.
To increase the generality of this study, we compared gametogenic phenology of oysters from three population sources. Choosing an oyster source for restoration typically requires evaluation of potential trade-offs with respect to restoration efficacy, involving logistics, availability, cost, genetic diversity, and relative fitness in the specific restoration environment Plough 2019, 2022). There currently are no relative performance data available on gametogenesis and phenology for oysters with distinct ancestries, even though these phenotypes can vary due to either local adaptations (genetic factors) or acclimation to distinct environments (epigenetic factors). Local native transplants (source 1: "native") are often recommended to maximize genetic diversity while minimizing maladaptive gene flow (Camara and Vadopalas 2009) and involve transplants from naturally recruiting populations within the same estuary or region. For simplicity, we will refer to naturally recruited oysters as "native," acknowledging that they are not necessarily pristine or genetically non-admixed with aquaculture lines at present or in the past. Native transplants can potentially be done at any post-settlement life stage, and do not involve hatchery propagation. Hatchery-propagated stocks for restoration or supplementation can derive from wild broodstock spawned in the hatchery (source 2: "hatchery") or from artificially selected lines used in aquaculture (source 3: "aquaculture"). Most US eastern seaboard hatcheries producing oysters sell primarily to aquaculture growers and do not typically work with wild broodstock, whereas restoration hatcheries and some research hatcheries spawn and propagate wild oysters for restoration (Hornick and Plough 2019). Thus, there can be many practical trade-offs affecting the decision to use one oyster source or another for restoration.

3
A recent survey of New York/New Jersey projects in which oysters were seeded on restored benthic habitat revealed that 19 out of 20 had obtained seed from an out-of-state hatchery, with at least partial use of local broodstock in only 3 cases (McCann 2019).
Here, we compare the timing of gametogenic development and among-oyster synchrony for these three oyster ancestries along the salinity gradient in a temperate estuary. By quantifying gonad development at locations ranging from 3 to 30 psu, we test for the effects of salinity and ancestry on reproductive phenology. Our first hypothesis is that gametogenesis is delayed in the early-reproductive season (June) at sites with lower salinity relative to oysters at the higher salinity sites (Butler 1949;Shafee and Daoudi 1991;Shumway 1996;Honig et al. 2014;Volety et al. 2017). Comparing oyster ancestries, our hypothesis is that the native line, from the local TZHB population, will demonstrate the most distinct phenology out of the three lines. More specifically, the long history of the TZHB population at an isolated, low salinity location (McFarland and Hare 2018) leads us to hypothesize that local adaptation or acclimation has reduced gametogenic suppression at low salinity for native TZHB, compared to the aquaculture and hatchery lines, both of which involved hatchery-produced seed (juveniles) from moderate salinity broodstock. Accordingly, the gametogenic delay predicted by hypothesis 1 is only expected in these latter two lines. In a second field season, only one of the oyster lines could be tested across the full range of salinities, and in addition, a fourth treatment consisted of TZHB mixed-age adults dredged and transplanted to cage sites, with 1 year of acclimation before measuring gonad index. The fourth treatment primarily informs oyster restoration planning by testing the gametogenic and spawning capacity of transplants across very different salinities.

Study Species Reproductive Biology
Eastern oysters (hereafter, oysters) are protandrous hermaphrodites, maturing first as males and later becoming females, and spending some period of the winter, depending on temperatures, with undifferentiated sex (Kennedy and Battle 1964). After sex differentiation, oysters begin gametogenesis, building, and developing gonad tissue and gamete cells (eggs or sperm). The oyster is a synchronous spawner, triggered by temperatures greater than 20 °C to release gametes into the water column as nearby individuals do the same (Loosanof and Davis 1952;Kennedy and Battle 1964;Barber et al. 1991;Barber 1996;Mann et al. 2014). In temperate waters like the HRE, the reproductive period (gametogenesis and spawning) begins in the spring and ends in the fall. As it gets colder, oysters sometimes achieve a smaller fall spawn (Hayes and Menzel 1981) and then enter a quiescent stage in the winter during which gonads are dormant and mostly undifferentiated until the next spring (Kennedy and Battle 1964;Hayes and Menzel 1981;Mann et al. 2014).

Outplant Locations
During a previous study initiated in 2016, juvenile oysters were deployed in cages along the estuarine salinity gradient to measure and compare survivorship and growth over 3 years. For this study, we sampled lines from six of these outplant locations in both 2018 and 2019 (Fig. 1). The locations at Irvington Boat Club (IRV), Hastings on the Hudson (HH), and Yonkers Science Barge (SB) represent the low salinity region of the HRE (3-13 psu). The Redhook (RED) location represents an intermediate salinity (15-25 psu) in New York Harbor. Paerdegat Basin (PGB) and Kingsborough Community College (KCC) are both in Jamaica Bay (connected to the lower HRE) and have the highest salinities (20-30 psu) ( Fig. 1). Temperature and salinity were recorded as point data during each sampling trip at each location, along with semi-continuous measurements collected using YSI 600 OMS (YSI Inc. Yellow Springs, Ohio) sondes deployed at KCC and RED. Environmental data for HH, IRV, and SB were obtained from Riverkeeper Incorporated's Water Quality Program in the Hudson River Estuary (https:// www. river keeper. org/ water-quali ty/ hudson-river/), the Hudson River Environmental Conditions Operating System (HRECOS) location at Piermont Pier (across the river from IRV), and the Center for the Urban River at Beczak (sonde deployed at SB) (http:// hudson. dl. steve ns-tech. edu/ hrecos/ d/ index. shtml). Water quality data for PGB were obtained from a Harbor Water Quality dataset provided by the New York City Department of Environmental Protection (https:// data. cityo fnewy ork. us/ Envir onment/ Harbor-Water-Quali ty/ 5uug-f49n).

Oyster Ancestry
At each location, three distinct oyster lines were compared: native, hatchery, and aquaculture. The 2018 native line was obtained in August/September 2016 as natural spat recruits on bivalve shell that was deployed in the vicinity of the TZHB remnant population, and then transplanted as 1-month-old spat to various cage grow-out locations (Fig. 1). The hatchery line represents wild genetic diversity from a moderate salinity population with standard hatchery procedures presumably imposing some genetic bottlenecking (Hornick and Plough 2019). The hatchery line was produced by strip-spawning 12 males and 12 females collected from a native population in Edgartown Great Pond, a moderate salinity (15-25 psu) lagoon in Martha's Vineyard, Massachusetts. The larvae were then cultured by the Martha's Vineyard Shellfish Group and sent to the Cornell Cooperative Extension hatchery in Southold, New York, to set on micro-clutch and grow out in a nursery upwelling system for 3 weeks during August 2016 before deploying to sites across the HRE (Fig. 1). Larval culture was at 27-28 psu and the nursery was 28 psu. Finally, a moderate salinity aquaculture line was acquired as seed oysters (clutch-less spat) in August 2016 for outplant in the same cages. We estimated that all three lines were within approximately 2-3 months of the same age and were thus approximately 2 years old when sampled in 2018.
In 2019, only the aquaculture line (now 3 years old) could be sampled across the full range of salinities and compared to 2018. The 2019 samples of TZHB native oysters consisted of a different treatment: mixed-age adults (mean length = 74.57 ± 1.13 mm) were dredged in June 2018 and immediately transplanted to experimental cages at the study sites ( Fig. 1). Unlike the 2018 TZHB native oysters that developed from spat in the experimental cages at each study site (nearly full post-settlement opportunity for developmental plasticity), 2019 TZHB native oysters had acclimated as adults for 1 year at each study site before summer sampling in 2019.

Oyster Sampling and Histology Preparation
Oysters were collected and sampled for histology three times during the summer of 2018 (June 15-21, July 8-13, and August 8-13). Logistical and financial constraints limited collections to two dates in 2019 (July 28-August 1 and September 1-4), and aquaculture individuals were not available at every location in 2019. The number of individuals that could be histologically analyzed was limited by oyster survivorship as well as funding. At each sampling site visit, shell height was measured to the nearest millimeter using Mitutoyo Calipers (Mitutoyo America Co., Aurora, IL, USA) from the hinge to the farthest edge of the shell prior to dissections.
In 2018, oysters were transported live in coolers with ambient HRE water to the Fort Totten Urban Field Station Laboratory in Queens, New York (US Forest Service) where dissection and fixation of tissue sections were performed. In 2019, oysters were transported back to Cornell University live in coolers with ice packs, location-specific HRE water, and aeration for processing the next day. After cleaning and shucking the oysters, a standardized 4-mm cross section was taken to include reproductive tissue at the intersection of the labial palps and gills (Morales-Alamo and Mann 1989; Howard et al. 2004). The tissue cross section was placed in a labeled plastic cassette and preserved in Davidson's Fixative for 7 days and then 70% ethanol (Fisher et al. 1996). Samples were then sent to the Cornell University Animal Health Diagnostic Center for histology slide preparation with hematoxylin and eosin staining.

Categorical Gametogenesis Analysis
A common approach to measuring oyster reproductive status is with qualitative categorical scoring of gametogenesis progression from histological sections (Fisher et al. 1996). We scored each sample using the 5 stages defined by Lango-Reynoso et al. (2000), Barber (1996), and Kennedy and Battle (1964) (Table 1). Histology slides were examined with a compound microscope under 100 × total magnification to determine sex and gonad developmental stage. Gametogenesis was scored by KMG. A subset of oysters with borderline category assignments were also scored independently by MPH until consistent agreement was obtained.

Dermo Disease Testing
All individuals sampled in 2019 for histology (n = 240) were tested for levels of P. marinus infection. To test for the presence of P. marinus and measure intensity of infection, a small section of rectum tissue and approximately 0.5 cm 2 from the posterior gill tissue were removed from the same individuals on which histology was performed and placed in a culture tube containing 9.5 mL of Ray's thioglycollate medium (RFTM; Ray 1952Ray , 1954. Then, 0.5 mL antibiotic solution (equivalent to 500 units of penicillin G and 500 units of streptomycin per mL of medium) and 50 µL of nystatin solution were added to the culture tube. The cultures were incubated at 25-26 °C in the dark for 5-7 days. After the incubation period, all tissues were removed from each tube and placed on a clean glass microscope slide. The tissue was macerated using a razor blade and 3-4 drops of Lugol's iodine were added before applying a cover slip. Each slide was immediately examined at 100-400 × magnification for Dermo cells using the Mackin scale of 0 to 5, with 0 indicating no infection and 5 heavily infected (Mackin 1962;Dittman et al. 2001). Because infections were so sparse, several slides were sent for confirmatory readings by Iris Burt at Rutgers Haskin Shellfish Research Lab, New Jersey.

Statistical Analyses
The proportions of each gametogenic stage were calculated for each unique combination of location, month, and line. A set of ordinal logistic models were created to test the relationship between proportions of each gonad development stage and location, month, sex, and line for the 2018 sampling year ( Table 2). Given that complete or nearly complete gametogenic data were available for sites representing two environmental extremes in the estuary, northern low salinity, and southern high salinity, location was used as a proxy for salinity and included in all models. Month was included in all models because observations were made at this time interval during the typical main gametogenic season for oysters at this latitude (Kennedy and Battle 1964). All resulting possible interaction effects were included for 2018 models. The 2019 model did not include a line variable because of limited data across ancestries, and because the Native line was based on adult transplants, not juvenile transplants like in 2018. Instead, 2019 aquaculture data are described qualitatively, and only native data were included in a location + month + location*month model to best mirror the top ranked model in the 2018 model selection (Table 2). Models were ranked in terms of AIC values with the criterion that well-supported models have a delta AIC value < 2 relative to the top model (Burnham and Anderson 2002). Pairwise post hoc contrasts were performed using the "emmeans" package (Lenth 2022). All statistical analyses were performed using R and RStudio (R Core Team 2018).
To quantify the synchrony of oysters, we calculated the Shannon-Wiener diversity index (hereafter, diversity index) for each unique combination of month, location, and line (Shannon and Weaver 1949;Tuomisto 2012). The formula, where p i is the proportion of the total sample found to have each of the five gametogenic categories, captures evenness of gametogenic stages at a given time and place (Eq. 1). Used this way, the index has an inverse relationship with synchrony such that as the index increases, the level of synchrony decreases. The diversity index was compared across locations and lines using generalized linear models, following the formulas created for models 1-7 of the 2018 ordinal logistic modeling. Models were ranked in terms of AIC values, but because this diversity metric lacks replication within sites Table 2 Sample sizes and sex ratios of oysters sampled for histology (# females: # males). Locations are listed by latitude, from north to south. Locality abbreviations explained in Fig. 1 legend and  and lines, post hoc analyses were not possible. These analyses were performed using R and RStudio (R Core Team 2018).

Environmental Conditions
Across locations, temperature varied from 0.9 to 30.1 °C throughout the year (Fig. 2). There was a seasonal flux in temperature, with only slight variation between years and among locations. Treated as two regions, the northern sites (IRV, HH, SB) and southern sites (PGB, KCC) saw little difference in average temperature during any sampling month (Fig. 2). Salinity variation (Fig. 2) was non-overlapping for northern (0.9-15.7 psu) and southern sites (18.4-30.1 psu), supporting our use of these locations as proxy variables for low and high salinity, respectively. For each site, there was little variation in salinity between years, except at KCC and SB. In 2018, there was a gradual decrease in salinity at KCC (29.9 to 19.0 psu) between April and October, but in 2019, the trend was reversed, with salinity increasing from 23.8 to 30.1 psu in those months. The SB site had much lower spring

Gametogenic Development Stage
The moderate salinity site (RED) had extremely low cumulative survivorship and was not included in model testing.
For other sites, mortality prevented fully balanced sample sizes (June 2018 average n = 13, all other sampling periods average n = 20; Table 2). In both years, sex ratio was either roughly even or male-biased, with larger individuals being females as expected for this protandrous hermaphrodite ( Table 2, Table S1). Based on model selection criteria using AIC values, the top ranked model for 2018 was model 7: Location + month + line + l oca tio n*month + location*line + month*line (Table 3). Model 11 was also supported, differing from model 7 only in the inclusion of sex as a parameter. Sex was determined to be an uninformative parameter based on the criteria in Arnold (2010), so the post hoc analysis was thus performed only for model 7. In August of 2018 and in all 2019 samples, the aquaculture line at PGB became 100% male (Table 2). For the linear models of the diversity index, the top model was model 4, location + month + line + month*line, indicating that location effects on synchrony were not interacting with line or month (Table 4).
In June 2018, only native and aquaculture could be compared between low and high salinity sites and in both cases gametogenic development had advanced more at low salinity. A greater proportion of individuals were mature or spawning at the low salinity site SB (75-80%) compared to the high salinity site PGB (0-20%; native contrast, combined n = 26, p < 0.0001; aquaculture, n = 26, p = 0.0001; hatchery could not be tested; Table 5, Fig. 3a, g). The low salinity site HH showed the same trend in pairwise contrasts with the high salinity site PGB (Table 5). For the native line at high salinity, the complete lack of mature or later stage individuals contrasted with the diverse gametogenic stages observed in aquaculture at the same site, while patterns were similar among all three lines at low salinity (Fig. 3). The diversity index was highest, and thus, synchrony was lowest, for the aquaculture line at the PGB high salinity site (Fig. 3, Table 6).
In July, gametogenesis had not advanced much in any line at low salinity sites, whereas at the moderate and higher salinity sites, the proportion of spawning and reabsorbing stages had increased, where comparisons could be made to June. Thus, the June pattern of relatively more advanced phenology at low salinity had reversed by July. Phenology again differed between low salinity SB and HH sites and the high salinity PGB site (averages across lines for spawning or reabsorbing individuals 0.9% vs. 21%, respectively), with 4 out of 6 post hoc site contrasts p < 0.05 among the three strains (Table 5). For aquaculture, an even greater proportion of spawning and/or reabsorbing stages was observed at the moderate salinity RED site than at high salinity PGB. In contrast, the two low salinity sites, SB and HH, both were dominated by late active and mature stages with almost no spawners (Table 5; Fig. 3b, e, h). Comparing the ancestry lines, the most striking difference was a slightly slower phenology in aquaculture as indicated by 18.6% early active stage oysters averaged across all sites compared with 4.5% in hatchery and native combined (Fig. 3b, e, h). In fact, aquaculture phenology went slightly backward from June to July (more early active stage at all sites except PGB), even while spawning increased (Fig. 3h). All ancestry lines contained all 5 stages when tallied over all locations in July, but the among-location average diversity index was 1.2 for the aquaculture line, indicating lower synchrony than native and hatchery with 0.84 and 0.62 diversity indices, respectively ( Table 6). The July hatchery line within-locale synchrony was higher (diversity index lower) than in any other month and line (Table 6).
In August, gametogenic progression to a later mix of stages was observed in all lines, at all sites. All three lines continued to show a higher proportion of later stages at higher salinity sites, dominated by spawning and reabsorption, relative to low salinity sites where late active and mature stages were abundant (all six contrasts p < 0.0026; Table 5). Relative to July, low salinity sites in August progressed less than higher salinity sites except for the aquaculture line. The only other line distinction of note in August was the native oysters having 44% pre-spawning stages at high salinity PGB, whereas hatchery and aquaculture lines were < 15% pre-spawning stages (Fig. 3c, f, i). The diversity index was again overall highest for the aquaculture line.
The 2019 data only included 2 months, and data collection was shifted 2 weeks later relative to 2018. Two sites were included for both the high salinity and low salinity regions, but complete data were only obtained for the native line (adult transplants that had acclimated for a year in this case, as opposed to 2 years growth from spat on site for 2018 data). No hatchery line oysters were available in 2019. At all sites, oysters demonstrated gametogenic progression between the two sampling dates (Fig. 4). The post hoc contrasts between locations found little difference in stage proportions for either month (Table 7; Fig. 4a, b). This spatial gametogenic synchrony was coupled with temporal synchrony within sites based on low average diversity index values of 0.38 and 0.44 for July/August and early September, respectively (Table 8). At the high salinity sites where the aquaculture and native lines could be compared, the former had substantially more advanced gametogenic stages at the end of July (+90 % spawning and reabsorption stage for aquaculture, but only 15-45% for native), but only modest differences with the same trend at the beginning of September.   Table 2). Locations are listed in order by latitude, from north to south. Locality abbreviations explained in Fig. 1

Dermo Disease
Only 3 of 240 individuals had P. marinus spores detected (1.3%). All detected infections were in the native line at PGB during the September sampling period. Two samples had Mackin intensity levels of 0.5, and another had a level of 1.

Discussion
In one of the first studies to measure salinity effects on gametogenic phenology in a bivalve, and compare oyster lines with different ancestries, we found that in July and August of 2018, all three oyster lines had delayed gametogenic development at low salinity sites relative to high salinity sites, confirming our first hypothesis that low salinity would cause a delay. Surprisingly, the reverse was true in June 2018 for the two lines that could be compared in different salinities up and down river. June spawners only occurred at low salinity for native oysters and contributed to a weighted average gonad index of 2.38-2.92 at low salinity sites, compared to only 1.46 at high salinity. The aquaculture line showed this trend less dramatically, with weighted averages of 2.62 vs. 2.08 at low and high salinity, respectively, partly because the few June spawners were observed only at high and moderate salinity sites. The June pattern of more advanced oyster gonad index proportions at low salinity relative to high salinity was reversed in July due to gonad index stasis at low salinity versus progression toward spawning and reabsorption at high salinity. Our second hypothesis, that the native line would have the most distinct phenology characterized by earlier (less delayed) gametogenesis at low salinity, was rejected. All Table 6 Shannon-Wiener diversity index values for each unique combination of month, location, and line for 2018. Locality abbreviations explained in Fig. 1 legend and    three lines had similar gonad index distributions at low salinity in June. Instead, the most striking distinction observed in the native line was the extent of gonad index delay at high salinity in June, as indicated by zero oysters with a "maturing" or later stage (weighted average gonad index 1.46), whereas at the same site, the aquaculture line had 25% maturing or spawning oysters (2.08). By July, the native oysters had compensated and had gonad index proportions similar to the other lines. The other line distinction of note in 2018 was the consistently higher diversity of gonad index stages in the aquaculture line throughout the summer, indicating 30-48% less synchrony at any particular time or place. Low variation in temperature among sites and minimal Dermo infections detected overall highlight regional salinity differences as likely key to phenology variation among sites in 2018, but causal inferences are tentative because oyster cage outplant sites differed in other unmeasured ways. The 2019 sampling provided fewer spatial and line comparisons. Despite similar environmental conditions in the 2 years, 2019 patterns indicated greater synchrony and only slight salinity effects relative to 2018. We discuss these findings with reference to environmental variation and future restoration goals.

Month
The model that best explains 2018 gametogenic stage proportions includes month, location, line, and their interactions. Given the summer temperature variation (Fig. 2), month is the most obvious factor of expected importance since the species is considered to have locally synchronous spawning triggered by rising temperatures and conspecific spawning (Thompson et al. 1996;Bernard et al. 2011;Aranda et al. 2014). While some spawning occurred in each month (average proportion of spawners and reabsorbing individuals equal to 3.5%, 7.5%, and 46.8% in June, July, and August 2018, and 58.63% and 97.5% in July/August and September 2019), the spawning period in the HRE peaked in August and continued into September for these outplanted, caged oysters. An August peak for native TZHB oysters is supported by spat monitoring in the HRE, where the most abundant recruitment of spat (after a 2 to 3-week larval stage) has been observed in September (McFarland and Hare 2018). Our results suggest an extension of the spawning season for the HRE relative to nearby Long Island Sound where spawning peaks in July with spent oysters by August (Ford et al. 1990;Barber et al. 1991). Cold Adirondack mountain spring melt waters in the Hudson River, or differences in phytoplankton bloom timing and concentration, could be responsible for this geographic difference (Hofmann et al. 1992;Bernard et al. 2011).

Salinity
Our hypothesis of delayed gametogenic progression at lower salinities was only observed after the month of June, suggesting that the initiation of gametogenesis is less constrained by salinity variation, but that low summer salinity at northern sites relative to southern sites could be causing a delay at later gonad maturation stages. By the August 8-13 sampling period in 2018, northern low salinity sites (HH, SB) had fallen to near 5 psu and a majority of oysters were still in pre-spawning stages, a trend we expected based on literature suggesting reproductive depression caused by low salinity (Butler 1949;Shafee and Daoudi 1991;Shumway 1996;Honig et al. 2014;Volety et al. 2017). In contrast, the majority of oysters at the high salinity PGB site (approximately 25 psu) were spawned or reabsorbing for all three lines in August. These June and later summer findings in Hudson River oysters are similar to those reported by McFarland et al. (2022) for eastern oysters in the Caloosahatchee River Estuary in Florida. They found over 15 years of monthly sampling, across 5 sites ranging from 12 to 28 average salinity, that average gonad index at sites had a negative relationship with salinity. The lowest salinity site had the earliest gametogenic development; it had a 15-year average March gonad index of 2 while all other sites had averages below 1.5 (Fig. 7 in McFarland et al. 2022). We did not measure spring gametogenesis in Hudson River oysters, but initiation much earlier than June in the northern TZHB population seems unlikely given that it often experiences extended periods below salinity 5 during spring . Instead, we speculate that the advanced June gametogenesis observed here results from a rapid compensation for stressful spring conditions. Whatever the mechanism, the early gametogenic progression that produced some spawning in northern low salinity sites, coupled with a small proportion of spawners in August relative to high salinity sites, indicates that oysters Contrary to salinity effects observed in 2018, native adult transplants after 1 year of acclimation showed no significant differences in gonad index among sites in 2019. In addition to different acclimation histories, summer variation in salinity at northern sites was less in 2018, lacking a 5 psu drop in August or September of 2019 (Fig. 2). Furthermore, sampling dates were offset in the 2 years, and factors such as food availability (Luna et al. 2000;Bernard et al. 2011) or local freshwater inputs (i.e. river discharge) (Baba et al. 1999;Wilson et al. 2005) could drive interannual variation. We recommend further study using temporal replication to determine the interannual factors that can affect reproductive phenology.
In terms of outplant method, it is notable that the 2019 native dredged oysters were able to develop gametes at a variety of salinities only 1 year after transplant. Furthermore, the fact that both outplant strategies (outplanted as spat, outplanted dredged adults) resulted in oysters that were able to undergo a full gametogenic cycle is encouraging for future restoration efforts considering the potential advantages associated with restoration using locally adapted species (Hofmann et al. 1992;Camara and Vadopalas 2009;Flanagan et al. 2018;Plough 2019, 2022). This phenotypic plasticity is characteristic of oysters, but presumably has limits that need to be accounted for in restoration planning (Li et al. 2017(Li et al. , 2018(Li et al. , 2021.

Ancestry
The prediction that the native line is adapted to low salinity and therefore would uniquely demonstrate a lack of reproductive depression at low salinity was only consistent with patterns observed in June 2018. Unfortunately, only comparisons between native and aquaculture strains were possible in June. The spatial contrast in gonad index between low and high salinity sites was much stronger in the native line. The average proportion of early active and mature individuals for natives in the north was 0.08 and 0.58 respectively, compared to 0.54 and 0.0 at high salinity PGB. The overall gonad index average for natives was higher at low salinity vs high salinity sites, 2.65 vs 1.46. This trend also was detectable in aquaculture individuals but was less extreme (gonad index 2.62 in the north vs 2.08 at high salinity PGB). This suggests that TZHB native oysters may have local adaptations that allow them to rapidly compensate for gametogenic delays experienced during very low spring salinity levels, perhaps in response to overall metabolic depression (Gurr et al. 2020).
In August, there was no single oyster line that stood out in terms of phenology or salinity effects, but there were trends in the gonad index progression with native slowest (3.37 among-site average gonad index), hatchery intermediate (3.54), and aquaculture most advanced (3.82). The same among-line trend was seen both at the low salinity sites (average gonad indices 3.03, 3.17, 3.50 across lines, respectively) where gametogenic progression lagged relative to the southern high salinity site PGB (4.06,4.30,4.45). This trend for 2018 indicates that under some conditions, the native line has a more extended reproductive season than the other oyster lines. A lengthening of the gametogenic period can be interpreted as a form of depressed gametogenesis (Butler 1949;Shafee and Daoudi 1991;Shumway 1996;Honig et al. 2014;Volety et al. 2017).
The 2019 samples allowed comparison of lines only at high salinity sites, albeit confounded by different acclimation treatments. At high salinity sites, the gonad index differences between native and aquaculture lines were not large, but the 2018 trend of a more extended reproductive season for the native line also was evident in 2019. At the end of July, gonad index averages for PGB and KCC high salinity sites were 3.15 and 3.50 for native, versus 4.50 and 4.82 for aquaculture, respectively. In early September, the native line was only slightly delayed, 4.71 and 4.68 mean gonad index versus 4.85 and 5.0 for aquaculture, respectively. In fact, all sites for both lines were almost 100% spawning or post-spawning in September, indicating a spawning cycle was successfully completed.
One of the biggest differences among lines was in synchrony as measured by the diversity index. When averaged across all sites and over the entire 2018 season, the stage diversity in native, hatchery, and aquaculture was 0.78, 0.81, and 1.05, respectively. Thus, aquaculture was the outlier with less synchrony. We are not aware of other comparable data and can only speculate that long-term culture and domestication led to relaxed selection (i.e., more random change under genetic drift) on sensory systems that contribute to synchrony under natural conditions. However, the low aquaculture synchrony was not observed in 2019, so this difference among strains probably relates to distinct expressions of genotype by environment plasticity (Eierman and Hare 2015) that require further study.

Conclusions and Context for Future Oyster Restoration
This study demonstrated differences in reproductive phenology along the HRE salinity gradient that were associated with locality (low vs. high salinity regions) and oyster ancestry. These novel results in a temperate estuary suggest that reproductive phenology within an eastern oyster metapopulation is not simply a function of temperature, but also varies along the salinity gradient. The 2018 patterns were consistent with salinity effects on phenology documented elsewhere, but our interannual results (one ancestry only) caution against generalizing from the 2018 patterns because genotype by environment interactions likely have strong effects. Additional replication is necessary to further understand genotype by environment effects. Future studies would be valuable that can experimentally link aspects of phenology with specific environmental drivers. In particular, future studies should sample in Spring as well as Summer because salinity effects seem to vary at different stages, and experimental designs to test for reproductive compensation would help evaluate the consequences of stressors.
Similar to previous reports in the HRE McFarland and Hare 2018), the near absence of Dermo infections in this study makes it reasonable to infer salinity effects separate from disease interactions. Furthermore, the lack of Dermo presence across all sites supports the conclusion that disease is not the reason for the native TZHB population's current isolation at low salinities, but further study would be necessary to confirm this.
The moderate salinity aquaculture line studied here succeeded in spawning across diverse environments, attesting to its plasticity with respect to salinity. However, it also showed lower gametogenic within-site synchrony compared to the other lines, perhaps indicating domestication effects lowering a fitness-related trait. Selective breeding for commercially valuable traits can cause reduced genetic diversity and potentially reduced lifetime fitness in a natural environment, making aquaculture lines often the least desirable option for restoring oyster populations (Baggett et al. 2014). Enhancing or restoring a population with stocks that do not possess adequate genetic variation, or that have maladaptive variation, could be detrimental to overall population fitness and diversity (Camara and Vadopalas 2009;Morvezen et al. 2016;Plough 2019, 2022;Hughes et al. 2019).
With the growing threat of climate change, species inhabiting estuaries will likely be exposed to more frequent and extreme storms that can rapidly shift salinity to the edge of their tolerances. Investigating the effect of salinity on reproductive phenology is only one aspect of research that is essential to understand and possibly mitigate the potential effects of climate change on estuarine oyster populations. The native TZHB population is at the edge of habitable salinity variation with little quality habitat to expand to upstream , and it is not known to occur in any habitat further south except as sparse spat (new recruits; Medley 2010) or adults on pilings in Hudson River Park, Manhattan (Fitzgerald et al. 2020). The ability of TZHB oysters to survive and reproduce at moderate salinities, as evidenced in this study, indicates that its continued isolation north of New York City is likely a function of hydrodynamics and/or stressful downstream water quality for the mobile larval stage. If poor water quality in the lower HRE amplifies larval mortality above typical levels, then further water quality improvements will be needed before restoration with any oyster line has a chance of establishing a self-sustaining metapopulation (i.e., promoting the full life cycle) in the lower estuary. Assuming adequate water quality, transplantation of adults or spat to seed hard substrates in moderate salinity habitats could help protect the population from climate change effects by reestablishing metapopulation dynamics (Lipcius et al. 2008). Here, results demonstrate that the TZHB population provides a valuable resource for transplant restoration strategies that could leverage and expand local genetic variation.

Acknowledgements
We would like to thank Harmony Borchardt-Wier, Iris Burt, Sarah Gisler, Hannah Hartung, and Tiffany Medley for their support in this research effort through field work assistance and/or methods and data analysis consults. We appreciate the careful oyster husbandry contributed by Amandine Hall at Martha's Vineyard Shellfish Group and by Gregg Rivara at the Suffolk County Marine Environmental Learning Center hatchery. Pete Malinowski generously provided some of the oysters analyzed in this study. We also thank Captain Jim Nickels, Monmouth University, for a successful oyster dredging trip in 2018, ultimately allowing for this experiment to continue into 2019. Advice provided by Cornell Statistical Consulting Unit on data analysis and display was greatly appreciated. This study would not have been possible without housing and lab space provided by the U.S. Forest Service Fort Totten Urban Field Station. This work was only possible through the generosity of location hosts that kept experimental oyster cages secure: Irvington Boat Club, Hastings on Hudson Neighborhood Association, Groundwork Hudson Valley Science Barge, Inwood Canoe Club, Red Hook Barge Museum, Sebago Canoe Club, and Kingsborough Community College. We also thank two anonymous reviewers who provided valuable feedback to improve the quality of this manuscript.
Funding Funding was provided by a Tibor T. Polgar Fellowship from the Hudson River Foundation, a grant from the NYS Water Resources Institute, and a grant from the Cornell University College of Agriculture and Life Sciences Charitable Trust.
Data Availability All data and relevant code can be found at https:// github. com/ kaili grego ry/ oyster. git.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.