Spatial and temporal variability in harbour porpoise density: implications for conservation in UK seas

“No good fish goes anywhere without a porpoise.”

-Lewis Carroll

  Abstract

Environmental change and anthropogenic impacts undoubtedly affect species, particularly populations in shelf-seas whose ranges overlap with areas of high recreational and commercial activity. Understanding the degree to which stressors impact a population relies primarily upon knowledge of baselines. Basic ecological data such as long-term changes in distribution and abundance is lacking for most cetaceans, ultimately hindering an in-depth assessment of impacts. Although the harbour porpoise Phocoena phocoena is a widespread predator in the Northeast Atlantic, knowledge gaps remain in their fine-scale, temporal density. This study combined 12 years of sightings data collected by citizen scientists on platforms of opportunity to estimate the relative abundance and density of harbour porpoises along crossings in the Celtic Sea, English Channel, and North Sea. A two-stage process was employed to account for imperfect detection and generate density surface models (DSMs) retaining static and dynamic covariates explaining patterns in distribution. This study shows how harbour porpoise population estimates are sensitive to both space and time. Distributions were influenced by distance to coast, chlorophyll-a, sea surface temperature (SST), and surface-seabed temperature differences. The study highlights persistent high-density areas from the available time series, and identifies a seasonal shift in distribution across the English Channel and North Sea.

Key words: harbour porpoises, density surface models, distance sampling, platforms of opportunity, citizen science

Introduction

Understanding how organisms interact with their environment and respond to change is the foundation of ecology, which provides information vital for management. This necessity poses a challenge for the study of highly mobile species, as their migrations may cross international boundaries and ocean basins, and ranges can exceed thousands of kilometres. Studies are further complicated by the behavior of the target species as well as the logistical and financial costs involved with conducting research at sea over the spatial and temporal scales necessary to evaluate population trends. Available resources are often allocated to research on endangered species. While common species have a less urgent conservation priority, knowledge gaps remain in the understanding of both how vulnerable and how resilient populations are in the face of increasing environmental change (Bearzi et al. 2003). Cetaceans are exposed to a spectrum of anthropogenic stressors including prey exploitation, subsistence and commercial hunting, noise and chemical pollution, ship-strikes, and destructive fishing practices. Information on the distributions and densities of cetaceans – both endangered and common – is necessary in order to assess the scope of these impacts on both regional and global scales, as well as to mitigate threats.

The harbour porpoise Phocoena phocoena is a small cetacean limited to cold temperate and subpolar waters in the Northern Hemisphere. They are the only porpoise species occurring in the Northeast Atlantic, and their expansive distribution and prevalence in coastal waters makes them both important predators and indicator species, as variation in their distribution and density may reflect underlying environmental change (Gilles et al. 2001). Although they occur commonly in shelf-seas, historical hunting pressure on some populations, habitat degradation, and pollution, as well as incidental fisheries mortality in European and UK waters in recent years has given reason for concern (Hammond et al. 2002, Viquerat et al. 2014, Peschko et al. 2016). In the UK alone, it was estimated that nearly 1,500 individuals were bycaught in gillnets during 2015 (Northridge et al. 2016). These continued threats are in spite of government protection, including listings under ASCOBANS (Agreement on the Conservation of Small Cetaceans of the Baltic, Northeast Atlantic, Irish, and North Seas) and the Convention for the Protection of the Marine Environment in the Northeast Atlantic (OSPAR). The European Union’s Directive on the Conservation of Natural Habitats and of Natural Fauna and Flora (Habitats Directive) was formed as the primary regulatory body for conservation and diversity, including the provision of protective actions for cetaceans. Harbour porpoises are listed as species of interest under Annex II and Annex IV, requiring Member States to establish Special Areas of Conservation (SACs) in addition to providing full protection of all life stages throughout their range.

Prioritizing areas of importance requires thorough knowledge on population sizes and trends in abundance; due to their wide range, inconspicuous nature, and small size, however, harbour porpoises are especially difficult to spot at sea (Prescott & Gaskin 1981, Embling et al. 2010, Oakley et al. 2012). The most accurate estimates of absolute abundance come from dedicated, purpose-designed surveys which allow extensive spatial coverage of a study area. In the Northeast Atlantic, the best available information on distributions and population sizes has been provided by SCANS (Small Cetacean Abundance in the North Sea and Adjacent Waters) which surveyed cetaceans in European and UK shelf-seas along randomly placed transect lines (Hammond et al. 2002, Hammond et al. 2013, Hammond et al. 2017). Because purpose-designed surveys are intensive and require employed research platforms (vessels and/or aircraft), they are conducted infrequently and data collection tends to be limited to favourable sighting conditions in summer months. These “snapshots” of abundance lack information on inter- and intra-annual variability and long-term trends; when analysed in isolation, they lack detail on population dynamics and stability (Evans & Hammond 2004).

Patterns in seasonal and fine-scale distributions and abundances of harbour porpoises are poorly understood (Bailey & Thomson 2009) and only a few candidate SACs have been designated in the UK (IAMMWG 2016). Where large-scale surveys of absolute abundance leave inter-annual gaps, long-term studies in more specific areas can be employed to assess variation. Vessels such as ferries and cruises, or platforms of opportunity, can serve as an inexpensive base for conducting research along consistent crossings over extensive time scales. In this way, data collected aboard such vessels increases temporal coverage, and this regular monitoring may provide more power for early detection of population change – in addition to its possible connection to anthropogenic threats – in the stretches between large-scale surveys (Viqueret et al. 2014). Fixed routes characteristic of platforms of opportunity are considered independent of the occurrence of the study species, and therefore allow abundance to be extrapolated along the surveyed transects from sightings and effort data (Isojunno et al. 2012, Williams et al. 2006).

By combining data on spatiotemporal distributions with multivariate regression modelling, abundance can be estimated to look at density patterns. Oceanographic and environmental variables affect productivity, drive the availability of prey species, and influence the distribution of predators (Bjørge 2001). As small, energetically demanding species, harbour porpoises must spend an increased amount of their time exploiting patchily distributed prey (Prescott & Gaskin 1981, Cox et al. 2017).  Environmental variables that influence primary productivity or serve as proxies for such processes have been shown to influence patterns of cetacean occurrence by shifting distributions of prey species (Kiszka et al. 2007, Embling et al. 2010, Isojunno et al. 2012, Booth et al. 2013, Gilles et al. 2014, McClellan et al. 2014). Density surface models (DSMs) provide a robust method for explaining patterns of relative abundance from non-designed surveys in complex, dynamic environments. The provided output can be used to predict areas of importance and identify critical habitats for mobile species by relating static and dynamic covariates to density estimates (Dransfield et al. 2014, McClellan et al. 2014, Harvey et al. 2017).

This study aims to derive valuable information related to harbour porpoise distribution and abundance in the Celtic Sea, English Channel, and North Sea from data collected by citizen scientists on platforms of opportunity by: (1) accounting for imperfect detection during line-transect surveys, and (2) generating spatially-explicit models of density. The data provides more detailed information on the fine-scale, long-term harbour porpoise population trends in UK and European waters, and in combination with findings from purpose-designed surveys, can be used to manage important areas and mitigate threats.

Methods

Study area & data collection

The oceanography of the Northeast Atlantic is dominated by the Gulf Stream and North Atlantic Drift. This brings relatively warm water to shallow shelf-seas, flowing north around the British Isles and into the South-Western Approaches of the Celtic Sea and English Channel (Hardisty 1990). ORCA’s trained citizen scientists conducted line-transect surveys aboard ferries in three regions surrounding the UK: in the Celtic Sea, across the English Channel, and across the North Sea (Figure 1). The Celtic Sea is a shallow embayment of the continental shelf between Ireland, Wales, and England in OSPAR region III. Survey effort in this sea was exclusive to a crossing from Penzance to St. Mary’s on the Isles of Scilly. Ferries departing from the south coast crossed the English Channel along five different routes to France and Spain. OSPAR region II, the North Sea, is exposed to land along most of its boundaries and influenced by circulation from the Strait of Dover in the south and the Atlantic Ocean and Norwegian Sea in the north. Surveys were conducted along five routes in this region from mainland UK to the Shetland Islands, Norway, Sweden, Denmark, and the Netherlands.

Figure 1: Study sites in the Celtic Sea, English Channel, and North Sea. Black lines show effort along ferry routes from 2006-2017 and grey lines outline the boundaries of OSPAR regions (base map made using GEBCO bathymetry data).

Data related to cetacean sightings and environmental conditions were collected between 2006 and 2017 using distance sampling methodology (Buckland et al. 2001). Two observers were stationed on either side of the platform at the bow of the vessel, monitoring on either side of the transect from 270° to 90°. Survey effort was tracked by GPS along with vessel-specific variables such as speed and platform height. Environmental conditions including visibility, Beaufort sea state, relative swell height, precipitation, and glare were noted throughout the survey. When a group was first detected, the sighting angle from the transect was measured using an angle board and reticle distance was estimated by binoculars. Species, behaviour (such as fast swimming, feeding, or bow-riding), group size, and cue (breaching, dorsal fin, blow) were also recorded.

Detection function estimation

The detection function fits the relationship between both the probability of detection and the perpendicular distance from the transect line. Reticle measurements and angles from the transect to each harbour porpoise sighting were used to calculate perpendicular distances to detected groups. Sighting cue and behaviour of the individuals were considered for relevant sources of bias, such as responsive movement. After preliminary data exploration, the data was stratified to fit models specific to the Celtic Sea, English Channel and North Sea routes as delineated by their respective OSPAR regions.

The package ‘Distance’ (Miller 2017) for R (R Core Team 2017) was used to fit a detection function unique for each of the three study regions by estimating the probability of detecting harbour porpoises from a given distance. The calculation assumes that all groups directly ahead of the vessel on the line of travel are detected with certainty (g(0)=1) (Buckland et al. 2001). Although such certainty is highly unlikely, especially for inconspicuous harbour porpoises, the true value of g(0) could not be calculated with the employed single platform survey technique (Berrow et al. 2014). The second assumption of distance sampling relies on recording distances based on an animal’s original location (Buckland et al. 2001). Harbour porpoises are not only highly mobile predators, but are also sensitive to noise (Bernd & Evans 2002); because of this, the possibility that individuals were detected after responsive movements, such as swimming away from or towards the vessel, could not be eliminated. To reduce bias in this data, exact measurements were smeared into distance intervals, or ‘bins’, by selecting cutpoint boundaries between 0 m and the truncation distance (w). Truncating outliers improves the fit of detection functions by removing sightings recorded at the greatest distances. Truncation distances and cutpoint locations were changed by improving goodness of fit. The resulting functions removed the farthest 5-10% of sightings, resulting in w=650 m in the Celtic Sea, w=405 m in the English Channel, and w=600 m in the North Sea.

Detection is influenced to some degree by complex conditions related to both survey methods as well as localized oceanographic and environmental conditions. Research has shown that there is a gradual reduction in detectability with increasing distance (Hammond et al. 2002) and a sharp decrease for surveys conducted in Beaufort sea states ≥ 3 (Teilmann 2003). Several model variations were run to test the significance of vessel speed, platform height, visibility, swell, Beaufort sea state, and group size on sightings. Using both hazard-rate and half-normal key terms, models with different bin sizes were built as functions of perpendicular distances and covariate values. A model was selected for each region by retaining those covariates which best explained detection (chi-square for binned data, p<0.05) and by minimizing Akaike’s Information Criterion (AIC). The effective strip width (ESW), which is the distance from the transect where as many porpoises were detected beyond it as were missed within, was used to estimate detection probability by

g(x) = probability of detection directly on the transect, and

ESW = effective strip width.

The ratio of recorded harbour porpoise sightings and detection probability provided an estimate of relative abundance (Buckland et al. 2001):

ddf.gof(ph.pod.speed_celtic$ddf)

summary(ph.pod.speed_celtic)

Satellite data: environmental and oceanographic data were downloaded from external sources as .nc files. The following code was used to join satellite data to relevant segments by x and y coordinates, and write the output into a .csv file for export into ArcGIS.

###import multiple netCDF files as rasters (Copernicus data requires two stacks: one each for the reanalysis and forecast products)

celtic1 <- read.csv(“celtic_reanalysis.csv”)

coords1 <- celtic1[c(“x”, “y”)]

celtic2 <- read.csv(“celtic_forecast.csv”)

coords2 <- celtic2[c(“x”, “y”)]

###merge multiple netCDF layers and save as one raster

#load packages

library(ncdf4)

library(raster)

#create a raster stack of temperature at the seabed for each month from each oceanographic model

#reanalysis (the majority of months are removed from code to save space)

x1<- stack(

raster(“BT_200603.nc”),

…

raster(“BT_201310.nc”))

##forecast model

x2 <- stack(

raster(“ST_201503.nc”),

…

raster(“ST_201710.nc”))

# Name each layer

##reanalysis model

names(x1) <- c(‘2006_Mar’, …,’2013_Oct’)

##forecast model

names(x2) <- c(‘2015_Mar’ ,…,’2017_Oct’)

#write rasters

writeRaster(x, filename = ‘bottomTemp_CELTIC_reanalysis.nc’, format =”CDF”, overwrite=TRUE)

writeRaster(x, filename = ‘bottomTemp_CELTIC_forecast.nc’, format =”CDF”, overwrite=TRUE)

#write .csv files to attach data to Celtic segments based on specified coords

write.csv(celtic1, “bottomTemp_reanalysis_celtic.csv”)

write.csv(celtic2, “bottomTemp_forecast_celtic.csv”)

Density surface modelling: the second stage of the DSM process involves fitting GAMs to segment data. These models include corrected abundances provided by the detection function and significant environmental covariates that influence the distribution of harbour porpoises. The three models were fit in the package ‘dsm’ (Miller 2013). The DSMs specific to the Celtic Sea, English Channel, and North Sea all included smooths of significant covariate using thin-plate regression splines. A subset of the code used in model selection is attached for the Celtic Sea.

#load packages

library(dsm)

library(rgdal)

library(knitr)

#attach segment data, observation data, and prediction data separately

segs <- read.csv(“segdata_celtic.csv”)

obvs <- read.csv(“obvsdata_celtic.csv”)

preddata <- read.csv(“preddata_celtic.csv”)

#check for collinearity (make separate spreadsheet for covariates)

covar <- read.csv(“covar_celtic.csv”)

covar <- read.csv(“covar_celtic.csv”)

pairs(covar, panel=panel.smooth)

cor(na.omit(covar))

###fit density surface models

#backward stepwise selection…build a model with all covariates, check significance, remove, refit, repeat

dsm.tp.tw <- dsm(abundance.est~s(coast, bs=”tp”)+s(log(SST), bs=”tp”)+s(slope, bs=”tp”), ddf.obj=df_celtic, segment.data = segs,

observation.data = obs, family=tw(), method=”REML”, gamma=1)

summary(dsm.tp.tw)

#check for concurvity in ‘good’ models

vis.concurvity(dsm.tp.tw)

plot(dsm.tp.tw, shade=TRUE, shade.col=”green”)

gam.check(dsm.tp.tw)

#run model variations, adjusting shrinkage, smoothing parameters, knots, maximum penalties, family, etc. and assess fit

#select model by REML, AIC, %deviance explained – maintain parsimony

dsm5.2 <- dsm(abundance.est~s(coast, bs=”tp”)+s(log(SST), bs=”tp”)+slope, ddf.obj=df_celtic, segment.data = segs,

observation.data = obs, family=tw(), method=”REML”, gamma=1)

summary(dsm5.2)

plot(dsm5.2, shade=TRUE, shade.col=”turquoise2″)

gam.check(dsm5.2)

Density predictions: the final models were used to predict density on each segment by using an offset of segment-specific area. The following code was used to estimate density. A subset is attached for the Celtic Sea.

preddata <- read.csv(“preddata_celtic.csv”)

#prepare preddata, using segment specific area (2wL) as the offset

celtic.pred <- predict(dsm5.2, preddata, preddata$area)

#load packages

library(viridis)

library(raster)

celtic.pred <- as.data.frame(segs2, xy=TRUE)

###offset by area of each prediction segment

celtic.pred$off.set <- preddata$area

#predict, find abundance per segment

pp <- predict(dsm5.2, celtic.pred)

##Sum these for total abundance in Celtic Sea

sum(pp, na.rm=TRUE)

#assign the predictions to the prediction grid data.frame

celtic.pred$abundance.est <- pp

write.csv(celtic.pred, “segment_density_celtic.csv”)

#look at abundance specific to year and/or months

sum(celtic.pred$abundance.est[celtic.pred$year == “2009”])

sum(celtic.pred$abundance.est[celtic.pred$year == “2010” & celtic.pred$month == “8”])

###calculate variance

var.dsm5.2 <- dsm.var.gam(dsm5.2, celtic.pred, off.set=preddata$area)

summary(var.dsm5.2)

Supplementary figures

Table S1: Effort for each survey month and year between 2009 and 2017 in the Celtic Sea.

Table S2: Effort for each survey month and year between 2006 and 2017 in the English Channel.

Table S3: Effort for each survey month and year between 2006 and 2017 in the North Sea.

Table S4: Total abundance for each survey month and year between 2009 and 2017 in the Celtic Sea.

Table S5: Total abundance for each survey month and year between 2006 and 2017 in the English Channel.

Table S6: Total abundance for each survey month and year between 2006 and 2017 in the English Channel.

Figure S1: Pairs plot for covariates in the Celtic Sea. A test of collinearity showed that no two variables were highly correlated. All Pearson correlation coefficient were < 0.70.

Figure S2: Pairs plot for covariates in the English Channel. A test of collinearity showed that no two variables were highly correlated. All Pearson correlation coefficient were < 0.70.

Figure S3: Pairs plot for covariates in the North Sea. A test of collinearity showed that no two variables were highly correlated. All Pearson correlation coefficient were < 0.70.

Table S7: Specifications and output for the final DSM in each region.

Figure S4: The ‘best’ model in the Celtic Sea retained slope as a covariate, but goodness-of-fit was improved by including it as linear term. REML scores could not be compared between the models with and without smooths, so a slight decrease in AIC was used for selecting between the two.

Figure S5: Visual results of a test of concurvity on smooths in the final Celtic Sea model. The colour scale represents the degree to which the apparent relationship in each smooth is explained by a different smooth in the same model. A value of 1 suggests a major problem with concurvity.

Figure S6: Visual results of a test of concurvity on smooths in the final English Channel model. The colour scale represents the degree to which the apparent relationship in each smooth is explained by a different smooth in the same model. A value of 1 suggests a major problem with concurvity.

Figure S7: Visual results of a test of concurvity on smooths in the final North Sea model. The colour scale represents the degree to which the apparent relationship in each smooth is explained by a different smooth in the same model. A value of 1 suggests a major problem with concurvity.