Environmental associations of cownose ray (Rhinoptera bonasus) seasonal presence along the U.S. Atlantic Coast

Identifying the mechanistic drivers of migration can be crucial in shaping conservation and management policies. The cownose ray (Rhinoptera bonasus) is a relatively poorly understood elasmobranch species that occurs along the U.S. Atlantic coast and undergoes large-scale seasonal migrations. To better understand the drivers and timing of cownose ray seasonal migration in order to inform potential management measures, we analyzed telemetry detections of 51 mature cownose rays (38 female, 13 male) tagged with acoustic transmitters in the Maryland and Virginia portions of Chesapeake Bay. Detections within their summer habitat in Chesapeake Bay and winter habitat in the vicinity of Cape Canaveral, Florida, were matched with publicly available sea surface temperature (SST) data recorded by data buoys near the areas of tag detections and with local photoperiod and day of year. These variables were used in boosted regression tree models of ray presence (all rays combined, females only, and males only) in each seasonal habitat. Models were developed for presence during the entire summer and winter season, and for the time periods of arrival and departure from both summer and winter habitats. Seasonal presence in both summer and winter habitats was associated with distinct temperature, photoperiod, and date ranges, with temperature as the most influential variable in seasonal models. In models of arrival and departure periods, southward migration (departure from Chesapeake Bay and arrival off Cape Canaveral) was strongly associated with SST for all rays and arrival in the Chesapeake Bay region after northward migration was most strongly associated with day of year. The most influential variable during the period of northward departure from Cape Canaveral differed between males (day of year) and females (SST). This suggests that mature female northward migration may be driven by temperature while male northward migration may be driven by endogenous cues. These findings provide detailed information on the timing of cownose ray arrival at, presence in, and departure from seasonal habitats and provide potential justification for including the species in cross-taxa comparative studies on migratory behavior.


INTRODUCTION
Mechanisms that signal the seasonal phenology of animal migrations have interested researchers for centuries. Debate on proximate influences was historically as diverse as animal migration itself, but was typically based on discriminating between exogenous and endogenous factors. An "internal clock" mechanism was proposed to influence bird movement as early as 1702 (reviewed in Gwinner 1996), while early research on fishes primarily focused on temperature (Goode 1879, Cunningham 1895, Gurley 1902, Chidester 1924. However, it is now widely accepted that round-trip, seasonal migration is adaptive in nature, involving both internal mechanisms and environmental factors (Secor 2015, Shaw 2016) and often generalizes across taxa (Dingle and Drake 2007).
Migration patterns in many groups have been shown to vary with sex, including in fishes (Dahl et al. 2004, Barnett et al. 2011, Lea et al. 2015, Sinnatamby et al. 2018) and birds (Maggini andBairlein 2012, Pedersen et al. 2019). For example, sex-specific differences in departure and arrival at breeding and feeding sites have been seen across taxa (reviewed in Morbey and Ydenberg 2008), with the potential for variation between fall and spring migrations (Sharma et al. 2018).
Understanding migration as an ecological phenomenon requires a broad focus that compares migration between taxonomic groups (Dingle andDrake 2007, Shaw 2016). Unfortunately, this is uncommon in traditional migration research. Shaw (2016) proposed a framework for classifying round-trip migration into three categories: "breeding," "refuge," and "tracking." In breeding migration, movement in one direction is related to breeding, while movement in the other direction is considered a return to preferential feeding habitat. In refuge migration, movement in one direction is related to escape from seasonally unfavorable conditions, and return is related to favorable feeding or breeding habitat. Lastly, in tracking migration, nomadic species follow food resources continuously (Shaw 2016).
Elasmobranch species (sharks and rays) may serve as useful models for the study of migration, ultimately providing new opportunities to compare migration across taxa. First, elasmobranchs share anatomical or life-history characteristics with a variety of taxonomic groups. Like other fishes, elasmobranchs have gills and exist in water; they have similar life-history characteristics to mammals in that they are long-lived and slow growing; and like birds, they participate in long-distance, round-trip, seasonal migrations (>1000s km, in many cases). Second, elasmobranchs collectively fit into all three round-trip migration types proposed by Shaw (2016). Finally, many elasmobranch species are distributed globally or have wide geographical distributions, allowing for the potential comparison of migration between conspecifics in different habitats. For example, tiger sharks (Galeocerdo cuvier) in the Northwest Atlantic Ocean have been shown to spend summers in highly productive high-latitude oceanic regions and overwinter on Caribbean coral reefs where they mate (Lea et al. 2015), a breeding migration in the context of the Shaw et al. framework. In contrast, tiger sharks tracked in Shark Bay, Australia (Heithaus et al. 2007), stayed in large home ranges of 1000s km 2 and movements were hypothesized to be related to prey density, a tracking migration. Even conspecifics in habitats that are in relatively close proximity sometimes display variation in migration patterns. For example, whitespotted eagle rays Aetobatus narinari in the Northwest Atlantic Ocean remained resident in the Indian River Lagoon, Florida, while conspecifics in the Gulf of Mexico displayed seasonal southern migrations thought to be associated with water temperature; a refuge migration (Degroot et al. 2021). Assessing differences in migratory patterns in conspecifics between habitats could be useful to understanding ultimate drivers of migration.
Understanding the movement patterns of elasmobranch species is also important for effective management across their migratory routes. Characterized by a slow life history including late maturity and low fecundity, elasmobranchs are vulnerable to exploitation and a thorough understanding of their migratory movements is crucial for the implementation of effective management measures (Hueter et al. 2005). For example, the cownose ray Rhinoptera bonasus is a benthicpelagic species that occurs in the Western Atlantic Ocean, from New England, USA to Argentina (Schwartz 1990, Last et al. 2016. Individuals migrate along the U.S. Atlantic Coast and enter estuaries to feed, pup, and mate. Cownose rays occupy estuaries in the northern part of their range during spring and summer and move south in the fall and winter (Last et al. 2016). Ogburn et al. (2018) found that cownose rays show consistent site fidelity to their summer habitats, which may indicate philopatry as these estuaries also function as primary nurseries. Consistent philopatric migrations may drive population structure in rays (Flowers et al. 2016). Ogburn et al. (2018) showed overwintering off Cape Canaveral and Smith and Merriner (1987) documented arrival in Chesapeake Bay by May. Pupping occurs in June or early July, and mating follows within a few weeks (Fisher 2010, Fisher et al. 2013, after which mature males may roam to other coastal areas while mature females remain within the estuary until the beginning of fall migration (Fisher 2010, Omori andFisher 2017). Although seasonal movements of Chesapeake Bay cownose rays are known, proximate cues have not yet been explored, so it is unclear which migration strategy from Shaw (2016) this species uses. Understanding why cownose rays migrate will be important for their future management, in addition to establishing the species as a model for the comparative study of migration patterns.
In the present study, we explore the relationship between sex-specific movement differences and migratory cues, with environmental variables among acoustically tagged cownose rays from the Chesapeake Bay population. Our objectives were to (1) identify differences in migratory timing over a three-year period (2014-2017) for male and female cownose rays tagged in the Chesapeake Bay and (2) determine the relative influence of environmental and behavioral variables on these movement differences. We focus on associations with the presence of male and female rays at their summer and winter seasonal migratory extents, where they are more likely to remain resident for the duration of the season (Ogburn et al. 2018).

Capture and tagging
For this study, we used acoustic telemetry data from 51 mature cownose rays (38 female, 13 male) captured and tagged within Chesapeake Bay between August 2014 and June 2016 as part of a larger study on migration and habitat selection detailed in Ogburn et al. (2018). Rays were tagged during two separate tagging efforts, with 15 (seven female, eight male) in the portion of Chesapeake Bay falling within the state of Maryland and 36 (31 female, five male) in the portion of the estuary falling within the state of Virginia. Details on tagging procedures can be found in Ogburn et al. (2018) but in general rays were captured in commercial pound net gear and internally tagged with Innovasea V13 or V16 69-kHz acoustic transmitters. Rays held for a 24-72h observation post-tagging retained 100% of their transmitters, and post-release survival was high, though approximately one-third of the tagged rays were not detected after the first year of tracking presumably due to natural mortality (Ogburn et al. 2018).

Acoustic telemetry and environmental data
Detection data for tagged rays were obtained from Innovasea VR2W passive acoustic receiver arrays within the Chesapeake Bay owned by the Smithsonian Environmental Research Center (SERC), Virginia Institute of Marine Science (VIMS), the US Navy, and NOAA's Chesapeake Bay Office Chesapeake Bay Interpretive Buoy System (CBIBS). Data from arrays elsewhere within the Chesapeake Bay region and in the vicinity of Cape Canaveral were also contributed by Atlantic Cooperative Telemetry (ACT) and FACT network researchers (Young et al. 2020). Tag detections recorded between 31 May 2014 and 18 December 2017 were used for our analysis. In areas where multiple tags are detected, occasionally signal collisions can result in a mixed transmission that matches an existing tag number, which can be recorded as a "false detection" of that transmitter. To account for this, prior to analysis tag detections were mapped and plotted against date to identify and discard tag detections that suggested unrealistic movement behaviors and locations (e.g., far up a freshwater river or extremely long distances traveled over less than a day). As in Ogburn et al. (2018), mean daily locations were calculated for each ray using the mean latitude and longitude of receiver locations that detected each transmitter that day. As with the tag detections, we mapped mean daily positions to ensure that none fell within locations suggesting unrealistic movement behavior.
v www.esajournals.org Environmental data were obtained from NOAA CBIBS and other data buoy stations located in or near the Chesapeake Bay and near Cape Canaveral (Fig. 1). Several buoys had oceanographic data, including salinity, temperature, and chlorophyll A levels, while other buoys provided only sea surface temperature (SST) (NOAA 2018). Because SST was the only environmental variable recorded by all available buoys and is frequently found to be significantly associated with movement behavior in marine species, we used it as the sole oceanographic variable in our analysis. Buoy station TPLM2 provided data for the upper Chesapeake Bay (UCB), YKTV2 provided data for the middle Chesapeake Bay (MCB), 44042 provided data for the lower Chesapeake Bay (LCB), and 44099 provided data for waters East of the Chesapeake Bay mouth (ECB). Buoy station 41009 provided data for offshore Cape Canaveral (OCC) and 41113 provided data for nearshore Cape Canaveral (NCC). These buoys recorded SST at time intervals ranging between every 6 min to every hour. Polygons covering areas expected to have similar oceanographic conditions to those recorded by the associated data buoys in the Chesapeake Bay and Cape Canaveral regions were manually delineated in ArcGIS (ESRI, Redlands, CA, USA) ( Fig. 1), and mean daily SST data were assigned to each of these polygons for each day within the study period (31 May 2014-18 December 2017). Individual rays were marked as either present or absent within each buoy-associated polygon based on whether mean daily detection locations fell within that polygon each day. The number of individual detection days (number of days each individual was present) was calculated, and the total number was used to quantify differences in ray presence and relative abundance between months and data buoy regions.

Statistical analysis
All explanatory variables were tested for multicollinearity using pairwise Pearson's correlation and variance inflation factor (VIF) analysis. Pearson's correlation analysis was conducted using functions from the R package psych (Revelle 2020), and VIF values were calculated using the "imcdiag" function in the R package mctest (Imdadulla and Aslam 2016). We used thresholds from Dormann et al. (2013) to assess whether collinearity was occurring. Multicollinearity tests were performed independently for variables in the Chesapeake Bay and Cape Canaveral regions.
We used binary boosted regression tree (BRT) modeling to evaluate the presence probability of rays in the Chesapeake Bay and near Cape Canaveral using the R package gbm.auto (Dedman et al. 2017). For this analysis, ray presence was defined as at least one individual being detected within a given data buoy region on a given day. BRT analysis is robust to most error distributions found in ecological data, insensitive to multicollinearity and outliers, and can account for pairwise interactions between variables (Elith et al. 2008, Dormann et al. 2013, Dedman et al. 2017. The procedures are briefly summarized here, but for a more detailed description of each step and the code functions, see Dedman et al. (2017) and see Elith et al. (2008) for a description of the underlying statistical theory. Regression tree models of presence probability from ray presence/absence and habitat variable data were replicated until variance ceased to change significantly between individual model runs. During this process, a predetermined proportion of the data, referred to as the bag fraction (bf), were used to train the model and the training data were cross-validated against the remaining data. Our models used SST, day of year, and photoperiod (obtained from NOAA 2018) as explanatory variables, providing a continuous oceanographic variable known to vary between years (SST), a temporal variable that differed between regions but could be assumed to be temporally consistent (photoperiod), and a temporal variable fixed across both regions that could be used as a proxy for endogenous timing (day of year). We calculated photoperiod as hours of daylight converted to integers with a factor of 100 rather than 60 (i.e., 9 h and 30 min was converted into 9.5 h). Models were run independently for presence/absence in the Chesapeake Bay and Cape Canaveral regions, and were run using detections for all rays and for males or females only. For each model, combinations of bag fraction and learning rate (lr) were tested until the combination yielding the greatest area under curve (AUC) and cross-validation (CV) score (Elith et al. 2008) was found. Tree complexity was set at 2 for all models to account for pairwise interactions between variables. Model results were used to plot the v www.esajournals.org marginal effects of each habitat variable on cownose ray presence likelihood in each region, and the percentage of tree splits accounted for by each variable, which we interpreted as a measure of that variable's relative influence on ray presence. Pairwise interactions between each variable were calculated by performing linear modeling between each pair of explanatory variables, and the residual variance of each model was used to measure the interaction strength.
Because conditions associated with migratory arrival and departure from seasonal habitats can differ from those associated with overall presence or absence, we also conducted BRT analyses using subsets of the data covering the periods of southward migration (departure from Chesapeake Bay and arrival at Cape Canaveral) and northward migration (departure from Cape Canaveral and arrival at Chesapeake Bay). These time frames began a week before the earliest arrival or departure of the first tagged ray and ended a week after the latest. The time frames selected for each data subset were, in order of occurrence post-tagging during mid-late summer, day of year 250-365 (Chesapeake Bay departure, "CB Leave"), 250-365 (Cape Canaveral arrival, "CC Enter"), 1-90 (Cape Canaveral departure, "CC Leave"), and 1-150 (Chesapeake Bay arrival, referred to as "CB Enter"). As in the full season models, the percentage of tree splits based on each habitat variable was interpreted as the relative influence of that variable during each of the migration periods and pairwise interactions were assessed.

RESULTS
We recorded 15,828 detections of the 51 tagged rays, which were summarized as 1298 daily presence/absence records. The tagged rays occurred within the Chesapeake Bay during the summer, arriving in early May and departing for the area off Cape Canaveral, Florida in late September or early October, with some individuals still detected as late as early November (Fig. 2). In the Chesapeake Bay, females were detected in all regions of the bay by early May ( Fig. 2A) and were detected in the MCB and UCB areas approximately two weeks before the males (Fig. 2B). Females remained in the MCB area of the bay throughout the summer and early fall ( Fig. 2A). The majority of males appeared to leave the Chesapeake Bay area by late July, and those still detected in the area were detected primarily in the LCB and ECB regions (Fig. 2B). Both females and males were detected off Cape Canaveral within the same approximate time period during the late fall, winter, and early spring, though more females were detected in early spring (Fig. 2C).
Presence was observed in SST as low as 14°C, but this was a rare occurrence and they were more commonly observed in SST above 15°C (Fig. 3A). The maximum SST observed with detections was 31°C. Tagged rays were found in areas with an average photoperiod of 13.3 h, with a maximum observed of 14.92 h and a minimum of 9.92 h (Fig. 3B). The temperature and photoperiod ranges within which tagged rays were detected were similar in both the Chesapeake Bay and the area around Cape Canaveral (Fig. 3).
Multicollinearity tests showed evidence of moderate pairwise collinearity. Pearson's correlations were moderate (0.4 < r > 0.7) for all variables in both the Chesapeake Bay and Cape Canaveral regions, but all correlations were significant (P < 0.001). The strongest correlations in both regions were between SST and photoperiod (Chesapeake Bay: r = 0.59, Cape Canaveral, r = 0.64). VIF values ranged from 2.41 to 3.54 across all explanatory variables. The results of both multicollinearity tests were below thresholds indicating strong collinearity (r > 0.7, VIF > 10) identified by Dormann et al. (2013).
The best-performing combination of parameters was the same for all BRT models (lr = 0.001 and bf = 0.6). Model performance was high overall, with CV scores over 0.6 and AUC values over 0.9 for all models (Table 1). SST was the most influential variable in all seasonal models and for all ray groups. Photoperiod was the second most influential variable in all seasonal models for the Chesapeake Bay. Photoperiod and day of year were nearly equal in influence for all rays combined and males in Cape Canaveral, and day of year was the second most influential variable for females in the overwintering habitat. The strongest pairwise interactions were between photoperiod and SST across all models (Table 1).
Marginal effects from seasonal presence models corroborated responses to SST and temporal variables observed in the tag detection data. For all rays combined, SST > 15°C, photoperiod > 11 h, and day of year between 120 and 300 showed positive effects on cownose ray presence likelihood in the Chesapeake Bay (Fig. 4A). These effects were similar for females, but the effect of day of year on presence likelihood in Chesapeake Bay showed a bimodal response with increased likelihood at days 110-180 and within a narrow range around day 300 (Fig. 4B). Similar marginal effect trends were observed among males with the exception of photoperiod in the Chesapeake Bay, in which a photoperiod > 13 h showed a positive effect on presence likelihood (Fig. 4C). In the Cape Canaveral region, presence likelihood for all rays combined declined as SST exceeded 20°C, between day of year 10 and 320, and as photoperiod exceeded 10 h (Fig. 5).
Boosted regression tree models of arrival and departure periods performed comparably to and male (blue) rays detected within areas associated with data buoys off Cape Canaveral by month. Chesapeake data buoy regions are Upper Chesapeake Bay (UCB), Middle Chesapeake Bay (MCB), Lower Chesapeake Bay (LCB), and East of Chesapeake Bay (ECB).
( Fig. 2. Continued) v www.esajournals.org   (Table 2). Model results showed clear differences in the relative variable influence between migration periods and between sexes during one of the migration periods. SST was most influential across all groups during departure from Chesapeake Bay and arrival at Cape Canaveral (Fig 6,  Table 2). Relative influence differed between sexes during departure from Cape Canaveral.
Day of year and SST were nearly equivalent in relative influence for all rays combined during this migration period, but SST was most influential for females and day of year was most influential for males. Day of year was the most influential habitat variable during arrival at Chesapeake Bay for all ray groups. Photoperiod was the second most influential variable during departure from Chesapeake Bay, but was the least influential variable during all other migration periods (Fig 6, Table 2). The strongest pairwise interaction was between day of year and SST for all ray groups during all migration periods except for departure from Chesapeake Bay for all rays combined and males, for which the strongest pairwise interaction was between photoperiod and SST (Table 2).

DISCUSSION
Departure of cownose rays from summer habitat in Chesapeake Bay and adjacent shelf waters was associated with declining SST across all groups, while northward migration from overwintering areas off Cape Canaveral, FL was associated with day of year for males and increasing SST for females. The differences observed for departure from overwintering habitat suggest male migration is driven more by endogenous cues than SST or by environmental cues not included in the analysis. These results suggest that migrations of Chesapeake Bay cownose rays fall into the "refuge" category from the framework proposed by Shaw (2016), whereby animals move away from one habitat with unfavorable environmental conditions and return to breed. Our findings are also consistent with sex-specific differences in migration timing, cues, and geographic extent in other fish species (Dahl et al. 2004, Sinnatamby et al. 2018, including other elasmobranchs (Barnett et al. 2011, Lea et al. 2015. Similar sex-specific migratory variations have also been found in studies on birds (Chandler andMulvihill 1990, Izhaki andMaitav 1998), which, like cownose rays, fit into the "refuge" migration category from Shaw (2016).
Our observation that the timing of southward migration for both sexes is triggered by falling temperatures may be related to behavioral thermoregulation to avoid cold winter temperatures. Previous research supports a lower lethal limit of 15°C for cownose rays (Schwartz 1964), which corresponds with the SST at which both males and females begin southward migration. Our data also suggest that southward migration likely varies across summer habitats and years based on the date at which local water SST falls Table 1. Boosted regression tree parameters (number of trees = n trees, learning rate (lr) = 0.001 and bag fraction (bf) = 0.6 for all models), performance values (cross-validation score = CV AE standard error, training data area under curve = AUC), and the percentage of tree splits attributed to each habitat variable (sea surface temperature = SST), variables included in the strongest pairwise interaction, and interaction strength for models of cownose ray (Rhinoptera boansus) seasonal presence by sex and region. v www.esajournals.org below this threshold. Ogburn et al. (2018) found latitudinal differences in the timing of southward migration of cownose rays from different summer habitats, with rays summering in more southern estuaries in Georgia initiating the southward migration later than rays from the Chesapeake Bay but all rays departing from the shared overwintering habitat at approximately the same time. The influence of temperature on fall migration timing is widespread among migratory fish species in the Northwest Atlantic Ocean. Henderson et al. (2017) found that fish distribution during fall surveys was significantly affected by temperature with most species occurring farther north during years with longer summer conditions. Interestingly, this timing co-occurs with a period of strong southward flowing currents along the southeast US (Weber and Blanton 1980) during late fall and winter, which likely reduces energetic costs of the migration. If temperature is the primary influence on the timing of southward migration and temperatures in summer habitats cease to fall below 15°C for the entire year due to climate change, this may remove the environmental cue for cownose rays to migrate south at all. There is already evidence that warming trends have affected migration phenology and nursery habitat selection in at least one elasmobranch species , and our results suggest this could also be possible for cownose rays.
Female departure from Cape Canaveral also seems to be associated with changing temperatures. The upper lethal limit for this species is unknown, but evidence suggests a preference for temperatures below 30°C (Smith and Merriner 1987, Collins et al. 2008, Ajemian and Powers 2012, Omori and Fisher 2017. However, female rays initiated northward migration at considerably lower temperatures, suggesting that rays are not leaving Cape Canaveral in direct avoidance of warm temperatures. Additionally, it is hypothesized that females prefer a narrower temperature range than males (Omori andFisher 2017, Ogburn et al. 2018) and prior studies suggest that migration into warmer waters by some elasmobranch species is tied to gestation (reviewed in Dudgeon et al. 2013). Many of the females that embark on this seasonal migration are pregnant and return to the Chesapeake Bay v www.esajournals.org to give birth (Fisher et al. 2013). Perhaps northward migration in females is delayed until temperatures are within the thermal range ideal for embryonic development. Temperature and other external conditions have been shown to have significant effects on embryonic development in other elasmobranchs (Wheeler et al. 2020), so maintenance of ideal conditions for development may be worthy of further study as a driver of migratory behavior.
Day of year was more influential than photoperiod or SST in predicting female arrival in Chesapeake Bay and for male departure from Cape Canaveral. This suggests that arrival in Chesapeake Bay by females and northward migration of males may be influenced by an endogenous timing mechanism, potentially related to the timing of parturition (females) or mating (males). Female cownose rays give birth in June or July and mate just a few weeks later (Fisher 2010, Fisher et al. 2013. In our study, both sexes arrived at the mouth of the Chesapeake Bay at approximately the same time in May, but female rays were abundant and consistently present throughout the upper portions of the bay approximately two weeks earlier than males. This movement into the upper bay coincides with the timing of parturition (Fisher et al. 2013). Cownose rays are not known to undertake a resting period between pregnancies (Poulakis 2013), and our data suggest males may travel into the upper regions of the Chesapeake Bay around the time of parturition by females, potentially improving access to mating opportunities. This behavior is in direct contrast to breedingrelated sex-specific migration behavior seen in birds, for which the earlier arrival of males at breeding grounds compared to females is much more prevalent (Pedersen et al. 2019). However, this is likely related to differences in physiology of reproduction between groups. The majority of males departed all but the lowest reaches of the bay by the end of July, shortly after the period of overlap with females in the upper bay, and their apparent lack of fidelity to the estuary during the remainder of the summer may maximize male foraging success in part by reducing competition with females and juveniles (Fisher 2010). In this way, male cownose rays may combine the "refuge" and "breeding" strategies set forth in Shaw (2016), since they have three distinct habitats. This has been seen in other fishes (Eiler and Bishop 2016), including other elasmobranchs (reviewed in Speed et al. 2010), but differs from migration strategies of migratory birds that typically follow a "refuge" migration strategy. The Table 2. Boosted regression tree parameters (number of trees = n trees, learning rate (lr) = 0.001 and bag fraction (bf) = 0.6 for all models), performance values (cross-validation score = CV AE standard error, area under curve = AUC), and the percentage of tree splits attributed to each habitat variable (sea surface temperature = SST), variables included in the strongest pairwise correlation, and interaction strength for models of cownose ray (Rhinoptera bonasus) arrival and departure by sex and region. importance of day of year is suggestive of consistent timing of annual northward migration. The ultimate cause of this consistent migration timing is likely a behavioral or physiological factor such as an annual cycle of hormones related to reproduction (Maruska and Gelschleichter 2011), though associations with the timing of prey availability have been suggested for at least one elasmobranch species (Barnett et al. 2011). Cownose rays show similar associations with temperature and photoperiod to those seen in other elasmobranchs, notably species that occupy a similar overwintering range. Kessel et al. (2014) found that the presence of lemon Fig. 6. Relative influence of each habitat variable (% of tree splits attributed to that variable in BRT models) during periods of departure from Chesapeake Bay ("CB Leave"), arrival at Cape Canaveral ("CC Enter"), departure from Cape Canaveral ("CC Leave"), and arrival at Chesapeake Bay ("CB Enter") for all tagged rays combined, females only, and males only. v www.esajournals.org 13 September 2021 v Volume 12(9) v Article e03743 sharks (Negaprion brevirostris) between Cape Canaveral and Delray Beach, FL, was most strongly associated with temperatures below 24°C, with photoperiod also having some effect. Kaijura and Tellman (2016) found that mass migrations of blacktip sharks (Carcharhinus limbatus) into the waters of Palm Beach County, FL, were associated with temperatures less than 25°C and also showed a significant correlation with day of year but not photoperiod. That multiple species of migratory elasmobranchs show similar associations with temperature off the Atlantic coast of Florida demonstrates the importance of the area as a thermal refuge during winter. The presence of migratory sharks in the area also correlates with the timing of spawning aggregations of prey fish species, suggesting that migrations may be timed to follow movements of prey (Kaijura and Tellman 2016). Cownose rays primarily feed on benthic invertebrates that are either non-migratory or completely immobile (Ajemian and Powers 2012, Bade et al. 2014) so following specific prey species is unlikely in this case. However, we cannot discount the possibility that migration timing may be associated with spawning behaviors or other seasonal changes in the abundance of benthic prey. Fishes and birds share many similarities in movement ecology that arise from traveling through fluid mediums (water for fishes, air for birds) rather than over land (Secor 2015). While migrations between summer and winter habitats are broadly shared between both groups, a key question is whether migration timing is driven primarily by an endogenous or exogenous mechanism and if comparing proximal migratory cues between taxonomic groups could allow us to better understand ultimate drivers of migration in the animal kingdom. Cownose rays in our study appear to fit into the "refuge" migration strategy presented by Shaw (2016), with males potentially displaying a combination of "refuge" and "breeding" strategies. Males and females migrate south to overwintering habitats in response to decreasing water temperatures and northward migration appears to be driven by an endogenous timing mechanism (especially for males). This is similar to the refuge migratory strategies of migratory birds. In addition, cownose rays may fit within the framework proposed by Winger et al. (2019) for birds, in which the ultimate driver of migration is seasonally variable environmental conditions in primary habitats. Based on prior studies that reported catches of individuals as small as 45-47 cm disk width, a potential primary nursery for cownose rays may exist in the Indian River Lagoon (Snelson andWilliams 1981, Roskar et al. 2020), which is adjacent to the overwintering habitat off Cape Canaveral used by rays from the Chesapeake Bay. Comparing migratory timing and extent between conspecifics born in the lagoon and other primary nurseries may be an effective way of testing this hypothesis for cownose rays. If this species fits within the framework of Winger et al. (2019), rays from the less environmentally variable nursery in the Indian River Lagoon may show significantly less migratory behavior than rays born in other estuaries.
Our analytical approach was similar to the historical approach for documenting bird presence at summer and winter habitats, although we used acoustic telemetry detections in place of visual observations. This approach was justified by the tendency of migratory cownose rays to only spend extended periods of residency within summer and winter habitats at the extreme ends of their migration even though other summer habitats exist where the timing of arrival and departure may differ (Ogburn et al. 2018). Therefore, our results should only be considered conclusive for rays summering in the Chesapeake Bay and used as a framework for comparison with other U.S. east coast estuaries serving as summer habitat. Our approach is easily replicable in other species and at other locations, and could be particularly powerful for comparisons across multiple species rather than the comparison between males and females conducted here. Our results cover a time period of less than four years and the influence of temperature on southern migration timing may introduce considerable interannual variability, especially with the effects of climate change on ocean temperatures (Henderson et al. 2017). Extending these observations to decadal scales represents a potentially powerful method for understanding the effects of climate change on coastal migratory species.
Collinearity between the explanatory variables is a potential caveat of our analysis. While multicollinearity test results fell below thresholds indicating strong collinearity (Dormann et al. 2013), all pairwise correlations between the explanatory variables were significant. However, the machine learning nature of BRT modeling makes it relatively insensitive to collinearity even when strong correlations are present (Elith et al. 2008, Dormann et al. 2013) Further, correlation and collinearity occur naturally between explanatory variables in ecological studies and strongly correlated variables can still be useful in analyses if the mechanistic differences in their relationships with the response variable are understood. Removal of strongly correlated groups of variables can also be detrimental to model performance (Dormann et al. 2013). Each of the three explanatory variables in our analysis was purposely chosen to compare external and potentially endogenous drivers. Temperature in particular can be highly variable between years, and the effects of climate change show that longterm changes in the relationship between temperature and day of year are already occurring. Photoperiod is much more consistent but can vary by location. In our study, photoperiod was approximately 20 min longer in the Cape Canaveral region than it would be in Chesapeake Bay during both winter and summer. There is precedent for including photoperiod independently of day of year in analyses of elasmobranch migratory cues (Kaijura and Tellman 2016). With day of year functioning as a proxy for endogenous behavioral cues, this selection of variables allowed us to assess the influence of a highly variable external cue (SST) and a relatively stable external cue (photoperiod).
Beyond contributions to theoretical studies on migration, our results have practical conservation implications. Increasing ocean temperatures resulting from climate change are likely to influence migration timing in cownose rays. Associations with temperature increase the likelihood of climate impacts on migration timing and extent among females in particular. For example, females may be more likely than males to occur within either the Chesapeake Bay or Cape Canaveral regions outside of the active period of a static time-area closure due to temperaturedriven changes in migratory behavior. Migration timing and environmental associations, including differences between sexes, should be taken into consideration as fishery management and conservation plans are developed for cownose rays.
The results of this study provide insight into the potential drivers and timing of cownose ray migration and contribute to the understanding of migration in general. Cownose ray migratory behavior exhibits characteristics seen in other taxa, including non-fish species, and fits within a variety of theoretical frameworks on migration ecology. While much of the debate concerning the evolutionary origins of migration has taken place in the context of avian research, cownose rays are clearly among the fish species with birdlike migratory behavior (Secor 2015), potentially making the species a viable non-avian model system for future comparative research.

ACKNOWLEDGMENTS
We thank Corey Corrick, Matthew Fisher, Michael Goodison, Margaret Kramer, Kristene Parsons, Cassidy Peterson, Robert Semmler, Christopher Sweetman, and several commercial watermen for technical assistance capturing and tagging rays. We also thank the coordinators of the ACT (Dewayne Fox and Lori Brown at the time the study was conducted) and FACT (Joy Young) networks and many researchers for sending tag detection data. The following acoustic arrays were particularly important in our analysis: the Cape Canaveral Array supported by the US Bureau of Ocean Energy Management (BOEM) and the NASA Environmental and Medical Contact #80KSC020D0023, lower Chesapeake Bay receiver arrays administered by Christian Hager at Chesapeake Scientific and supported by the US Navy and BOEM, mid-upper Chesapeake Bay receiver arrays administered by David Secor's laboratory at the University of Maryland Chesapeake Biological Laboratory and supported by the Atlantic States Marine Fisheries Commission and the NOAA Chesapeake Bay Office, and other receiver arrays within Chesapeake Bay administered by Charles Stence at the Maryland Department of Natural Resources and supported by an award from the NOAA Species Recovery Grants to States. We offer special thanks to Autumn-Lynn Harrison for providing useful feedback on this manuscript. Personnel and tagging efforts were supported by the Smithsonian Institution Office of the Undersecretary for Science, Aramco Services Company, and the Smithsonian Women's Committee.