Original Article Identifying the distribution of Atlantic cod spawning using multiple ﬁxed and glider-mounted acoustic technologies using

April 2019. Effective ﬁshery management measures to protect ﬁsh spawning aggregations require reliable information on the spatio-temporal distribution of spawning. Spawning closures have been part of a suite of ﬁshery management actions to rebuild the Gulf of Maine stock of Atlantic cod ( Gadus morhua ), but difﬁculties remain with managing rebuilding. The objective of this study was to identify the spatial and temporal distribution of cod spawning during winter in Massachusetts Bay to improve our understanding of cod spawning dynamics and inform ﬁsheries management. Spawning was investigated in collaboration with commercial ﬁshermen during three winter spawning seasons (October 2013– March 2016) using acoustic telemetry and passive acoustic monitoring equipment deployed in ﬁxed-station arrays and mounted on mobile autonomous gliders. Tagged cod exhibited spawning site ﬁdelity and spawning primarily occurred from early November through January with a mid-December peak and some inter-annual variability. The spatial distribution of spawning was generally consistent among years with multiple hotspots in areas > 50 m depth. Current closures encompass most of spawning, but important areas are recommended for potential modiﬁcations. Utilizing multiple complementary technologies and deployment strategies in collaboration with commercial ﬁshermen enabled a comprehensive description of spawning and provides a valuable model for future studies.


Introduction
The spatial and temporal distribution of the spawning of marine fishes plays a major role in their life history, recruitment success (Cushing, 1990), and interactions with fishing ( van Overzee and Rijnsdorp, 2015;Erisman et al., 2017). Spawning closures have been frequently implemented to prevent the overexploitation or disruption of fish while aggregated to spawn (e.g. Murawski et al., 2000;Russell et al., 2011). However, the effectiveness of spawning closures has varied, with some cases providing clear conservation benefits (e.g. Burton et al. 2005;Nemeth, 2005;Hamilton et al., 2011) and others yielding few or none (e.g. Clarke et al., 2015;Grüss and Robinson, 2015). The performance of spawning closures for meeting their objectives is influenced by multiple interacting factors, such as understanding of spawning behaviour, the size and timing of closures, fishing pressure in open locations and seasons, and the timeliness of management actions (Grüss et al., 2014). Therefore, reliable information on a species' spawning dynamics is critical for designing effective fishery management measures.
Atlantic cod (Gadus morhua) exhibit a lekking mating system (Nordeide and Folstad, 2000) and form large, dense spawning aggregations in locations and seasons that are often predictable (International Council for the Exploration of the Sea (ICES), 2005; Zemeckis et al., 2014a). Cod spawning sites are commonly located close to shore, making them easily accessible to fishing pressure (Armstrong et al., 2013), which can disrupt their complex mating behaviours (e.g. Morgan et al., 1997;Dean et al., 2012) and reduce the chances for successful reproduction (Brawn, 1969). These factors make Atlantic cod spawning components vulnerable to extirpation and prime candidates for the application of spawning closures as part of a multifaceted approach to fishery management (Grüss et al., 2014;Zemeckis et al., 2014a;Sadovy de Mitcheson, 2016). Fishery managers have implemented closures to protect spawning aggregations for many cod stocks (e.g. Baltic Sea, International Council for the Exploration of the Sea (ICES), 2004; Iceland, Jaworski et al., 2006), including in the Gulf of Maine (Armstrong et al., 2013) where efforts to eliminate overfishing and achieve rebuilding have not been effective (Rothschild et al., 2014;Northeast Fisheries Science Center, 2017). Declines in abundance of the Gulf of Maine cod stock have been associated with the extirpation of historical spawning components (Ames, 2004;Zemeckis et al., 2014b) and reductions in stock productivity and stability (Reich and DeAlteris, 2009;Kerr et al., 2014).
An improved understanding of the spatial and temporal distribution of spawning is required to inform fishery management and support rebuilding of the Gulf of Maine cod stock, which is comprised of genetically distinct spring-and winter-spawning subpopulations (Kovach et al., 2010;Zemeckis et al., 2014b). Previous studies have provided insights into the regional spawning dynamics of both subpopulations (e.g. Berrien and Sibunka, 1999;Hoffman et al., 2012), but research at finer spatial scales (i.e. spawning sites) has thus far focused on the spring-spawning subpopulation (Dean et al., 2012Gurshin et al., 2013;Siceloff and Howell, 2013;Zemeckis et al., 2014cZemeckis et al., , 2017. These studies demonstrated that spring-spawning Gulf of Maine cod exhibit multi-year spawning site fidelity, connectivity among inshore spawning sites, and complex sex-specific spawning behaviours that are vulnerable to disruption. Massachusetts Bay has long been a major winter-spawning ground (Fish, 1928;Rich, 1929), and this area has been the focus of multiple iterations of cod spawning closures since the 17th Century when harvest was prohibited in December and January (Charters and General Laws of the Colony and Province of Massachusetts Bay, 1814). From the mid-1990's until 2010, groundfish "Rolling Closures" prohibited commercial fishing in locations and seasons of high cod abundance, offering partial protection to spawning cod during October and November (Figure 1a; Murawski et al., 2005). In 2003, the Winter Cod Conservation Zone (WCCZ) was implemented within Massachusetts state waters (Figure 1a; Armstrong et al., 2013), and since 2005 has prohibited commercial and recreational fishing from 15 November through 31 January. The Gulf of Maine groundfishery transitioned to a catch share management system in 2010 and several effort control measures were removed, including the "Rolling Closures" and daily trip limits (New England Fishery Management Council (NEFMC), 2009). Therefore, spawning aggregations in federal waters were unprotected and this area was subjected to increased fishing pressure (Brewer, 2014). Given concerns about the sustainability of this expanding fishery, the group of commercial fishermen with whom we partnered to complete this project approached our research team to initiate a study that would investigate the spawning dynamics of this subpopulation and inform fishery management measures to prevent the extirpation of the remnant spawning components while still allowing fishing in adjacent areas not utilized by spawning cod.
Fishery managers implemented additional closures to protect spawning cod after this study began in 2013. For example, in 2014, emergency measures closed broad areas for November through February (Department of Commerce, 2014), which were modified in 2015 (i.e. Framework 53) to permit fishing for other more abundant species while still offering protection for spawning cod from November through January (Figure 1a; Department of Commerce, 2015). The critically depleted status of the resource created a sense of urgency and forced fishery managers to act on coarse spatial information and an incomplete description of cod spawning. The objective of this study was to identify the spatial and temporal distribution of cod spawning during winter in Massachusetts Bay to improve our understanding of cod spawning dynamics and inform the seasonal fishery closures intended to protect this group of spawning cod.

Methods
The distribution of cod spawning was investigated by combining acoustic telemetry data from individual cod that were tagged when in spawning condition with recordings of cod grunts produced by males as a part of their courtship rituals (Brawn, 1961;Rowe and Hutchings, 2006;Fudge and Rose, 2009) throughout Massachusetts Bay and on Stellwagen Bank (Figure 1b), including portions of the Stellwagen Bank National Marine Sanctuary (NOAA National Centers for Coastal Ocean Science (NCCOS), 2006). Combining data from both approaches provides a holistic description of the spatial and temporal distribution of cod spawning: acoustic telemetry identifies when and where tagged fish returning to the spawning ground aggregate, whereas the passive acoustic recordings provide indications of actual spawning events based on their mating system.
Putative spawning sites were monitored during three consecutive winter spawning seasons from October 2013 through March 2016. In order to establish a study area for the first project year, Identifying the distribution of Atlantic cod spawning multiple data sets were explored to identify putative spawning sites and seasons, including ichthyoplankton surveys (Berrien and Sibunka, 1999), an industry-based trawl survey , areas of high catch-per-unit-effort from commercial fishery observer data, reports from collaborating commercial fishermen, and existing passive acoustics data (van Parijs et al., 2015). Preliminary results from each year were used to annually modify the extent of the study area and timeframe of monitoring to best capture the distribution of spawning while balancing the need to maintain consistent stations and time periods across years.

Fixed station monitoring Acoustic telemetry
Spawning cod were tagged with acoustic transmitters during the first two project years (i.e. Year 1 ¼ winter of 2013-2014 and Year 2 ¼ winter of [2014][2015]. Cod were collected during 10-30 min tows aboard research charters with commercial bottom trawl fishing vessels. Total length (nearest cm) was measured for all fish and their sex and maturity stage were determined via visual inspection of eggs or milt extracted by cannulation or external pressure on the abdomen. Each fish was assigned a macroscopic maturity stage based on guidelines from Burnett et al. (1989). Cod were considered suitable for tagging with acoustic transmitters if they were in excellent physical condition, 55 cm total length, and either developing, ripe, ripe and running, or spent. Following the procedures of Dean et al. (2012Dean et al. ( , 2014, coded 69-kHz acoustic transmitters with a 60 s mean transmission interval (model V16-6H; 95 mm, 34 g in air, battery life > 1300 d: Vemco Division, AMIRIX Systems, Inc., Nova Scotia, Canada) were surgically implanted into spawning cod via small incisions (<4 cm) in the abdomen that were closed using braided silk sutures. The surgical procedures took <2 min and fish were monitored in holding tanks with fresh flowing seawater for up to 15 min to confirm full recovery from capture and the surgical procedure before release at the sea surface as close to the capture location as possible (<2 km).
Acoustic telemetry receiver deployments, maintenance, and retrievals were completed aboard commercial fishing vessels or Massachusetts Division of Marine Fisheries' research vessels. The deployment duration (i.e. receiver days) varied among stations and years ( Figure 1b, Table 1), because of lost receivers or vessel and equipment limitations. In Year 1, 34 Vemco VR2W acoustic receivers were deployed to monitor putative spawning sites from November 2013 through February 2014. The array was expanded to 42 receivers in Year 2 and the acoustic receivers were deployed from October 2014 to March 2015. The array was expanded to 56 receivers in Year 3, with the receivers deployed between mid-September and mid-October of 2015 and recovered March 2016. Fixed-station acoustic telemetry receivers were attached to mooring systems that positioned the receivers 18 m above the seafloor to effectively detect tag transmissions from spawning cod that typically remain within 3 m of the seafloor in Massachusetts Bay . Results from a range test in Massachusetts Bay indicated that the receivers had a detection radius of 1 km throughout a variety of weather conditions. The acoustic telemetry dataset was edited to remove tagged fish that were determined to have died based on a consistent lack of movement over the last month or more of detection. Given the detection area surrounding each receiver, it is difficult to determine the precise timing of death. Therefore, the entire final series of acoustic telemetry detections at the station where mortality was presumed to occur were omitted from analyses. Returning fish that were determined to have died mid-season were included in the calculation of mean arrival dates, but excluded from the calculation of departure dates and residence times. The sources of Array of fixed-station deployments for all three years of monitoring. Different shape symbols indicate station type and if an acoustic telemetry receiver was deployed at a station (circles) or a collocated acoustic receiver and MARU for passive acoustics (PA) monitoring (squares). Dots below each symbol indicate the years in which equipment was deployed at each station. The symbols are colour-coded based on assigned groupings. mortality could not be identified, but they could have included natural mortality, fishing mortality, or mortality associated with the capture and tagging process.
The spatial distribution of tagged fish was described using a Brownian bridge movement model. This technique was developed for datasets in which individual animal trajectories are comprised of a sequence of unique positions with known or estimable precision, such as GPS transmitters (e.g. Horne et al., 2007) or acoustic telemetry positioning systems (e.g. Dean et al., 2014). Our receivers were spaced far enough apart (2.8 km) to prevent simultaneous detection at multiple receivers. Therefore, the Brownian bridge movement model was applicable to our dataset, because any sequence of observations at different receivers represents movement between the receiver detection radii. The Brownian bridge movement model requires two parameters: mean location error (d) and Brownian motion variance (r 2 m ), which is related to animal mobility. The d parameter was estimated to be 740 m using the mean distance between estimated tag positions and a known receiver position (i.e. when within range of that receiver) from a previous study that employed a high-resolution positioning system (with identical tags and receivers to the present study) to track spring-spawning cod in Massachusetts Bay . The Brownian motion parameter was then estimated from the raw tag trajectory data via the maximum likelihood approach of Horne et al. (2007). This was done separately for males (r 2 m ¼ 1.90) and females (r 2 m ¼ 1.01), because the movement and space use of spawning cod is strongly influenced by sex and diel period . Separate estimates for day and night were not pursued, because the spatial resolution of the array was too coarse to reliably resolve sub-daily patterns.
After estimates for the two parameters were obtained, the Brownian bridge movement model was used to predict utilization distributions, which are two-dimensional probability distributions for the position of a tagged individual over a given time period (i.e. an entire spawning season, in this case). Utilization distributions were estimated only for fish that returned to the spawning ground in years after their release year, because this Table 1. Number of deployment days at each fixed station by year (i.e. Y1 ¼ Year 1) and receiver type (i.e. acoustic telemetry or passive acoustics), which are divided by region (see Figure 1).
permitted monitoring of movements throughout the entire spawning season from arrival until departure. Patterns in group space use were described by averaging the utilization distributions from individual fish to form a composite utilization distribution. Given that all tags were released in Years 1 and 2, and because the spatial extent of the receiver array expanded each year, the composite utilization distribution in Year 3 provides the most comprehensive description of spawning. However, to better evaluate inter-annual variability in space use, annual composite utilization distributions were also constructed using only those stations that remained consistent between Years 2 and 3.

Passive acoustic monitoring
Passive acoustic recorders were deployed to describe the spatial and temporal presence of spawning events based on the male courtship rituals as a part of the cod mating system. The grunting activity of spawning cod is distinguishable from the calls of other species in the study site and has been previously investigated in Massachusetts Bay (Hernandez et al., 2013;Stanley et al., 2017). There were 17 deployments of Marine Autonomous Recording Units (MARUs, Calupca et al., 2000), all of which were programmed to record continuously at a sampling rate of 2000 Hz and were moored so that they were 2 m above the seafloor in order to record grunts of cod that typically remain near the seafloor. The overall detection range of cod grunts on MARUs is unknown, but Stanley et al.  Table 1). All MARU's deployed in Year 2 and Year 3 were collocated with an acoustic telemetry receiver. MARU sound files were processed using a custom-built cod grunt detector (Urazghildiiev and Van Parijs, 2016) that was developed as a part of this project to expedite data processing and the analyses were executed in MATLAB (R2014b: MathWorks, Inc., Natick, MA). Detections automatically classified by the detector as cod grunts were manually verified using Raven Pro 1.5 (Bioacoustics Research Program, 2014) and the analyst removed false detections. Individual grunt detections were viewed in a 5x5 spectrogram grid adjacent to a context spectrogram. Detections in the grid were viewed from 10 to 400 Hz using a fast Fourier transform (FFT) of 256 points and 75% overlap, a 1 s time pad, and brightness of 55 and contrast of 74. The context spectrogram was viewed at a 10 s time window with a FFT of 1024 points and 75% overlap from 0 Hz to 500 Hz. Brightness was set at 61 and contrast at 60. Atlantic cod grunts in the Gulf of Maine have a fundamental peak frequency of 49.7 Hz and can have 2-8 harmonics depending on fish proximity to a MARU (Hernandez et al., 2013). Of the total number of potential grunts identified by the automatic detector (G a ), up to 2000 per station and recording day were individually examined by an analyst (G e ) and classified as either a true cod grunt (G v ) or not. In cases where G a exceeded 2000 per station/day, a random sample of 2000 automatic detections were examined and the number of true cod grunts (Ĝ v ) was estimated using the fraction of those examined that were verified as true using the following equation: Temporal patterns in cod grunting activity were examined by first summarizing all data from fixed stations by hour. Generalized linear models (GLMs) were then fit to both the presence of cod grunts and the number of cod grunts per hour (i.e. grunt rate). The grunt presence model assumed a binomial error distribution, while a zero-inflated negative binomial GLM was used for the grunt rate model. Best-fitting forms of both models were selected from a set of candidate predictors using Akaike's Information Criterion (AIC). Candidate predictors included: station, multiple natural cycles (diel, lunar, semi-lunar, annual), and an interaction between station and annual cycle. Circular variables were used to represent each natural cycle, which involved converting time to radians according to the period of the cycle and applying sine and cosine functions (Zar, 1999). Thus, each circular variable required two parameters and enabled estimation of the magnitude of the effect of the cycle and where in the cycle the peak response occurred. Both linear models were applied only to stations and years that had >5 grunts per day on average, because data from other stations were deemed too sparse to inform the model.

Autonomous glider surveys
Mobile Slocum autonomous gliders (Teledyne Webb Research Inc., Rudnick et al., 2004), which are being increasingly used to study the spatial ecology of marine fishes (e.g. Oliver et al., 2013Oliver et al., , 2017, were deployed in Years 2 and 3 to expand the coverage of the fixed-station deployments. All gliders were equipped with an externally mounted Vemco VR2W acoustic telemetry receiver, two digital acoustic monitoring instruments to record cod grunts (Johnson and Hurst, 2007;Baumgartner et al., 2013;continuous 2000 Hz sampling rate), and oceanographic sensors measuring conductivity, temperature, pressure, chlorophyll fluorescence, and turbidity (see Fratantoni, 2008 andBaumgartner et al., 2013 for further details on gliders).
In Year 2, two gliders (glider ID # we04 ¼ Year 2-1 deployment; glider ID # we10 ¼ Year 2-2 deployment) were simultaneously deployed in December 2014 at the same start location, but with a 12 h offset. These gliders were programmed to follow a track consisting of twelve parallel east-west transects encompassing a 1000 km 2 area that included the fixed-station array but covered a three to four times larger area. The 12 h offset was intended to survey the same areas during day and night to account for potential diel patterns in space use. In Year 3, three autonomous gliders (glider ID # we04 ¼ Year 3-1 deployment; glider ID # we10 ¼ Year 3-2 deployment; glider ID # we03 ¼ Year 3-3 deployment) were deployed serially to provide five consecutive surveys of Massachusetts Bay from 2 November 2015 through 1 March 2016. Each glider in Year 3 followed approximately the same track, consisting of ten parallel east-west transects and covering the same general areas as in Year 2. On average, a single glider survey covered 400 km, which resulted in an acoustic telemetry detection area of 800 km 2 assuming a detection range of 1 km.
The acoustic telemetry detection data from the glider surveys were analysed using a similar approach as the fixed-station data: detections from dead fish were removed and the Brownian bridge movement model was used to predict individual fish utilization distributions, which were then assembled to create composite utilization distributions by year. In addition, individual and composite utilization distributions were constructed using the combined fixed station and glider dataset.
It is possible that glider depth could influence the detectability of cod grunts given the limited effective detection radius of grunts (Stanley et al., 2017), the range of depths surveyed (Figure 1), and the roving nature of the gliders. Despite this, there was no apparent pattern to the grunt detections and glider depth. The passive acoustic recordings from the glider surveys were therefore analysed like the recordings from the fixed-station MARU deployments by using the custom-built cod grunt detector and reviewing them in Raven Pro 1.5. The glider data review deviated from the MARU data review only in brightness (context: 50, thumbnails: 47 or 40 [for we04 in 2015]) and contrast (context: 60, grid: 68 or 71 [for we04 in 2015]) values. Diel periodicity in cod grunting activity was examined with the glider-based recordings. The data were assigned to a diel period based on the sunrise and sunset times for Boston, MA, USA on each day of recording (day: between sunrise and sunset; night: between sunset and sunrise). The rate of cod grunting activity (estimated grunts per hour) was compared between day and night using a Wilcoxon rank-sum test given that the distribution of values did not conform to a normal distribution.

Acoustic telemetry
A total of 317 cod in spawning condition were tagged with acoustic transmitters (168 males, 149 females; mean ¼ 67 6 8 cm SD), including 155 cod in Year 1 (83 males, 72 females; mean ¼ 67 6 8 cm SD) and 162 cod in Year 2 (85 males, 77 females; mean ¼ 66 6 8 cm SD). A total of 23 tagged cod (7% of tagged fish) were determined to have died within the fixed-station array. Ten of these fish exhibited no signs of movement and were continuously detected at a single receiver for the duration of the study. The remaining 13 mortalities initially showed movement but were later detected continuously at a single receiver in the final month or more of the study. A total of 51% of the fish tagged in Year 1 were never detected during that year. Most of these undetected fish were tagged and released on Lower Stellwagen Bank (Figure 1b) on 10 January 2014, which was 15 km east of the fixed-station array and near the end of the spawning season. The expansion of the array and different tagging sites in Year 2 improved the detection efficiency with 17% of tagged fish in Year 2 being undetected.
Of the fish tagged in Year 1 that were not presumed dead, 20% were detected returning to the fixed-station array in subsequent years (18% Year 2; 10% Year 3). For the fish tagged in Year 2, 39% were detected returning to the fixed-station array in Year 3. The higher return rate between Year 2 and Year 3 was likely influenced by the larger array size and higher survival resulting from the implementation of emergency management measures in November 2014. The mean arrival date (i.e. date of first detection) of returning fish was similar between Years 2 and 3 (t-test: p ¼ 0.115; df ¼ 100) (Figure 2). However, some fish were already present when the receiver array was deployed in Year 2 (n ¼ 3) and Year 3 (n ¼ 1). When including data from the release years, the mean departure date (i.e. date of last detection) from the array was 2 weeks earlier in  (Figure 2). The mean residence time (i.e. departure date minus arrival date) of returning fish in Year 3 was more than double that of Year 2 (X Y2 ¼ 28 d; X Y3 ¼ 62 d; t-test: p < 0.001, df ¼ 100), which is likely related to the increase in array size in Year 3. Returning fish were present over the entire 6-month period, but 58 of the 75 tagged fish (77%) returning in Year 3 were detected outside of the current fishery closure period (i.e. November through January). However, only 12% of Year 3 detections occurred outside of the timing of the current Framework 53 closure, indicating that most of the spawning was encompassed by the closure.
In Year 1 and Year 2, the stations with the greatest number of detections were located near the centre of the array, south of the liquefied natural gas (LNG) terminals ( Figure 3). However, the stations that detected the greatest number of unique fish were located at the northeastern edge of the array in both of these years (Figure 3). There was a sixfold increase in the number of detections per fish with the expanded array in Year 3. However, tagged fish periodically went undetected while in the study area, with the typical fish having no detections on approximately half of the days during its residency period, which was likely related to incomplete coverage between acoustic receivers. During Year 2 and Year 3, there were consistent spatial and temporal patterns among the detections of fish within the fixed-station array (Figure 4).

Identifying the distribution of Atlantic cod spawning
Returning spawners arrived at the spawning ground primarily from two directions: from deeper waters to the south or across Lower Stellwagen. For most of the spawning season, the focal point for many fish was in the northern portion of the array, primarily among the "north deep" and "Lower Stellwagen" receiver groups, but the WCCZ was also important for some fish. The southern portion of the array was utilized less towards the end of the season, with most fish leaving the fixed-station array across Lower Stellwagen. In all years, only a small percentage of tagged fish and detections came from <50 m depth (Figure 4). There appears to be substantial inter-annual variability in space use, but this was most likely a function of the differences in array design across years because there was little difference between years when examining data from only those stations that remained consistent across years ( Figure 5).

Passive acoustic monitoring
Multiple MARUs failed to record due to hardware issues or prematurely detaching from their mooring, which resulted in variable deployment durations (Table 1). A total of 5594 out of 25429 automatic grunt detections (22%) on the MARUs were positively identified as cod grunts in Year 1. In Year 2, a total of 3142 of 45541 detections (7%) were positively identified as cod grunts. The increased proportion of false positives in Year 2 was largely due to the recording at Station 41 having a strumming noise from vibrations on the mooring system that caused the detector to falsely identify cod grunts. Year 3 had approximately four times the number of automatic grunt detections as the previous two years (n ¼ 165576). Consequently, a random subsampling of 2000 detections was required within several days (all at Station 64) prior to verification. Of the 106647 automatic grunt detections from Year 3 that were manually reviewed, 48999 (46%) were positively identified as cod. After accounting for subsampling, the total number of cod grunts in Year 3 was estimated to be 104956.   . Average utilization distributions (UD) of returning fish from acoustic telemetry. The 50th and 95th UD contours are identified, representing the minimum area inside of which there is a 50% and 95% probability of locating a tagged cod, respectively. The top row uses data from all fixed stations, the second row uses data from only those stations that were consistent between Years 2 and 3, the third row shows only those data collected from gliders, and the fourth row uses the entire dataset from all fixed stations and gliders combined. The 50 m isobath is shown in blue.

Identifying the distribution of Atlantic cod spawning
In Year 1 and Year 2, most of the grunts were recorded at Station 53 (Figure 6), with a maximum of 592 grunts in 1 d on 23 October 2013. In contrast, the maximum estimated number of grunts in 1 d at Station 64 in Year 3 was 16519 on 28 November 2015. While sub-sampling was required to make this estimate, 96% of the 2000 sub-sampled detections recorded on this day were verified as true cod grunts. Over the 9 d that were subsampled, 92-98% of the 2000 sub-sampled detections from each day were found to be true grunts. In total, 94% of the positively identified grunts recorded in Year 3 came from Station 64. The peak in cod grunting activity was consistently detected earlier in the fall at Station 53 (i.e. September-November) than at Station 64 in Year 3 (i.e. November-December) ( Figure 6).
The best-fitting model for grunt presence included all candidate variables, except for the semi-lunar cycle ( Table 2). The highest probability of grunt presence occurred during night approximately midway between sunset and midnight, between the full and waning moon, and from November through January (Figure 7). The one exception to this annual cycle was at Station 53, where the highest probability of grunts occurred 2 months earlier than elsewhere, hence the significant interaction between station and year. All predictor variables in the best-fitting model The models listed in bold font with the lowest AIC values were selected as the best fitting. S, station; D, diel cycle; L1, lunar cycle; L2, sem-lunar cycle; A, annual cycle.
had p-values of <0.0001. The best-fitting model for grunt rate had the same variables as the grunt presence model, but with the addition of the semi-lunar cycle. However, the magnitude of this semi-lunar effect was minor. In general, the effect of each natural cycle was stronger under the grunt rate model, with more defined peak times. The highest grunt rates occurred at 20:00 EST (i.e. 4-6 h after sunset), just after the full moon, and between late November and early December.

Autonomous gliders
The broad-scale coverage of the gliders proved effective at detecting the presence of tagged fish within the study site (Figure 8). In most cases, more tagged cod were detected by the gliders than the fixed array, despite the gliders yielding only 1-2% of the detections per fish recorded by the fixed array. The gliders also detected 19 of the 23 cod that were presumed to be dead based on analysis of fixed-station data. These fish were detected at similar locations to their last fixed array detections, thereby providing additional information supporting our mortality determinations. In total, 25 fish were detected by the gliders that were never recorded by the fixed array (n ¼ 13 from Year 1 releases; n ¼ 12 from Year 2 releases). Although some of these fish could be mortalities, the mobile gliders do not provide the continuous time series of observations necessary to infer mortality from a lack of movement. Across both telemetry datasets (i.e. fixed stations and gliders), 36% of putatively alive Year 1 fish were detected returning in subsequent years (29% Year 2, 27% Year 3) and 54% of Year 2 releases were detected returning in Year 3. The composite utilization distributions generated from glider telemetry data show a similar pattern in space use to that of the fixed array, but with a broader geographic extent ( Figure 5). Three areas of intensive space use were identified: along the northwest corner of Stellwagen Bank, the southeast part of and just outside of the WCCZ, and near the LNG terminals. Where and when they overlapped, the gliders and fixed telemetry array detected tagged cod in similar areas ( Figure 5). In Year 2, most glider acoustic telemetry detections occurred outside of the fixed array, particularly to the northeast. When the Year 3 fixed array was expanded in this direction, the spatial similarity between the glider and fixed array detections increased.
Most of the verified cod grunts recorded by the gliders were in areas frequented by tagged cod, except for an area just west of central Stellwagen Bank (Figure 8). A total of 125 of 19331 passive acoustic detections were verified as cod grunts across both glider surveys in Year 2 (Year 2-1, n ¼ 81; Year 2-2, n ¼ 44). The high frequency of false positives was likely due to rudder noise produced by the glider that triggered the detector. By examining each Year 3 glider survey, it appears that the number of true cod grunts correlates with the activity of tagged fish both spatially and Figure 7. Effects of circular variables on the presence and rate of grunting activity. For each radial plot in the top two rows, the portion of the line farthest from the centre indicates where in the cycle peak activity occurs. A perfectly centred circle would indicate a variable with no effect. For the plots in the bottom row, each line represents one of the stations where a marine autonomous recording unit (MARU) was deployed and they point towards the month(s) with the highest probability of grunt presence or the highest grunt rate.
temporally (Figure 8). In Year 3, 911 of 54052 automatic detections were verified as cod grunts across all five glider surveys (Year 3-1 ¼ 40; Year 3-2 ¼ 289; Year 3-3 ¼ 548; Year 3-4 ¼ 33; Year 3-5 ¼ 1). More than 60% of the verified grunts recorded by gliders in Year 3 came from the third survey (7-28 December 2015). This survey also saw the most widespread distribution of grunting activity, with grunts recorded throughout the study area, but most were concentrated between the northernmost transect and the western edge of Stellwagen Bank. The second glider survey in Year 3 (19 November-7 December 2015) also recorded many grunts (32% of Year 3 total), primarily near Lower Stellwagen and the northwest corner of Stellwagen Bank (Figure 8), which is a similar spatio-temporal pattern to the fixed-station passive acoustic data. No diel periodicity was found in the number of grunts recorded per hour by the gliders (Y2 Wilcoxon test: p ¼ 0.230; Y3 Wilcoxon test: p ¼ 0.071).

Discussion
This study identified the spatial and temporal distribution of cod spawning during winter in Massachusetts Bay by utilizing multiple acoustic technologies and the local ecological knowledge of collaborating commercial fishermen. The acoustic telemetry results indicated that there was some inter-annual variability in the timing of spawning, but most of the tagged fish were present from October through January during all 3 years with the greatest number of fish detected between early November and late January in Years 1 and 3. In contrast, most of the tagged fish left the array by the end of December in Year 2, which is believed to be due to colder water temperatures and prevalent winter storms. Passive acoustics data from both the fixed stations and glider deployments corroborated this timing. The observed temporal distribution of spawning was largely consistent among years and technologies, which is reassuring because acoustic telemetry indicates when and where fish are aggregating and the passive acoustic monitoring data provide long-term temporal information on the occurrence of spawning events. The agreement among different data types is within the typical monthly temporal scale at which fishery management measures are implemented in the Gulf of Maine.
The observed seasonality in spawning corresponds well with, but is not completely covered by, the current Framework 53 fishery closure period (November-January). For example, more than three quarters of returning tagged fish were detected before or after the closure. However, it is possible that some fish that were present for prolonged periods (i.e. >3 months) or outside of the peak spawning times were resident inshore pre-or post-spawning rather than actively spawning, and these are a very low proportion of the total number of detections. Therefore, our results suggest that the current Framework 53 fishery closure period is appropriately aligned with the timing of most of the cod spawning during winter in Massachusetts Bay.
The acoustic telemetry data identified multiple aggregation areas, including the area just west of the northwest corner of Stellwagen Bank, the southeast portion of the WCCZ and just east of the WCCZ's boundary, and the areas surrounding the LNG terminals ( Figure 1a). There was some inter-annual variability in the spatial distribution of spawning observed from the acoustic telemetry results (Figures 3-5), but this variability is within the typical spatial scale (i.e. 30 min squares) of fishery management in the Gulf of Maine. Both acoustic telemetry and passive acoustic monitoring results from Year 3 demonstrate that there was a hotspot of spawning along the northeastern portion of our study site, which is the area just northwest of Stellwagen Bank that is also within the western boundary of the Stellwagen Bank National Marine Sanctuary. It is possible that fishing activity along the boundary of the Framework 53 closure disrupted spawning and influenced the fish to aggregate just inside of the closure. However, fishing effort was severely reduced due to the low allocations of Gulf of Maine cod, so it is unlikely that fishing influenced the distribution of spawning during Year 3. The gliders transected an area 3-4 times larger than the fixed-station array and identified two areas of spawning beyond the fixed-station deployments, including near the northwest corner of Stellwagen Bank and along the glider transect just north of the current Framework 53 closure (Figures 5 and 8).
Only intermittent detections were recorded south of 42 15 0 N ( Figure 5), thereby indicating that this area was not utilized by many spawning cod. In future modifications to the system of closed areas, this area could be considered for re-opening, if further evidence confirms a lack of cod spawning. Conversely, fishery managers should consider including the area surrounding the northwest corner of Stellwagen Bank in the seasonal closed area. This location was identified as important for winter-spawning cod in Massachusetts Bay based on the acoustic telemetry data from fixed stations and gliders ( Figure 5), as well as passive acoustic monitoring data (Figure 8).
While the present study focused on describing the overall spatial and seasonal distribution of cod spawning aggregations along the Massachusetts coast in winter, there are ample details yet to be explored within the acoustic telemetry dataset. Future analyses could further describe finer-scale spawning behaviour and patterns among individual fish to investigate whether size or sex influence space use and behaviour, thereby building on previous research (e.g. Dean et al., 2014). Additional monitoring of the area should be continued in order to evaluate whether the distribution of spawning observed during this study remains consistent and that these areas warrant being closed. Exploitation of this subpopulation is still possible as fish migrate to or from their inshore spawning sites, because the consistency in the spatial and temporal patterns observed during this study can make these fish a predictable target at the margins of the closure. Contrary to historical accounts and some of the information available at the beginning of this study, we found little to no occupancy by tagged cod along the west side of the fixed-station array and in waters <50 m depth. It is possible that local environmental changes (Pershing et al., 2015) influenced a shift to spawning in deeper waters, that the spawning components which historically spawned in these shallower areas have been extirpated, or cod that represent separate, untagged spawning components may still be aggregating to spawn closer to shore.
Many of the acoustically tagged cod (up to 54%) exhibited spawning site fidelity by returning to the array over multiple consecutive winter spawning seasons. This behaviour by cod has been previously documented within other stocks (e.g. Robichaud and Rose, 2001;Skjaeraasen et al., 2011) and spring-spawning cod in Massachusetts Bay (Zemeckis et al., 2014c), but this is the first evidence of spawning site fidelity among winter-spawning cod in Massachusetts Bay. The mean residence time of fish in Year 2 (x¼28 d) was similar to that of spring-spawning cod in Massachusetts Bay (x¼38 d: Zemeckis et al., 2014c) and Ipswich Bay (x¼ 30 d: Siceloff and Howell, 2013). However, the mean residence time of fish in Year 3 (x¼ 62 d) was approximately twice as long as these other observations, which was likely due to the larger spatial scale over which fish were monitored in Year 3.
The MARU deployment at the southernmost site (Station 53, Figure 1b) reliably detected cod grunts during the 3-year study period ( Figure 6), albeit several weeks earlier than all other locations. This was historically an important cod fishing ground known as "Fishing Ledge" and spawning cod were caught there during previous trawl surveys . However, few tagged cod were detected in this area (Figure 3), suggesting that the grunts recorded at this location could have been from a separate and untagged spawning component, or that the detected grunts could have been agonistic or served another communicative function. Other than at Station 53, the spatial and temporal patterns between acoustic telemetry detections and cod grunts were similar.
The largest number of cod grunts recorded was at Station 64 in Year 3. The number of grunts and their high overlap in time could loosely be termed a "chorus" (Cato, 1978;Nordeide and Kjellsby, 1999). However, chorusing was not directly assessed here, because Massachusetts Bay is an area of high anthropogenic activity and ambient noise is affected by passing vessels (Stanley et al., 2017), which make standard mechanisms for measuring chorusing difficult to apply. Cod calling activity was associated with multiple variables, whereas, the highest probability of grunt presence occurred at night between sunset and midnight, during the full and waning moon, and from November through January (Figure 7). Therefore, our results are consistent with Rowe and Hutchings (2006) who studied cod in holding tanks without interference from ambient noise and also found that grunt presence was highest at night. Nonetheless, the detectability of cod grunts in the wild could be influenced by diel variation in ambient noise, but a cursory review of ambient noise and vessel traffic data did not show a clear diel relationship in Massachusetts Bay. Future research should further examine the interactions among grunt detectability, ambient noise levels, and vessel traffic to assist with the interpretations of behavioural information from passive acoustic monitoring data. The highest grunt presence was observed between the full and waning moon, thereby suggesting a connection between cod spawning and the lunar cycle. The lunar cycle has been documented as a driver of fish calling activity across many species and geographic regions (e.g. Parsons et al., 2016;Monczak et al., 2017), including Atlantic cod off Iceland where spawning was tied to the lunar cycle, but there was no evidence for a connection with diel cycles (Grabowski et al., 2015).
The Gulf of Maine cod stock has continued to decline in abundance (Northeast Fisheries Science Center, 2017) despite all of the management measures to protect spawning cod. It is possible that these measures have not been successful, because they were "too little and too late" (e.g. Clarke et al., 2015), or due to the inability to end overfishing in other locations and seasons (Rothschild et al., 2014;Northeast Fisheries Science Center, 2017), because it has been suggested that spawning closures are most likely to be effective as a multifaceted approach that concurrently reduces overall fishing mortality (Grüss et al., 2014;Zemeckis et al., 2014a;Sadovy de Mitcheson, 2016). Nonetheless, the spawning protection measures for Gulf of Maine cod are also intended to reduce the likelihood of extirpating semi-discrete spawning components and disrupting natural spawning behaviours (Dean et al, 2012). Therefore, although stock rebuilding has not yet been recognized, it is likely that these spawning closures have been successful at achieving these objectives. These are important factors to consider, particularly under a catch share management system in the Gulf of Maine based on annual catch limits.
The impetus for this study came from local commercial fishermen who predicted the need for more robust scientific information to inform fishery management. Their local ecological knowledge helped to design and execute this project which employed multiple acoustic technologies and deployment strategies, each of which have their respective advantages and disadvantages. Acoustic telemetry data are only received from tagged fish when within the detection range of acoustic receivers and establishing and maintaining a large array can be challenging. Furthermore, the loss of acoustic receivers prevents the recovery of valuable data and equipment. However, fixed-station acoustic telemetry receiver deployments provide very valuable long-term data from known locations. Passive acoustic monitoring of sound produced by fish allows information to be received from any individual in the population and provides behavioural context, but it is believed that male cod produce the most sounds. When applied to cod, the low source level of their grunts in a noisy ocean (see Stanley et al. 2017) means that moored receivers must be very close to the spawning. The moored MARU's were all recovered, but 25% (4 of 16) failed to capture data for a full season, which makes equipment reliability another limitation. Autonomous gliders can carry both acoustic telemetry receivers and passive acoustic recorders, which enables the monitoring of animals over a broader area than typically possible by maintaining arrays of equipment moored at fixed locations (e.g. Wall et al., 2012;Mann et al., 2014;Schofield et al., 2015). The gliders performed well (i.e. were able to follow most of their predetermined tracks), are less susceptible to loss, and provide broad spatial coverage. However, the lack of continuous sampling at any one location makes gliders ineffective for identifying mortalities of tagged fish. Few studies have undertaken a similar multitechnology approach by utilizing acoustic telemetry and passive acoustic monitoring to investigate fish spawning dynamics (e.g. Rowell et al., 2015), and the incorporation of fixed stations and glider surveys with each technology was a valuable advancement that compensated for the advantages and disadvantages of each technology.
In conclusion, this study involved the collaboration of local commercial fishermen and utilized multiple acoustic technologies to identify the spatial and temporal distribution of cod spawning during winter in Massachusetts Bay. Based on a combined synthesis of the acoustic telemetry and passive acoustic monitoring data, from both fixed stations and glider surveys, the temporal distribution of cod spawning during winter in Massachusetts Bay was shown to have some inter-annual variability but spawning primarily occurred during early November through January with a peak in mid-December. The spatial distribution of spawning was generally consistent among years and concentrated in areas >50 m depth. Our results documented multiple hotspots of spawning, including just west of the northwest corner of Stellwagen Bank as the primary focal point of spawning, with other lesser focal points inside and near the WCCZ and LNG terminals. Our results are largely consistent with the spatial and temporal extent of the current Framework 53 closures. However, important areas are recommended for further evaluation by fishery scientists and managers to potentially modify current closures to optimally balance the protection of spawning cod, scientific uncertainty, and access to other more abundant species. The utilization of multiple complementary acoustic technologies and deployment strategies in conjunction with the local ecological knowledge of commercial fishermen proved critical for providing a comprehensive description of cod spawning dynamics and serves as a useful model for future studies with cod and other species.