Evidence for seasonal cycles in deep-sea fish abundances: A great migration in the deep SE Atlantic?

1. Animal migrations are of global ecological significance, providing mechanisms for


| INTRODUC TI ON
The deep oceans (>200 m depth) are temporally dynamic and spatially heterogeneous environments that together comprise the largest living space on Earth (Smith, Leo, Bernardino, Sweetman, & Arbizu, 2008). The deep sea provides a variety of goods and services to humankind, ranging from alternative mineral or food resources to carbon sequestration, all of which are becoming increasingly important as human activities in the oceans expand to greater depths (Mengerink et al., 2014). Nonetheless, most deep-sea ecosystems remain relatively poorly understood, particularly at the temporal and spatial scales of relevance to human impacts (Ruhl et al., 2011;Webb, Vanden Berghe, & O'Dor, 2010).
The oil and gas industry is moving operations into increasingly deep-ocean regions, for which there are very few baseline datasets available (Cordes et al., 2016;Peterson et al., 2012). Such a lack of ecological data can present a huge challenge to predicting and managing the extent and severity of potential impacts even in relatively well-studied regions. This was highlighted following the 2010 Deepwater Horizon oil spill for example, which occurred at c. 1,500 m water depth in the Gulf of Mexico and was an unprecedented deep-ocean event (Peterson et al., 2012), with long-term consequences that remain poorly understood. In other areas of major deep-sea oil development such as the Brazilian and West African continental slopes (Cordes et al., 2016) even fewer baseline data are available, and the need for long-term monitoring studies conducted at appropriate temporal and spatial scales is clear. Owing largely to logistical and technical challenges, short-term, localized studies have been the norm in deep-sea research and long-term biological surveys remain rare, though paleoecological studies can provide valuable retrospectives on historic conditions and temporal change at decadal and greater scales (e.g. Yasuhara, Cronin, deMenocal, Okahashi, & Linsley, 2008;Yasuhara, Okahashi, Cronin, Rasmussen, & Hunt, 2014). While there has been increasing interest and investment in long-term studies in recent years, biological observations remain under-represented, and there are notably fewer long-term monitoring stations in the southern hemisphere, compared to the northern hemisphere (Levin et al., 2019). Those long-term studies that do exist are proving highly valuable in understanding the extent of coupling between surface-ocean processes and the deep oceans (Glover et al., 2010;Ruhl et al., 2011).
The majority of demersal deep-sea fauna are entirely reliant on detrital flux from the surface waters for energy. Most of this energy enters the deep sea as a diffuse flux of sinking phytodetritus, comprised largely of dead plankton, faecal matter and other biological waste products . An extensive literature has demonstrated clear associations between temporal and spatial changes in surface-ocean production and the quantity of detrital flux reaching the deep-sea floor (Lampitt et al., 2010;Smith et al., 2006). Similarly, associations between the biomass, abundance and biodiversity of deep-sea benthic invertebrate fauna and overlying productivity have been demonstrated throughout the deep oceans (Billett, Bett, Reid, Boorman, & Priede, 2010;Johnson et al., 2007;McClain, Allen, Tittensor, & Rex, 2012;Smith, Ruhl, Kahru, Huffard, & Sherman, 2013). In the Arctic, data from the Hausgarten observatory (which has been operational since 1999; Soltwedel et al., 2005) suggest a similar link between the quantity of detritus reaching the seafloor (which is reliant on the extent of sea ice cover in this region) and the benthic community composition at bathyal depths (Bergmann, Soltwedel, & Klages, 2011). However, the majority of studies conducted to date are focused either on relatively short time frames, or on inter-annual patterns, rather than intra-annual (seasonal) patterns. One recent study made use of 2.5 years of data from Berkeley Canyon (985 m water depth; NE Pacific Ocean) and reported significant temporal changes in the epibenthic, megafaunal invertebrate assemblage at 'broad' time-scales (months-years) which was correlated weakly with both local and epipelagic environmental variables (Chauvet, Metaxas, Hay, & Matabos, 2018).
Fewer data are available for demersal fishes, although abundances have been correlated with spatial gradients in surface productivity (Armstrong, Bagley, & Priede, 1992;Henriques, Priede, & Bagley, 2002;Thurston, Bett, & Rice, 1995), and with inter-annual variation in invertebrate and carrion availability (Drazen, Bailey, Ruhl, & Smith, 2012). Since demersal fishes within deep-sea ecosystems are typically predators or scavengers of large carrion (Drazen & Sutton, 2017), such correlations are suggestive of strong bottom-up trophic control within deep-sea food webs. In each case, the data highlight the potential for human activities to significantly alter the structure and functioning of trophic control could lead to changes in the abundance of fishes across the seafloor by affecting secondary production of prey species and/or carrion availability for example. 5. In summary, we present the first evidence for seasonally recurring patterns in deep-sea demersal fish abundances over a 7-year period, and demonstrate a previously unobserved level of dynamism in the deep sea, potentially mirroring the great migrations so well characterized in terrestrial systems. deep-sea ecosystems either directly or through the indirect effects propagating from the surface ocean (Bergmann et al., 2011;Ruhl, Ellena, & Smith, 2008;Smith et al., 2008Smith et al., , 2013. Here, we quantify temporal patterns in demersal fish abundance at two sites on the Angolan continental slope (c. 1,400 m water depth) using 7.5 years of time-lapse photographic data from the Deep-ocean Environmental Long-term Observatory Systems (DELOS; Vardaro et al., 2013). These patterns are compared to corresponding temporal variations in the flux of particulate organic carbon (POC) to the seafloor using methods described in Behrenfeld and Falkowski (1997) and Pace, Knauer, Karl, and Martin (1987).

| Study area and data sources
Time-lapse photographs were collected from the 'Near Field' (NF) and 'Far Field' (FF) DELOS observatories (Vardaro et al., 2013)  Angolan continental shelf, with NF located c. 50 m from an active oil well (7.90°S, 12.14°E), and FF located at a reference site c. 14 km to the south-east (7.95°S, 12.28°E; Figure 1). Each camera was programmed to photograph at 3-hr time intervals between February 2009 and August 2011 and at 9-hr intervals thereafter to better conserve battery life. Intermittent camera or flash failures resulted in a non-continuous time series of images over the study period (Table S1). In cases where the camera flash failed to fully illuminate the field of view (i.e. where the flash and camera were not synchronized correctly), images were processed manually to increase brightness. Given the comparatively large size of the observed fishes relative to the cameras' field of view, we do not believe these technical problems would have biased the findings presented here.
Net primary production (NPP) values were extracted from the vertically generalized production model (resolution: 1/6 × 1/6º), described by Behrenfeld and Falkowski (1997), based on MODIS R2014 data as downloaded from http://www.scien ce.orego nstate.edu/ocean.produ ctivity. NPP estimates were taken from the cell located closest to the NF site. A time series of NPP data was collated for July 2002 to July 2016. Estimated POC flux to the seafloor was calculated for a water depth of 1,400 m, using the algorithm in Pace et al. (1987).

| Time-series analysis: Fish counts
The number of fishes per photograph was counted from each observatory and binned by day to produce two irregular time-series spanning approximately 7.5 years (Table S1). Since the NF data provided the longest time series for analysis and contained the largest number of fishes (Table S1), they were used to identify any inter-and intra-annual temporal patterns in the data. The NF data were divided into training (c. 60%) and test (c. 40%) datasets (Table S1), and generalized additive models (GAMs; Wood, 2006) fitted to the total fish counts observed each day, modelled against Julian day (to test for intra-annual patterns) and the elapsed number of days from 1 January 2009 (to test for long-term temporal patterns). The test data were selected to cover periods close to the beginning, middle and end of the 7-year time series to ensure that they were as representative of the study period as possible.
A GAM was then fitted to the NF training data to test for the presence of both inter-and intra-annual patterns (Equation 1) using the gamlss (Rigby & Stasinopoulos, 2005)  fitted as the number of days since 1 January 2009; the intra-annual smoother was fitted as the number of days since 1 January within each year. The intra-annual smoother was fitted using cyclic P-splines and the best fit selected automatically during the initial fitting process, while the inter-annual smoother was initially fitted using cubic splines and constrained to three degrees of freedom to avoid over-fitting. In each case, ≤50 iterations of the Rigby and Stasinopoulos algorithm (RS method; Rigby & Stasinopoulos, 2005) and ≤200 iterations of the Cole and Green algorithm (CG method; Rigby & Stasinopoulos, 2005) were used to fit the models.
Term selection was conducted by comparing the full model to two models in which either the inter-or intra-annual smoother had been removed, and an intercept-only model. The global deviances (GD) for each model as fitted to the test data are presented in Table S2. A sensitivity analysis was then conducted on the bestfitting model by manually altering the degrees of freedom (knots) of the fitted smoother(s) and comparing the goodness-of-fit to the test data as before ( Figure S1).
In all cases, model selection was conducted by comparing the GD and mean squared error of prediction scores of each training model, based on how well they each fitted the test data. Model validation was conducted by visual examination of the normalized quantile residuals and worm plots (van Buuren & Fredriks, 2001) plotted against the fitted values and fitted variables. Semi-variograms were produced using the gstat (Pebesma, 2004) package in r software to confirm that no temporal autocorrelation was evident in the residuals ( Figure S2).
To test whether a similar intra-annual pattern existed in the FF data, a time series of predicted data were generated from the NF model and correlated against the FF count data using the gamlss function in R software, to ensure consistency with the NF results.
This model was also fitted using a negative binomial distribution with log-link function.

| Time-series analysis: POC flux to the seafloor
A GAM was fitted to the log-transformed mean monthly POC flux to the seafloor estimates to test for the presence of both interand intra-annual patterns (as in Equation 1, where Y = ln POC flux) using the gamlss package in r software. This was done to identify any overall trend in the data and extract the seasonal component of the time series. The inter-annual smoother was fitted using cubic splines as the number of months since January 2009; the intraannual smoother was fitted using cyclic P-splines as the number of months since January within each year. Models were fitted using a Gaussian distribution with identity link function. Term selection was conducted using AICc scores, and model validation was conducted as described previously.
Comparison of the POC flux data to the fish densities was conducted in two stages, by first correlating the predicted fits and then the raw data. The strength of the correlation between the temporally lagged, predicted fits was calculated using Pearson's correlation coefficient and the 'cor' function in r software. To compare the raw data, the fish abundances from the NF observatory were converted to densities per calendar month (i.e. summed fish counts divided by the number of images taken per month). Since the raw data could not be cross-correlated directly as a result of a number of data gaps in the NF time series, a series of time lags between 0 and 12 months were individually calculated using the 'cor.test' function in r software to estimate Pearson's product moment correlation (Table S3).

| Multivariate analyses
To determine whether any spatial or temporal changes occurred in the fish assemblage composition, fish abundances were summed per calendar month for each family. Both square root and presenceabsence (PA) transformations were applied to the count data to weight the analyses towards the more common and the rare taxa, respectively. While it was possible to identify some taxa to genus or species, the fish fauna from this region are currently poorly described, meaning that other taxa (particularly the Macrouridae) could not be reliably identified beyond family from photographs. All analyses were therefore conducted at the Family level of taxonomic resolution. The monthly time period was selected as a compromise between temporal resolution and appropriate sampling unit size.
Months with no fish observations were excluded from these analyses. The data were converted to a similarity matrix using Bray-Curtis similarity index (for square root transformed data) or Sørensen's index (for PA data) prior to analysis using primer 7 software with PERMANOVA (Clarke & Gorley, 2015).
Three temporal models were generated using dummy data. The first modelled a linear trend, and was coded as increasing integer values increasing from 1 (January 2009) to 91 (July 2016). The second modelled an annual cycle, in which June to December were coded as integers from 1 to 7; January to May were coded from 6 to 2. The third modelled a biannual cycle, in which September to December were coded from 1 to 4; January to March from 3 to 1; April to June from 2 to 4; and July and August as 3 and 2, respectively. The dummy values were selected to correspond to the peaks in fish abundance that were detected by the time-series analyses.
Differences between the observatories were modelled using binary codes (i.e. NF = 0; FF = 1). A fourth model was also fitted using estimated POC flux to the seafloor as the explanatory variable at a 4-month temporal lag to match observations from the time-series analyses. Species contributing most to significant temporal changes were identified using BEST (BVSTEP) analysis in primer 7 software.
The statistical significance of the temporal patterns was tested using distance-based linear models (Legendre & Anderson, 1999).
(1) Y ∼ s 1 + s 2 , Models were fitted using the BEST selection procedure and AICc as the selection criterion.

| Time-series analysis: Fish counts
Between February 2009 and July 2016, 502 fishes belonging to at least 10 families were observed at the NF, and 217 at the FF from a total of 6,783 and 5,920 photographs, respectively ( Figure 2; Table 1).
The temporal GAM that best fit the test data from the NF indicated a recurrent, bimodal intra-annual pattern (Figure 3a; Table S2), with the highest fish abundances occurring in late November (Julian day 327) and a smaller, secondary peak occurring in June (Julian day 163).
Inter-annual variation in counts was not well predicted by the training data, and the residual variability remained high (Figure 3a). Despite the lower numbers of fishes observed at the FF, a positive correlation was identified between the FF fish counts and the predicted values from the NF GAM (t-value = 3.04, p = 0.002).

F I G U R E 2
Example images of the fish morphotypes observed from each family during the present study.

| Time-series analysis: POC flux to the seafloor
The best-fitting model of POC flux to the seafloor indicated a statistically significant bimodal seasonal pattern, with a primary peak occurring in August and a secondary peak occurring in February/ March (Figure 3b). Significant positive correlations between monthly averaged fish density at NF and estimated POC flux to the seafloor occurred at lags of 4 and 5 months, with the strongest relationship at a 4-month lag (r = 0.383, p < 0.001; Figure 3c; Table S3).

| Multivariate analysis
The best-fitting multivariate models describing the composition of fish fauna at the NF and FF based on both PA and square root

| D ISCUSS I ON
Understanding seasonal changes in the abundances of deep-sea fishes provides important insight into how deep-water ecosystems are structured, their inherent variability and a better understanding of the role fishes may play in the transfer and redistribution of energy within the deep-ocean realm. Here, we identified a total of 502 demersal fishes from at least 10 families over a 7.5-year period on the Angolan margin and identified a recurring intra-annual (seasonal), bimodal pattern in the total density of fish at both DELOS observatories. Inter-annual changes (between years) were not predicted well, but this was likely a consequence of the relatively short time series available for the present study. Some inherent stochasticity is also inevitable when modelling a highly mobile, low-density animal assemblage, and additional variability will also come from other, unmeasured drivers operating at different spatial, temporal or taxonomic scales than those considered by any single study. For the purposes of this study, we focus only on the observed intra-annual, seasonal patterns.
The observed and predicted fish densities were significantly correlated with satellite-derived estimates of primary production overlying the study area, with the strongest correlation occurring and discharge from the Congo River by approximately 2 months (Biemans et al., 2009). Both these observations lend weight to the hypothesis that the underlying driver of fish densities is a general one.
An indirect, mechanistic relationship between POC flux to the seafloor and demersal fish abundance could occur if temporal variations in the POC flux are propagated through the food web as a bottom-up process, following time lags that allow for POC sinking rates (c. days-weeks; e.g. Berelson, 2002) and secondary production by prey (e.g. Rowe, 2013). This hypothesis is supported by both empirical (Johnson et al., 2007;Smith et al., 2008) and theoretical studies (Marczak, Thompson, & Richardson, 2007;Polis, Anderson, & Holt, 1997) that indicate strong bottom-up trophic control within food webs that rely on detritus, such as those of the deep-sea floor.
While these precise rates are not presently known for this region, a 4-month lag between POC production and the arrival of mobile fish predators seems a plausible time frame.
The fact that no intra-annual variation in species composition was detected over the study period lends further weight to the hypothesis that the fish assemblage was likely responding to a general driver capable of simultaneously influencing multiple species.
While the present analyses were restricted to the family level, the observed fish taxa belonged to a range of trophic guilds including scavengers, and benthic and bentho-pelagic predators (Drazen & Sutton, 2017). The lower density of fishes observed at the FF and the difference in faunal composition between observatories indicate that the fish assemblage is influenced by different processes at these spatial scales, with proximity to an active oil field being one obvious potential candidate. The absence of an increasing temporal trend in the observed density of fishes at either the NF or FF suggests that neither the oil production infrastructure nor the DELOS observatories promote strong reefing behaviours in the observed deep-water fishes. Similar conclusions have been previously drawn with regard to fishes inhabiting the abyssal NE Pacific (4,100 m water depth; Vardaro, Parmley, & Smith, 2007).
Annual migrations are well documented among shallow-living marine fishes (Beamish, McFarlane, & King, 2005), and although similar migrations are expected to occur among deep-sea fishes, there is limited empirical evidence to support this idea (Campana, Valentin, Sevigny, & Power, 2007). Where previous studies have typically inferred seasonal patterns within individual species indirectly from otolith microchemistry (Campana et al., 2007), from fisheries landings data (Kuo & Tanaka, 1984) or from activity rates and body sizes observed at baited cameras (Jones, Collins, Bagley, Addison, & Priede, 1998;Priede, Deary, Bailey, & Smith, 2003;Smith, Priede, Bagley, & Addison, 1997), the present study provides evidence of repeating seasonal patterns of abundance and suggests recurrent seasonal movements of the deep-sea fishes as a potential explanation.
While the life-history traits of deep-sea fish vary (Drazen & Haedrich, 2012), the longevity of most demersal species means that processes such as recruitment or mortality can be excluded as causal drivers, and the observed intra-annual patterns must reflect fish movements into and out of the observed area. Any causal mechanism is unlikely to be directly related to oil production activities, since such activities would not be expected to extend to the FF. Similarly, it is difficult to imagine sources of sampling biases that would have given rise to the bimodal patterns observed here. The operation of a work-class ROV during the annual servicing process is potentially the most likely and impactful source of bias, but the servicing periods do not correspond consistently with peaks or troughs in the time series (Figure 3; Table S1), and would not account for the bimodal pattern observed even if they do cause a temporary disturbance. Variations in swimming speeds are also unlikely to account for the observed patterns, since the increased chance of passing through the field of view at greater speeds would likely be offset by a lower time within the observed field of view, and therefore the same overall chance of being photographed. Consequently, the most likely explanation for the observed pattern is that it reflects real differences in the density of fishes within the study areas.
If the observed variations in fish abundances were driven by relocations, rather than recruitment or mortality, then three nonexclusive hypotheses might be considered: (a) movement up or down the continental slope; (b) movement across the slope and (c) movement into the water column. POC flux to the seafloor is influenced by a complex range of physical (e.g. temperature, oxygen concentration or water depth), chemical (e.g. POC composition) and biological process (e.g. microbial remineralization) that influence the particle size, biomass and 'quality' of POC that is exported from the surface ocean and ultimately becomes available to the benthos (e.g. Belcher Ocean several hundreds of metres above and below the depth of the DELOS observatories (Jamieson, Linley, & Craig, 2017;Priede et al., 2010).
Along-slope movements might be cued by recurring spatial variations in surface ocean primary production, such as zonal variations in phenology or localized enhanced productivity associated with upper-ocean boundary fronts or eddies and subsequent POC flux to the seafloor (e.g. Cram et al., 2018;Gaube, McGillicuddy, Chelton, Behrenfeld, & Strutton, 2014), or with variations in Congo River outflow (Biemans et al., 2009). While the deep seas are characterized by relatively stable abiotic conditions at very broad spatial and long temporal scales (e.g. Yashuhara & Danovaro, 2016), such mesoscale variability may be more relevant localized, short-term drivers of behavioural responses in mobile fauna. The option for a general synchronous shift from a benthic habit to a pelagic habit (i.e. movement in to the water column) seems less likely given the broad spectrum of species involved. While these three hypotheses are highly simplified, they nonetheless provide a basis for further mechanistic studies in the region.
An alternative mechanism to direct or indirect responses to POC flux to the seafloor and subsequent benthic production may be the response to temporal and/or spatial variations in the supply of carrion to the deep-sea floor. A long-term study in the California Current region of the NE Pacific suggested that the abundance of deep-sea demersal macrourids (Coryphaenoides spp.) covaried with the abundance of migratory Pacific hake in the epipelagic zone, which provide an important, transient carrion resource at depth . A similar situation may exist in Angolan waters, which support large stocks of sardinella, sardine and horse mackerel and whose distributions vary in response to both seasonal changes in plankton availability, temperature and upwelling events (Angelini & Vaz-Velho, 2011;Binet, Gobert, & Maloueki, 2001), as well as to longer-term variations driven by the Atlantic El Niño that influences precipitation rates and outflow from the Congo River (Binet et al., 2001). Similarly, large elasmobranch carcasses (i.e. planktivorous mobulid rays and a whale shark) have been observed on the deepsea floor in Angolan waters (Higgs, Gates, & Jones, 2014), which the authors attributed to close linkages between elasmobranch abundance and seasonal variations in sea surface temperature and prey availability. Humpback whale migrations occur offshore along the Angolan margin in spring and winter each year (Weir, 2011)  Further trophic studies would be useful in addressing the relative importance of benthic production, bentho-pelagic production, and carrion supply to the demersal fishes in this region.
While seasonal changes in primary productivity appear the most likely candidate driver, it is possible that the fish fauna may be responding to other environmental cues that vary over similar time-scales, or to multiple drivers that occur out-of-phase with each other. One candidate, albeit unlikely, alternative driver may be the West African oxygen minimum zone (OMZ), which occurs at approximately 350 m depth (25/75% quartile boundaries: c. 293-416 m) and varies seasonally in location and extent within the area of the DELOS observatories (Chapman & Shannon, 1987;Helly & Levin, 2004). Although probably unlikely given the relatively thin OMZ in the West African region, and the relatively limited extent of hypoxia (Helly & Levin, 2004), oxygen is nonetheless vital to the survival of all animals and widespread hypoxia meets the criteria of a general driver that would force organisms to move away from the affected area if they could not tolerate the conditions. It is not clear that fishes would be unable to tolerate the conditions present in the study region (Gallo & Levin, 2016). Additionally, to be observed at the depths of the DELOS observatories, the effects of the OMZ would need to be propagated approximately 900 m downslope, and potentially some distance across the slope (depending on the distance from the DELOS locations to the edge of the OMZ). The mechanism for such an effect is unknown, but could feasibly involve a cascading densitydependent faunal response if mobile organisms were forced into deeper waters to avoid increased rates of competition or predation, for example, caused by an expanding OMZ. As with changes in primary productivity, such an effect would not need to impact the demersal fauna directly and could also operate via the redistribution of prey resources for example. Alternatively, the presence of an OMZ could affect remineralization of POC through the water column and influence the quality of resources reaching the benthos (e.g. Cram et al., 2018), and potentially influence the fish fauna through changes in secondary production as discussed previously.
Regardless of the underlying mechanism, the consequences of widespread movements of fish in the deep sea could be highly important in the redistribution of carbon and energy across the pelagic-benthic boundary (Trueman, Johnston, O'Hea, & MacKenzie, 2014), and across the seafloor (e.g. by excretion while travelling away from a foraging area). Our observations of a recurrent, intra-annual change in the densities of demersal fishes on the Angola Margin reveal a level of dynamism that has not previously been observed at these scales across the deep-sea floor and raise new questions regarding the functional role of demersal fishes in structuring deepsea ecosystems. For example, large, mobile, generalist consumers and migratory animals have been shown to convey stability in terrestrial and freshwater trophic networks by their ability to adapt to changing resource availability and mute local responses to resource pulses over a range of spatial and temporal scales (Bauer & Hoye, 2014;McCann & Rooney, 2009). Better understanding the temporal dynamics and typical foraging ranges of large, long-lived demersal fishes could therefore be extremely important in understanding patterns and drivers of connectivity between deep-sea assemblages at ocean-basin scales, and the rates and spatial extents over which surface-derived energy can be dispersed and cycled through deepocean ecosystems.

ACK N OWLED G EM ENTS
The authors would like to thank Mr John Polanski and Dr Stewart Chalmers for their assistance with data collection and processing, as well as the captains and crews of the field support vessels for their expertise in deploying and recovering the DELOS observatory modules. This research was funded by BP Exploration and BP Angola, and

DATA AVA I L A B I L I T Y S TAT E M E N T
The fish density data used in the analyses of this paper are available