Long-term passive acoustic recordings track the changing distribution of North Atlantic right whales (Eubalaena glacialis) from 2004 to 2014

Given new distribution patterns of the endangered North Atlantic right whale (NARW; Eubalaena glacialis) population in recent years, an improved understanding of spatio-temporal movements are imperative for the conservation of this species. While so far visual data have provided most information on NARW movements, passive acoustic monitoring (PAM) was used in this study in order to better capture year-round NARW presence. This project used PAM data from 2004 to 2014 collected by 19 organizations throughout the western North Atlantic Ocean. Overall, data from 324 recorders (35,600 days) were processed and analyzed using a classification and detection system. Results highlight almost year-round habitat use of the western North Atlantic Ocean, with a decrease in detections in waters off Cape Hatteras, North Carolina in summer and fall. Data collected post 2010 showed an increased NARW presence in the mid-Atlantic region and a simultaneous decrease in the northern Gulf of Maine. In addition, NARWs were widely distributed across most regions throughout winter months. This study demonstrates that a large-scale analysis of PAM data provides significant value to understanding and tracking shifts in large whale movements over long time scales.


Results
A total of 35,600 days of acoustic recordings were processed with the Low Frequency Detection and Classification System (LFDCS) 43 , and subsequent NARW upcall detections were manually reviewed, resulting in 2,527 days (7%) with confirmed NARW acoustic presence. NARWs were acoustically present along the entire eastern seaboard of North America from the western Scotian Shelf (region 3) to the waters off Jacksonville, Florida (region 10) throughout the winter months from late October through early April, with the exception of no detections on the Scotian Shelf from December through February (Fig. 1). A decrease in detections was seen in summer months in southern regions, reflecting the known movement of breeding individuals towards northern feeding grounds. NARWs were also detected near Iceland and Greenland (region 2) from July-October (see previous analysis by Mellinger, et al. 22 ). Davis Strait (region 1) contained one day with possible NARW detections in both December Scientific RepoRts | 7: 13460 | DOI: 10.1038/s41598-017-13359-3 and March. However, with high bowhead whale (Balaena mysticetus) and humpback whale calling also present, NARW presence could not be confirmed with confidence. There were no NARW detections in any recordings from Bermuda or the Caribbean (region 11), confirming the likelihood that NARWs do not currently venture into the waters surrounding Bermuda and the Caribbean.
Average weekly acoustic presence was broken up into two time periods representing before (2004-2010) and after (2011-2014) the described distribution shift starting in 2010 ( Fig. 2) 33 . To test whether the occurrence of right whales in regions differed over the two time periods, we ran a Generalized Linear Model (GLM; see methods) with the number of days in which whale calls were detected as the dependent variable, and the time periods (2004-2010; 2011-2014) and regions as independent variables, with their interaction effects included in the model. There were too few data (in some time*region cells) available for Davis Strait, Iceland and Greenland, Georges Bank, Cape Hatteras, and Bermuda and the Caribbean (regions 1, 2, 6, 9 and 11) for them to be included in the model.
For all other regions, both factors and their interactions were significant (Table 1). Pairwise comparisons of time periods across individual regions (run using phia::testInteractions) demonstrated differences between the two time periods (Table 2), except for in Massachusetts Bay (region 5). Northern regions (3 and 4; Scotian Shelf and Gulf of Maine) saw a reduction in calls in 2011-2014, but mid-Atlantic regions (7 and 8; southern New England and the mid-Atlantic), and the Southeastern U.S. (region 10) saw increases (Table 3, produced using phia:interactionMeans). We used the False Discovery Rate 44 to adjust for alpha-value inflation.
Spatial distribution of NARW acoustic occurrence was summarized from 2004-2014 by seasons in Fig. 3. Seasons were defined based on Roberts, et al. 45 with November to February as Winter, March to April as Spring, May to July as Summer, and August to October as Fall. Seasonal acoustic occurrence of NARWs reflected patterns seen in the daily presence plots (Figs 1 and 2), with NARW presence along the entire coast in both winter and spring seasons. Across all seasons, NARWs were detected from Cape Hatteras (region 9) to Nova Scotia (region 3), highlighting the expansive habitat of NARWs for most of the year (Fig. 3).
The LFDCS, with a Mahalanobis Distance (MD) threshold of 3.0, missed days with NARW upcalls at an estimated rate of 31%; if we used a threshold of one true upcall detection per day for daily presence, the rate of missed days with NARW upcalls would have been 25%. Manually tallied upcalls from the Gulf of Maine analysis 21 revealed that days where LFDCS detected true NARW upcalls had a median of 259 calls per day and ranged from 20-2770 calls per day, and days where LFDCS missed NARW upcalls had a median of 7 calls per day and ranged from 1-66 calls per day. It is important to note that detection rates, and therefore missed detection rates, will be highly variable due to variability in ambient noise levels, recorder types, habitat, bathymetry, and calling behavior  Fig. 4 and for all years of the study (2004)(2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014). Horizontal lines within the boxes indicate the median, box boundaries indicate the 25 th (lower boundary) and 75 th (upper boundary) percentiles, vertical lines indicate minimum and maximum values, and black dots represent outliers. Grey blocks indicate weeks where no data were available for that region. of the animals. However, this analysis indicates that at greater calling rates, indicative of higher calling activity and potentially reflective of multiple calling individuals, daily NARW acoustic presence was most likely captured in areas where they were vocalizing. Sporadic calls, which may be produced by lone individuals passing through an area, were less likely to be detected in this analysis (note that lone animals can be missed by visual surveys, particularly aerial surveys, as well). Given the main goals of capturing broad-scale movements of the entire population with this analysis, failing to detect the presence of some individuals does not compromise the overall results.

Discussion
This study on acoustic presence of NARWs demonstrates this species' high mobility and broad geographic range. It also shows that the habitats NARWs frequent can change over time and that they are often distributed widely across their entire range. Previous understanding of NARW movements assumed the majority of the population migrated between the calving grounds in winter months and northern feeding grounds in summer months, where visual survey effort was therefore concentrated 27,46 . Little was known about the whereabouts of whales that did not frequent these habitats, or the exact timing and location of other important areas. This study demonstrates nearly continuous year-round NARW presence across their entire habitat range, particularly north of Cape Hatteras, suggesting that not all of the population undergoes a consistent annual migration. It is possible that the non-migrating whales could be mobile individuals occupying broader areas throughout the year, like some subspecies of blue whales in the Antarctic 47 , or individuals that do not migrate annually, like a portion of the east Australian population of humpback whales 5 . Our data clearly demonstrate that NARWs occur along the entire eastern seaboard of North America for most of the year, even if that distribution has shifted within the past decade.
The acoustic data supports NARW distributional changes that have been recently observed during visual surveys. Furthermore, the acoustic data provide additional insights into where NARWs are located at times of the year when the poor weather and lack of light make visual surveys highly restricted (i.e., late fall to early spring). NARWs appear to have shifted from previously prevalent northern grounds, such as the Bay of Fundy and greater Gulf of Maine (regions 3 and 4), to spending more time in mid-Atlantic regions year-round (regions 7 and 8) 25,48 . In addition, improved visual survey effort in Canadian waters shows that NARWs have also been sighted more frequently further north in areas such as the Gulf of Saint Lawrence 35 . Acoustic data are currently being processed to allow further evaluation of this region, but were not available for this study. This is an area of extreme importance, with at least 10 linked NARW deaths this year alone 49 . Moving forward, it is crucial to have better available coverage for past and future monitoring efforts in such important regions where entanglement and ship strike threats are prominent. Likewise, we know of past NARW occurrence in historical habitat areas such as Greenland and Iceland, with acoustic presence in 2007 22 , and it would be valuable to understand current usage of this region and confirm the northern extent of their range with more recent data. Therefore, this type of long-term, large-scale approach using PAM is invaluable for helping to track presence, and changes in presence, of key species to understand (1) how they have shifted their distribution and (2) whether these shifts puts them at increased risk from anthropogenic threats. It remains unclear if these observed distribution shifts are due to environmental or anthropogenic effects, if they are a response to short-term changes in the environment, or part of a longer-term cycle in which NARWs shift their distribution. With recent studies finding the Gulf of Maine is the fastest warming body of water in the world 50 , it is not surprising to see distributional changes across marine species. We suspect further changes in distributions will occur as water temperatures continue to rise, forcing movements towards both favorable oceanographic conditions and food sources elsewhere. Regardless of the factors influencing these changes in distribution, it is critical for management strategies to reflect new threats that may arise for this species as they move into regions outside of existing management areas 51 .
The purpose of this study was to provide baseline information of NARW distributions across their current range. Here, PAM data was used for distributional analysis and assumes homogeneity in the probability of upcalls across all regions. We caution that this is not necessarily the case, as known calling rates vary among demographic groups, such as quieter mother-calf pairs 52 , resulting in lower detection probabilities. Call rates are also lower with certain behaviors, such as foraging and logging 53 , and therefore it is reasonable to suspect different calling rates over the different habitats and seasons summarized in this study. Additional variation exists across all the regions, with detection probabilities likely varying with different acoustic habitats, bathymetry and recorders.       Table 1 shows the overall GLM results between the two time periods (timePeriod) and region (Region). (***) indicates significance of P < 0.001.
Risch, et al. 12 show that ambient noise levels varied across three sites included in this study (Massachusetts Bay, mid-Atlantic, and Southeastern U.S.; regions 5, 7, and 10) over each season, affecting detection probabilities and ranges; while Rice, et al. 54 provide detailed results on the varying acoustic environment at 10 different sites within the NARW range. These studies highlight the challenges with large-scale datasets, however, further analysis into detection ranges were beyond the scope and goal of this study. Thus, this work used an average range of detection distances based on published studies 55,56 . As such, it provides a minimum estimate of NARW presence, with the understanding that in some regions the detection ranges may have been slightly more constricted or wider than the average used. This study integrates data from a suite of smaller-scale studies focused on the fine-scale occurrence of NARWs and other species. Comparing our results with these smaller-scale studies such as Hodge, et al. 48 , Bort, et al. 21 , Salisbury, et al. 25 , Mellinger, et al. 22 , Morano, et al. 23 , and Whitt, et al. 57 , there are time periods with discrepancies between the NARW acoustic detections presented in our study and these previous studies, with the previous studies finding right whale calls in some months that we did not find confirmed detections. This is likely due to a difference in detectors and analysis methodology. In order for this study to include and process such a large acoustic dataset, it was necessary to use an automated detector, LFDCS, and manually review these detections at a coarser level than may have been done in the smaller-scale studies. Our review of the LFDCS performance with a MD threshold of 3.0 at three sites found that 31% of days with right whales present were missed when compared to a full manual evaluation of the data, and that days with lower calling rates tended to be missed more often. Therefore, our results represent a minimum presence compared with these more detailed studies, and it is possible that NARWs have higher occurrence in some areas than is reported here. The missed detection rate could be reduced by lowering our classification threshold (i.e., increasing the maximum MD), but at the cost of increasing the time required to manually screen each detection, since increasing the maximum MD will result in more detections that must be reviewed by an analyst. As with all detection systems, there is a balance between accuracy and processing time that must be considered when choosing a detector within the context and scope of any study. This is one of the first comprehensive, long-term passive acoustic studies to investigate an entire habitat range for a marine mammal at this temporal and spatial scale, and is made possible only by the cooperation and collaboration of an extensive research community. Even in areas where data were collected for alternate purposes, the combined contributed recordings provided crucial information to assess both the acoustic occurrence and changing distribution of the NARW population. This analysis demonstrates what can be accomplished for other poorly understood species, and encourages broad research collaborations in the future. All contributing data sets were combined for the common goal of understanding the distribution of a critically endangered species facing extreme threats from anthropogenic and environmental influences. In planning future large-scale studies, standardization of acoustic recorders and methods should be considered to improve the quality of datasets.
PAM is a powerful, cost-effective, long-term monitoring tool that can give a better understanding of temporal trends and reveal range expansion, decline, or distribution shifts in populations, as well as interannual changes. This information can be used to direct science and management to focal areas of interest. Most importantly, in an ocean where conditions are changing rapidly, adaptive management is needed to identify and protect areas that are crucial for species on the brink of extinction. Potential ways forward include setting up real-time passive acoustic monitoring systems (see NEPAN 39 ), or thinking beyond the traditional means of classifying NARW critical habitat as static, confined areas. This is especially relevant when considering the mobile nature of this species whose distribution patterns may still be changing. It is imperative to continue effective surveys and timely conservation efforts to ensure the recovery of this endangered species.

Methods
Data collection. Passive acoustic data collected for a multitude of different research projects and goals were combined to examine a decade-long acoustic record of the spatial and temporal occurrence of NARWs throughout the western North Atlantic Ocean. This large area was divided into 11 regions, covering the main historical areas where the existing NARW population has been found since the late 1600s, ranging from Florida, USA to Cape Farewell, Greenland 58 . Recorders were assigned a region based on the biological and geographical importance of the area (Fig. 4). Region 3 was broken up to include subregion 3A, representing the Gulf of St. Lawrence: an area having an increase in NARW sightings and research effort over the last few years, but insufficient acoustic data available to contribute to this study. Data were examined across two time periods (2004-2010 and 2011-2014) representing before and after the observed distribution shift in 2010. Due to the ad hoc nature of this largescale collaborative project, available data are patchy, a few recorders were duty cycled, and there are some regions or time-periods with no available data (Fig. 5).
Five types of bottom-mounted passive acoustic recorders were deployed from 2004 through 2014 (See Supplementary Table 1): the High-frequency Acoustic Recording Package (HARP) 59 , the Marine Autonomous Recording Unit (MARU) 55 , the Autonomous Multichannel Acoustic Recorder (AMAR) 60 , the National Oceanic and Atmospheric Administration's (NOAA) Pacific Marine Environmental Laboratory's (PMEL) Moored Autonomous Hydrophones (HARU) 61 , and the Guardbuoy (http://geospectrum.ca/guard-buoy). Data collected from each of these recorders varied from a minimum of 25 days to a maximum of 2 years (Supplementary Table 1). Some recorders (59 out of 324) were duty cycled, ranging in recording from 12-95% of the time, while most (265 out of 324) recorded continuously. The spatial configuration in which they were deployed varied from single units to lines and arrays of recorders; the configuration was determined by the original goal of the specific research project in question (Fig. 4). The majority of recordings were sampled at 2 kHz, with some ranging up to 250 kHz. All recordings were low-pass filtered and decimated to 2 kHz for analytical consistency across data in order to make them comparable.
Maximum detection ranges for NARWs can vary considerably, depending on recording equipment, location, and environmental conditions, as well as call type and behavioral context 62 , but are estimated to range from 8 to 16 km 55,56 . Consequently, single recorders were selected for analysis from array configurations with units spaced less than 8 km apart, to ensure full coverage of the area while minimizing duplicative detections across recorders. We focused our analyses on data collected between January 2006 to December 2014, with the exception of data collected in 2004 and 2005 in the Bay of Fundy, Emerald Basin, and Roseway Basin, Canada, since these were the only long term recordings available for these areas that had previous well-known occurrence of NARWs. Data from a total of 324 recorders were analyzed, comprising 35,600 recording days of data. All government funded acoustic data are publicly available upon request from the data owner.
Detection and classification of NARW calls. All acoustic data were processed using the Low Frequency Detection and Classification System (LFDCS) 43 which creates conditioned spectrograms (Fig. 6) using the short-time Fourier transform with a data frame of 512 samples and 75% overlap resulting in a time step of 64 ms and frequency resolution of 3.9 Hz. After tracing contour lines, or "pitch tracks", through tonal sounds, the program uses multivariate discriminant analysis to classify the pitch tracks into call types. Calls were classified based on a user-developed call library; our library included four North Atlantic baleen whale species: NARW, fin (Balaenoptera physalus), sei (B. borealis), and humpback whales. Here, we focused only on the detections classified as NARW calls, specifically the low-frequency modulated upsweep known as the upcall. The upcall is a contact call used throughout the NARW range, produced by all ages and both sex classes, and is therefore the most reliable call to use for determining right whale presence 53,55 .
The call library described in Baumgartner and Mussoline 43 was expanded and improved for this analysis to include a wider variety of examples of NARW upcalls and increase detection probability. Each detection was assigned a MD, which measures the deviation of a detection from the assigned call type (see Baumgartner and Mussoline 43 for a more complete description). A lower MD indicates a closer match to the assigned call type. All NARW upcall detections with a MD less than or equal to 3.0 (after Baumgartner, et al. 63 ) were manually screened by experienced analysts to determine which were correctly classified. For an ideal call type in the LFDCS (i.e., the seven attributes used in the discriminant function analysis are multivariate normal), 75% of actual calls will have a MD of 3.0 or less 63 ; we chose this threshold to make the laborious process of manual screening manageable at the expense of sometimes missing genuine right whale upcalls (see below). This approach ensured that false detections were eliminated. The high degree of variability in NARW upcalls and the overlap with other species' vocalizations, such as upsweeps produced by humpback whales, necessitated this extra manual step in data processing 63 .
For continuous data, a given day was marked as having NARWs present if three or more true upcall detections were found. Three upcalls were used to establish presence in order to be conservative and confident in stating NARW presence (we also conducted all analyses using a criterion of one upcall per day to indicate daily right whale presence, but neither our results nor conclusions changed). For duty-cycled data, the criteria was dropped to one true upcall detection signifying NARW presence so that presence was not underestimated due to a lower probability of recording vocalizing animals 64 . Weekly NARW acoustic presence per recorder was then summarized as the number of days per calendar week with daily presence.
To test whether the occurrence of right whales in regions differed over the two time periods, we ran a Generalized Linear Model (GLM) in R 3.4.1 65 , using the libraries ggplot2 66 , MASS 67 , car 68 , and phia 69 . This had the number of days in which whale calls were detected as the dependent variable, and the time periods (2004-2010; 2011-2014) and regions as independent variables, with their interaction effects included in the model. As the call data were counts, we ran the GLM with a Poisson distribution with log-link. The number of recording days was multiplied by the duty-cycle to correct for non-continuous data. Because effort (the number of days during which recorders were present) varied across time and region, we included the log of the number of days during which recorders were present plus 1 (as for some time*region cells, there were no recorders present) as an offset in the model.
The model formula in R was: * = = + nDaysWithWhales timePeriod Region, family 'poisson', offset log(nDaysRecordings 1) Detector evaluation/missed detection rate. Detector performance was quantified to evaluate whether the missed detection rate for upcalls resulted in underestimating full days of NARW presence in our analysis. Three MARUs, manually reviewed for upcall presence in previous studies, were selected and used as ground-truth datasets to compare to our findings for days with NARW acoustic presence. These sites included year-round recordings from the Gulf of Maine (region 1); southern North Carolina (region 10); and Georgia (region 10) 21,48 . Selected days from each recorder were manually screened for NARW upcalls. Every third day of the Gulf of Maine recording was viewed and all upcalls for those days were counted (for detailed results on this analysis, see Bort, et al. 21 ). Output from a detector 70 run on North Carolina and South Carolina units was manually reviewed for daily presence of at least one NARW upcall (for detailed results on this analysis, see Hodge, et al. 48 ). The resulting data were combined and used to generate a ground-truthed dataset of days with NARW presence; this was then compared to the number of days estimated to have NARW presence based on the manually screened output of the LFDCS system.