A Novel, Unbiased Analysis Approach for Investigating Population Dynamics: A Case Study on Calanus finmarchicus and Its Decline in the North Sea.

Marine populations are controlled by a series of drivers, pertaining to both the physical environment and the biological environment (trophic predator-prey interactions). There is heated debate over drivers, especially when trying to understand the causes of major ecosystem events termed regime shifts. In this work, we have researched and developed a novel methodology based on Genetic Programming (GP) for distinguishing which drivers can influence species abundance. This methodology benefits of having no a priori assumptions either on the ecological parameters used or on the underlying mathematical relationships among them. We have validated this methodology applying it to the North Sea pelagic ecosystem. We use the target species Calanus finmarchicus, a key copepod in temperate and subarctic ecosystems, along with 86 biological, hydrographical and climatic time series, ranging from local water nutrients and fish predation, to large scale climate pressure patterns. The chosen study area is the central North Sea, from 1972 to 2011, during which period there was an ecological regime shift. The GP based analysis identified 3 likely drivers of C. finmarchicus abundance, which highlights the importance of considering both physical and trophic drivers: temperature, North Sea circulation (net flow into the North Atlantic), and predation (herring). No large scale climate patterns were selected, suggesting that when there is availability of both data types, local drivers are more important. The results produced by the GP based procedure are consistent with the literature published to date, and validate the use of GP for interpreting species dynamics. We propose that this methodology holds promises for the highly non-linear field of ecology.

Introduction

Ecosystems dynamics are an integrated response of the ecosystem’s biological components (species/groups) to drivers, which act independently, synergistically or even antagonistically [1, 2]. These drivers are defined as any natural or human-induced factor that directly or indirectly causes a change in an ecosystem or population [3].

A long, unsolved ecological question is whether top-down trophic drivers (i.e. predation) or bottom-up drivers (often intended as climate, hydrography, food) control populations and ecosystem states [4]. With regard to these terms, Conversi et al [2] note that the usage of climate variables as bottom-up, although widely used, may not be entirely correct for the marine environment: in fact, climate related variables, such as temperature and other physical factors are likely to simultaneously affect several trophic levels in the marine food chain [5], or may skip some levels [6], hence are not operating in a strictly “bottom-up” manner. Hence, in this article we use Conversi et al [2] definitions of drivers: “physical drivers”, which are related to the physical environment, and “trophic” (or “biological”) drivers, which are related to predator-prey interactions. Only the latter are subdivided in bottom-up and top-down drivers.

Several studies have suggested that physical processes such as climate-induced temperature or circulation changes are the main drivers for planktonic populations [7-23], whilst others suggest that human induced top-down changes, such as overfishing, which leads to changes in trophic level structure and to trophic cascading, are the main drivers [24-31]. These theories are then confronted by the idea of synergistic relationships, where a combination of drivers control a population and a change in the population depends on how drivers affect each other and/or on the resilience of the ecosystem [1, 2, 4–6, 32–43], and the debate is still open.

In this work, we address the question of the drivers of planktonic populations using a novel, holistic approach, in which multiple, potential drivers are analysed without assumptions on either their relative importance or the mathematical relationship between them and the target species.

These potential drivers are analysed using a symbolic regression methodology based on genetic programming (GP). GP [44-46] is a domain-independent evolutionary computation methodology, capable of generating solutions to a given problem without any strong a-priori knowledge or assumption on the problems solution (see S1 File for details). GP has been used successfully in a wide range of applications [46], such as robotics [47, 48], physics [49], stock market analysis [50], and medicine [51]. Its use for understanding environmental patterns and drivers is very recent, for example, investigating the causes of copepod variability in the English Channel [52], predicting harmful algal blooms [53], predicting early warning of cyanobacterial blooms in freshwater ecosystems [54], and the energy output of wind farms based on weather data [55]. As well as predictions, genetic algorithms have been used to parameterise coupled biological–physical copepod population-dynamics computations [56]. In particular, the GP-based approaches in [54] have been compared with other methodologies for multivariate function fitting and knowledge discovery and have shown superior quantitative performance.

The proposed GP-based analysis is conceptually opposite to the mechanistic approaches where modelling is based on a-priori knowledge and the model parametrizations are based on literature values and laboratory experiments (e.g., [57, 58]). Quite the reverse, the symbolic regression methodology based on GP is capable of extracting mathematical models from the data, whose interpretation can provide useful information on the investigated context. These properties make it worth investigating how this methodology can be applied for analysing ecosystems, in particular for detecting variables that can be relevant for (i.e., explain ecosystem variability of) a target variable, and for identifying the mathematical models governing them [45, 59, 60]. The GP approach can be of particular service in the analysis of population dynamics, since species may respond to environmental change in a nonlinear manner (as in the case of ecological regime shifts), and the usual linear analyses may not be appropriate.

In this work, we aim to validate a novel methodology based on GP at distinguishing which drivers can influence species abundance. We use the target species Calanus finmarchicus in the North Sea, a key copepod in temperate and subarctic ecosystems, along with 86 biological, hydrographical and climatic time series, ranging from local water nutrients and fish predation, to large scale climate pressure patterns.

Calanus finmarchicus in the North Sea

The North Sea is situated on the continental shelf of northwest Europe. It opens into the Norwegian Sea and the Atlantic Ocean to the north, into the Atlantic Ocean to the southwest via the English Channel, and into the Baltic Sea to the east through the Skagerrak Strait. The North Sea is often divided into the shallow southern North Sea (<50 m depth), the central North Sea, the northern North Sea (~200 m depth), the Norwegian Trench (up to 700 metres deep) and the Skagerrak. In the North Sea, North Atlantic water mixes with freshwater runoff and river discharges in a predominantly anti-clockwise circulation. Smaller currents move southwards along the east coast of the UK and northwards along the continental western European coast. Shallow areas of the North Sea (<30 m) are generally fully mixed by tidal action, whilst deeper areas have a surface mixed layer (upper 30 m) which is usually mixed by wind action [61, 62].

The North Sea is an area of high importance for commercial fisheries and wildlife, and once provided 5% of the global fish harvest [41, 63]. However, the North Sea is under increasing pressure from the effects of climate change and associated temperature increase [64] and human populations; with increased coastal development, pollution, nutrient input, maritime transport, as well as food and energy demands with around 30 different commercial fish stocks still being exploited, and interests from gas, oil and renewable energy industries [61, 62]. It has been an area of focus and debate on the drivers of pelagic populations due to a series of abrupt ecosystem transitions, termed regime shifts [2, 65, 66]. These events affected the entire North Sea ecosystem, however, the exact timings and causes are inconclusive, with physical and trophic possible drivers, as well as synergistic combinations, being identified using a variety of methods [21, 32, 66–76].

Copepods constitute a key trophic group, transferring energy from phytoplankton to higher trophic levels [77]. Members of the genus Calanus are amongst the largest copepods and can comprise as much as 90% of the dry weight of mezozooplankton in regions of the North and Celtic Seas [77].

The calanoid copepod species Calanus finmarchicus is one of the most important large zooplankton species in the subarctic waters of the North Atlantic [78] and dominates the dry weight of the mesozooplankton in the northern regions. C. finmarchicus does not reside year-round in the North Sea (although it did so in the 1960s [57]), and its population is replenished each spring by advection of late-stage individuals mainly originating from an overwintering stock located beyond the shelf edge, in the deep Atlantic and Norwegian Sea, particularly the Faroe-Shetland Channel [79, 80], which makes this species presence in this area very susceptible to circulation changes.

Over the past four decades there has been a pronounced decline in its abundance in the North Sea: in 1962, the species represented 80% of the total Calanus, and it used to dominate the spring time biomass, whereas it represented only 20% of the genus by the beginning of the 2000s [8, 11, 12, 81, 82].

Calanus finmarchicus is subjected to a range of drivers, both physical and trophic (top-down and bottom-up) (Fig 1) [81]. From eggs to adults, it is a primary source of food for commercially important fish species, such as cod (larvae), herring, lesser sandeel, mackerel, blue whiting, anchovy and chaetognaths, and as such its dynamics have been extensively studied [28, 63, 83–89]. C. finmarchicus feeds on a range of microplankton including diatoms, dinoflagellates, ciliates, coccolithophores and rotifers [90, 91]. Sea surface temperature determines the geographical, vertical, seasonal distributions, the development, and physiology of calanoid copepods, and can be particularly relevant for C. finmarchicus, as this species is at the edge of its thermal niche in the North Sea [11, 23, 57, 58, 71, 73, 77, 92–95]. Because this species reproduces in the Norwegian Sea and individuals are transported from there into the North Sea, marine circulation and the atmospheric patterns driving it are crucial for its geographical distribution [11, 22, 79, 80, 96].

Fig 1

Variables identified as potential drivers of the abundance of C. finmarchicus.

Potential drivers have been sectioned into physical, on the left, and biological, on the right. Biological variables have subsequently been divided into two further groups, top-down and bottom-up, which are positioned above and below C. finmarchicus respectively. The data used in this research article are listed in Table 1.

Table 1

Timeseries used in this study.

The links to the data sets used are shown in S3 File.

Dataset		Units	Period	Gaps	Frequency	Area	Source
Physical variables	North Atlantic Oscillation (NAO)		1950–2012	No	Monthly	North Eastern Atlantic	NOAA-CPC
	East Atlantic Pattern (EA)		1950–2012	No	Monthly	North Eastern Atlantic	NOAA-CPC
	East Atlantic West Russia Pattern (EAWR)		1950–2012	No	Monthly	North Eastern Atlantic	NOAA-CPC
	Scandinavian Pattern (SCA)		1950–2012	No	Monthly	North Eastern Atlantic	NOAA-CPC
	Polar Eurasia Pattern (POL)		1950–2012	No	Monthly	North Eastern Atlantic	NOAA-CPC
	Atlantic Multidecadal Oscillation (AMO)		1948–2012	No	Monthly	North Atlantic	NOAA-PSD
	Northern Hemisphere Temperature (NHT)		1850–2012	No	Monthly	Northern Hemisphere	NOAA-anomalies
	N. Atlantic Southward Flow	Sverdrup	1970–2012	No	Monthly	Orkney-Norway	NORWegian ECOlogical Model System (NORWECOM)
	N. Atlantic Northward Flow	Sverdrup	1970–2012	No	Monthly	Orkney-Norway	NORWegian ECOlogical Model System (NORWECOM)
	N. Atlantic Net Flow	Sverdrup	1970–2012	No	Monthly	Orkney-Norway	NORWegian ECOlogical Model System (NORWECOM)
	English Channel Eastward Flow	Sverdrup	1970–2012	No	Monthly	Dover Strait	NORWegian ECOlogical Model System (NORWECOM)
	English Channel Westward Flow	Sverdrup	1970–2012	No	Monthly	Dover Strait	NORWegian ECOlogical Model System (NORWECOM)
	English Channel Net Flow	Sverdrup	1970–2012	No	Monthly	Dover Strait	NORWegian ECOlogical Model System (NORWECOM)
	Sea Surface Temperature (SST)	°C	1891–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (surface data)
	Sea Surface Salinity (SSS)	PSU	1891–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (surface data)
Biological variables
Bottom-up drivers	Total Nitrogen (N)	μmol/l	1970–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (Bottle data)
	Total Phosphorus (P)	μmol/l	1969–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (Bottle data)
	Silicate (SiO4)	μmol/l	1958–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (Bottle data)
	Chlorophyll-a (Chl-a)	μg/l	1961–2012	Yes	Monthly	55 to 60°N and -2.5 to 9°E	ICES (Bottle data)
	Phytoplankton Colour Index (PCI)	1–5 scale	1958–2011	Yes	Monthly	CPR areas C1 and C255 to 58°N and -3 to 11°E	Sir Alister Hardy Foundation for Ocean Science
Top-down drivers	Chaetognaths Eyecount	Mean	1958–2011	Yes	Monthly	CPR areas C1 and C255 to 58°N and -3 to 11°E	Sir Alister Hardy Foundation for Ocean Science
	Total Fish Larvae	abundance	1958–2011	Yes	Monthly	CPR areas C1 and C255 to 58°N and -3 to 11°E	Sir Alister Hardy Foundation for Ocean Science
	Herring (Claupea harengus) Total Stock Abundance	estimate	1947–2011	No	Annual	Subarea IV	ICES HAWG
	Herring (Claupea harengus)Total Stock Biomass	tonnes	1947–2011	No	Annual	Subarea IV	ICES HAWG
	Cod (Gadus morhua) Aged 1 year	‘000’s	1963–2011	No	Annual	Subarea IV and Divisions IIIa and VIId	ICES WGNSSK
	Cod (Gadus morhua) Spawning Stock Biomass	tonnes	1963–2011	No	Annual	Subarea IV and Divisions IIIa and VIId	ICES WGNSSK
	Target Variable:
	Calanus finmarchicus	Mean	1958–2011	Yes	Monthly	CPR areas C1 and C255 to 58°N and -3 to 11°E	Sir Alister Hardy Foundation for Ocean Science

While the literature cited above indicates that this species can potentially be influenced by a variety of drivers, many of the studies mentioned use only part of the possible drivers, which might result in a partial picture [97]. Whilst most research traditionally focuses on a single driver (e.g. NAO, SST, or fishery impacts), in this work we use the GP-based methodology on a large (86) collection of variables, encompassing both physical and trophic, bottom-up and top-down, potential drivers of C. finmarchicus (Table 1), without any a priori assumptions on the population’s drivers.

Materials and Methods

Overall analysis approach

The statistical approach used in this work is based on Genetic Programming, complemented by additional statistical analyses (Fig 2). The Genetic Programming approach is a data driven methodology, which has no a priori assumptions on the relationship between the data—these relationships are selected through an evolutionary process lasting hundreds of generations. This evolutionary/selection process is validated through a cross validation procedure. The Relevance Analysis complements the selection of the variables by identifying those that occur with a frequency superior to chance, which are therefore deemed relevant for the approximation of the target variable. The Gradient Analysis is specific to the identification of the drivers of the target variable since it provides essential information for causal associations, i.e., the direction (direct or inverse) between target and relevant variables. The CuSUM analysis identifies the year of a regime shifts, and is specific to the North Sea case study, as this sea was involved in an ecological regime shift in the 1980s.

Fig 2

The analysis approach used in this work.

The figure summarizes the proposed analysis approach, where TS1, …, TSn are the time series shown in Table 1, Norm is the data normalisation step needed to analyse time series with different magnitudes, the Genetic Programming and the Cross-Validation respectively generate and validate the functions that approximate the target variable, the Relevant Analysis identifies the relevant ecological variables and the modelling functions capable to express the target variable, and finally the Gradient Analysis identifies the role of the relevant variables in relation to the target variable. CuSUM is the cumulative sums analysis used to identify the starting year of the regime shift in the target and relevant variables.

The details of each analysis are given in the next sections. The mathematical in depth information is provided in S1 File.

Genetic programming

Genetic Programming (GP) is an evolutionary computation approach that generates solutions starting from an initial population of randomly generated functions, based on a set of variables, mathematical primitives, and constants [44-46]. As the name suggests, the initial solutions are improved by miming the selection processes that occur naturally in biological systems. This is done through an iterative process encompassing multiple generations, during which the Selection, Crossover, Clonation, and/or Mutation GP operators are applied.

This process is explained in detail in S1 File. The initial population consisted of 1000 individuals (i.e., functions) and the evolutionary process lasted for at most 500 generations, when the stop criterion was met. To evaluate the fitness of the evolved individuals we used the Root Mean Squared Error (RMSE) [98] between the target and the evolved individuals.

Similarly to the approach described in [99], the GP based procedure was coupled with a Cross-Validation framework (CV) [98, 100, 101], with the aim of selecting a set of approximating functions that generalise the target variable’s dynamics.

The CV process randomly splits each time series into two disjoint subsets: (i) the training subset used to run the GP procedure in order to evolve the functions which approximate the target variable; (ii) the validation subset used to evaluate the approximation capability of the evolved functions. Since the time indices of the validation subset are not involved in the training phase (training and validation subsets being disjoint), the generalization capability of an evolved function is defined as the error obtained by predicting the values of the target variable on the time indices of the validation subset. The smaller the prediction error, the more the function generalises and is therefore more suitable for explaining the dynamics of the target variable. The best generalising functions are selected to become part of a group of functions approximating the target variable, named population-pool. The heterogeneous mathematical forms of the population-pools individuals are normal in the context of GP, as the diversity of the evolved genetic material (genotype) captures the complexity of the target variable [102]. In fact, functions with different genotypes (function syntax) can produce similar phenotype effects (function semantic), as genotype diversity, together with phenotype convergence, indicates the existence of multiple solutions for the same problem [102].

All the population-pool functions were subsequently analysed for identifying the variables relevant to the target time series.

Relevance Analysis

A variable is deemed if it appears in the individuals of the population-pool more times than by chance. On the contrary, not relevant variables should be considered randomly selected by the GP procedure and for this reason they are not suitable for explaining the behaviour of the target variable.

The problem of identifying whether a variable is relevant or not can be formalised as a Bernoulli trial [103] defined on the number of functions of the population-pool. In the proposed experiments, the assumption that all the variables have the same probability to appear among the individuals of the population-pool is rejected with p-value equal to 0.001 according to Johnson [104]. According to the proposed Bernoulli trial, this p-value corresponds to a variable occurrence greater than or equal to a threshold variable occurrence. This means that in our study all the variables whose occurrence among the individuals of the population-pool is greater than or equal to the threshold are deemed relevant (for details see S1 File).

Some approximating functions, or individuals, in the population-pool contain only relevant variables; other functions contain both relevant and not relevant variables. The former individuals should be considered a reliable set of functions for modelling the target variable, since they do not contain random components; for this reason, in the remainder of this article the individuals that contain only relevant variables are named (listed in S1 File). In contrast, the individuals containing random components (not relevant variables) are deemed not valid for generalising the dynamics of the target variable [52].

The relevant variables and the modelling functions extracted from the population-pool through the k-fold cross-validation framework and the proposed statistical analysis reduce the possibility that random information contained in the evolved functions produces over-fitting of the investigated time-series.

Gradient Analysis

The GP procedure, as all multivariate or correlation analyses, can identify the relationships but cannot provide information on the corresponding causal direction, i.e., it does not have the ability to distinguish the internal dynamics in trophic chains, in particular drivers from “drivees”. This distinction needs to be made through the analysis of the relationships between approximating functions, original time series, and target variable [105].

This is a common problem in the analysis of populations within food webs. To solve this we have added to the results of the GP procedure a gradient analysis with the purpose to identify the role of each relevant variable with respect to the target variable dynamics. In this work, the partial derivatives of all the modelling functions have been analysed based on the finite difference approach in order to estimate the function gradient, and from this, to infer whether a relevant variable is directly or inversely proportional to the target variable.

Cumulative Sums Analysis

The relevant variables were tested for a change point in their mean to further compare their significance to the target variable. Cumulative Sums analysis (CuSUM) was used to determine the year of the most prominent shift in the time series. This method consists of plotting the cumulative sum of standardised values over time. Each value of the series is subtracted from the mean of the time series, resulting in a new time series of residuals which are used for the calculation of the cumulative sum. The interpretation is based on the slope of the line on the chart: a constant deviation from the mean of the time series shows a constant slope. Persistent changes from the mean of the times series cause a persistent change of the slope [66, 106].

Case Study: Calanus finmarchicus in the North Sea

Data

A total of 26 environmental variables representing large and local scale, trophic or physical potential drivers of C. finmarchicus were used in this study (Table 1). Climatic and hydrographic variables were further divided into seasonal averages because we wanted to evaluate the effect different seasons can have on zooplankton abundance [95]. Therefore, a total of 86 time series were used in this experiment. The multiplication of time series over several seasons is not considered as a duplication of information, as the GP based approach simply identifies the variables that are relevant in relation to C. finmarchicus and ignores all other variables, independently from the number of variables used.

Plankton data were obtained from the Sir Alister Hardy Foundation for Ocean Science (SAHFOS) database for the Continuous Plankton Recorder (CPR), survey regions C1 and C2, located in the central North Sea. The CPR survey method is described in detail in [107-110]. The CPR plankton data (abundances) were separated into the target variable (Calanus finmarchicus), bottom-up (phytoplankton colour index (PCI)—proxy for food availability), and top-down (chaetognaths and total fish larvae—proxies for predation) trophic variables.

Herring (Clupea harengus) and cod (Gadus morhua) estimates were chosen because they are major predators of C. finmarchicus [63, 84, 111]. Yearly fisheries stock assessment data for herring total stock biomass (TSB), herring abundance (sum of all age groups), cod spawning stock biomass (SSB), and 1 year old cod (cod1) were retrieved from the annual stock-assessment reports for the North Sea, ICES area IV [112, 113], which corresponds to the area selected for C. finmarchicus. Cod spawning stock biomass values were used as a proxy for the amount of larvae (which feed on C. finmarchicus) which could be produced by the adult cod that year, with a higher SSB inferring more cod larvae production. Cod stock assessments do not estimate year 0 age groups, therefore estimates of 1-year-old cod were also used as a proxy for the abundance of larvae present the year before. These two variables were used to represent the predation pressure of larval cod, as that is the key life stage which feeds on plankton.

Oceanographic data were downloaded from the International Council for the Exploration of the Sea (ICES) Oceanographic Database for the same area as CPR survey regions C1 and C2. Chemical and biological variables, total nitrogen, total phosphorus, silicate (SiO4), all as μmol/l, and chlorophyll-a (μg/l) were downloaded from the CTD and Bottle Data database (Table 1). Total nitrogen, phosphorus and silicate are considered proxies for primary productivity.

Seawater temperature and hydrography affect C. finmarchicus biogeography [39, 64]. SST (°C) and surface salinity (PSU) were downloaded from the ICES Surface Data database, which consists of CTD, Bottle and Underway/Pump data collected at depths < 10m (Table 1).

Additional temperature measurements included the Northern Hemisphere Temperature (NHT) index, which is a large scale temperature index, defined as the combination of land and sea surface temperature anomalies over the entire northern hemisphere and the Atlantic Multidecadal Oscillation (AMO), a regional index that describes long-term variations in sea surface temperature of the Atlantic Ocean [114] (Table 1).

Large-scale climatic patterns influence circulation in the North Atlantic and North Sea and therefore need to be considered as potential drivers of C. finmarchicus abundance. For this work, the following climate patterns where chosen: North Atlantic Oscillation (NAO), East Atlantic Pattern (EA), East Atlantic/ West Russia Pattern (EAWR), Scandinavian Pattern (SCA), The Polar/ Eurasian pattern (POL). These data were downloaded from http://www.cpc.ncep.noaa.gov/data/teledoc/telecontents.shtml (Table 1).

North Sea circulation and flow data are from the NORWegian ECOlogical Model system (NORWECOM). NORWECOM is a coupled 3D physical, chemical and biological model, validated for the North Sea and the Skagerrak [115]. Average monthly transports through an east-west section from Utsira (Norway) to the Orkney Islands along 59°17’N in the northern North Sea (N. Atlantic inflow) and a longitudinal section through the English Channel in the Dover Straits along 0°E (English Channel inflow) were computed, and values for inflow, outflow and net flow were used. Inflow represents southwards flow through the Utsira to Orkney transect and eastward through the Dover Strait, whilst outflow represents northward and westward flow respectively.

The data used shows some bias which should be mentioned. It has been suggested that the CPR shallow sampling depth may not be adequate at sampling species absolute numbers [116, 117], and those which migrate beneath the thermocline, such as C. finmarchicus [118]. However, analysis has shown the CPR is a consistent semi quantitative index of phyto- and zoo-plankton abundance in comparison to satellite data [119] and to other planktonic sampling techniques [107, 109, 120, 121], and in particular that it provides an accurate representation of both spatial and temporal patterns of C. finmarchicus [122].

Fishery stock assessment data are based on model estimates from all available data for that species and are reviewed extensively before publishing; they are favoured over total landing values because they take into account fishing effort [123]. The yearly-averaged nutrient data may not contain 12 monthly values, meaning some months may be underrepresented, especially at the start of the series. Cod SSB and cod1 are not actually feeding on C. finmarchicus larvae (see Discussion), and they were used as proxies for larval production. These two variables were chosen because they were the best larval cod proxy data that we could find.

Data manipulation

Due to missing data and different begin-end dates in the time series, this study used the period in common, 39 years from 1972 to 2011. Yearly averages which consisted of 3 or less monthly values were excluded from the GP-based model to improve accuracy. The removal of these years is the cause of the gaps in the approximating functions output.

The 39 time indices were split into k = 13 folds corresponding to 3 years each according to the k-fold Cross-Validation design, and the GP procedure was run 50 times per fold, each time randomly sampling 85% of the time indices.

ICES chemical and biological data recorded below the depth of 20m were excluded from analysis, due to the towed depth of the CPR and the productive surface layer above the thermocline generally forming at 15–20 m in the summer months in the North Sea [68]. Fisheries data were available at an annual frequency whilst all other variables were downloaded at or averaged to a monthly frequency. Because of this, all monthly times series were yearly or seasonally averaged prior to analysis. Seasonally averaged data means that the averages were calculated over the 3 months of each season, as follows. Climate seasons for NAO, EA, EA.WR, SCA, POL, AMO and NHT, were defined as winter (December-February), spring (March-May), summer (June-August) and autumn (September-November), whilst seasons for hydrographic variables (Atlantic flow, English channel flow, SST and salinity) lag one month behind climatic, i.e. winter is January-March, etc. [124]. Nutrient and biological data were averaged annually with no seasons implied.

The time series presented in Table 1 refer to the observed data with their own units and scales. All the time series were normalised by dividing each value by the maximum value of that time series [125]. The time series having only positive values were normalised in the range [0,1], while the time series having both positive and negative values were normalised in the range [-1,1]. Without this normalisation approach, time series whose values are in the order of 107 (e.g., cod and herring abundance) could not be compared to time series whose values were in the order of 101 (e.g., hydrographic parameters, climate indices). This normalisation approach discards the scale of the time series, and allows focusing only on the time series dynamics. In this way, the proposed analysis approach can select the relevant variables capable to describe the dynamics of the target variable.

Results

Application of the GP-based analysis to Calanus finmarchicus abundance in the North Sea

1. Identification of the relevant variables

In this work, we have used GP to identify the variables that are relevant for approximating the abundance of Calanus finmarchicus in the North Sea, using 86 biological and physical variables, chosen because they are potential drivers of C. finmarchicus. The GP parameters and the Cross validation design used in the case study are detailed in S1 File.

At the end of all the GP/cross-validation iterations, the population-pool consisted of 104 functions. The average validation error of these functions, resulting from the cross-validation process, was 0.102 with a standard deviation equal to 0.082.

The relevance analysis of the population pool identified the variables that occurred with a frequency superior to chance among the approximating functions of the population-pool. According to the statistic test based on the Bernoulli trial (see S1 File for details), the variables whose occurrence frequency was larger than 7 were deemed relevant. Out of the original 86 potential drivers, 9 variables, statistically relevant for the approximation of C. finmarchicus, where thus identified. These 9 relevant variables (6, once seasons and multiple datasets are combined, e.g., herring biomass and abundance represent 1 potential driver, herring) are, in order of frequency of occurrence in the final population-pool: herring (total stock biomass and total abundance estimates), cod (spawning stock biomass and abundance at age 1), phytoplankton (as represented by the PCI), SST (spring and winter), and two circulation variables, N. Atlantic net flow in winter (wNAtlNET) and English Channel eastward flow in summer (smEnglChanE) (Table 2, Fig 3). It is worth noticing that a) the relevant variables encompass both top-down/ bottom-up (predators, prey proxies) biological and physical (climate, temperature, circulation,) potential drivers, and b) they do not encompass large-scale climate indices.

Table 2

Relevant variables selected by the Genetic Programming based methodology combined with the relevance analysis, the abbreviations used in this article, and the frequency of occurrence of each variable in the 104 approximating functions in the population pool.

The last column indicates the type of relationship between the relevant variables and C. finmarchicus, identified with the gradient analysis.

Variable:	Short name	Occurrence	Direction of relationship
Herring Total Stock Biomass	HerringTSB	40	inverse
Cod Spawning Stock Biomass	CodSSB	38	variable
Phytoplankton Colour Index	PCI	35	inverse
Herring Total Abundance	HerringTAE	15	inverse
Cod Aged 1	Cod1	13	direct
Winter North Atlantic Net Flow	wNAtlNET	9	direct
Spring SST	spSST	8	inverse
Winter SST	wSST	7	inverse
Summer English Channel Eastwards Flow	smEnglChanE	7	inverse

Fig 3

Time series of Calanus finmarchicus and the 9 relevant variables identified by the GP procedure.

The time series are ordered from left to right of most frequently occurring. All time series were normalised by dividing each value by the time series maximum value before the GP process.

2. Temporal shifts in the relevant variables

In order to further evaluate the selection of the 9 variables by the GP procedure, each variable was tested for an abrupt shift in the mean using Cumulative Sums analysis. All 9 time series showed a shift during the 1980s (Fig 4) (see S2 File for each variables CuSUM plot).

Fig 4

Sequences of abrupt shifts in the North Sea.

Time series of C. finmarchicus average annual abundance, with arrows indicating the years of the regime shifts, detected using cumulative sums analysis both on this species and in the 9 variables identified as relevant by GP. The table in the insert specifies the year of the shift and shows its direction: + meaning an increase,—meaning a decrease.

A shift towards a lower mean abundance of C. finmarchicus was detected in 1986. This happened after shifts towards higher means in herring (1982/85), summer English Channel eastward flow (1984) and winter N. Atlantic net flow (1985). After the downward shift in C. finmarchicus, PCI saw an increase in 1987, while cod (both age 1 and spawning stock biomass) concurrently decreased. Winter and spring SST shifted upwards in 1988 and 1989 respectively.

3. Gradient analysis of the relevant variables

In order to understand the nature of the relationship of the 9 relevant variables with C. finmarchicus, gradient analysis was used. This analysis provides information on the type of relationship (direct, inverse, variable) between each relevant variable and C. finmarchicus, which gives important information for identifying the likely drivers.

The results of this analysis are presented in the last column of Table 2. All relevant variables are inversely related to C. finmarchicus, with the exception of N. Atlantic net flow (positively related), and cod, which is either positively related (cod1) or variable (SSB).

This analysis indicated that some of the relevant variables are very unlikely drivers of C. finmarchicus.

PCI showed an indirect relationship with C. finmarchicus. In the studies period PCI (a prey proxy) increased, yet C. finmarchicus decreased. Further, the CuSUM analysis showed a (positive) shift in PCI after a (negative) shift in C. finmarchicus. Both analyses suggest that PCI is not driving, but, if anything, is driven by C. finmarchicus, but see the Discussion. Hence, PCI has been eliminated as a likely driver.

Cod age 1 showed a direct relationship with C. finmarchicus. This predator (proxy) decreased in the studied period, yet C. finmarchicus decreased as well, indicating that it was not driven by it (but see the Discussion for more information). Hence Cod age 1 was eliminated as a likely driver.

Cod SSB had a variable relationship showing no consistent relationship with C. finmarchicus, hence it was also eliminated as a likely driver, but see the Discussion for more details on this variable.

The summer English Channel eastward flow was eliminated as a potential driver over magnitude considerations. The overall flow exchange with the N. Atlantic through the Strait of Dover (English Channel flow) is about 1/10th of the flow exchange through the Faroe-Shetland transect (N. Atlantic flow) [126]. However, due to the normalisation process applied in this work, the differences in the variables magnitude are lost, and all variables belong to the same range. While this variable resulted relevant, its small magnitude indicates that it is not a likely driver.

All the above variables remain nevertheless relevant and can be used to approximate C. finmarchicus abundance and for forecasting experiments,

With the use of gradient analysis, supported by a literature review, 5 (3 once combined) of the 9 relevant variables selected by the relevance analysis can be identified as likely drivers: herring (total abundance estimate and total stock biomass), SST (winter and spring) and N. Atlantic net flow in winter, shown in Fig 5.

Fig 5

The new conceptual model of potential drivers of Calanus finmarchicus abundance deduced from the GP method supplemented by the relevance and gradient analyses.

The new model is composed of two physical variables, North Atlantic net flow and Sea Surface Temperature, and one biological variable, Herring, which is recognised as a top-down driver.

Discussion

1. Evaluation of the GP-based methodology

The aim of the proposed analysis methodology is the identification of the relevant information needed to explain the dynamics of the investigated ecosystem. This is obtained by evolving a set of mathematical functions capable to generalise the behaviour of the target variable (i.e. C. finmarchicus abundance in the case study).

In the design of the statistical methodology used in this work we took special care in producing generalized (vs overfitting) results, in significantly selecting the relevant variables, and in correctly identifying potential drivers.

The generalising capability of the evolved functions was obtained by coupling the Genetic Programming procedure with a k-fold Cross-Validation (CV) design, which provided a prediction error for each evolved function. The individuals with the smallest prediction error (i.e., the greatest generalising capability) were selected for the final population pool.

The relevance analysis provided the tool to statistically identify within the population pool relevant (not occurring by chance) variables for the target species. Additional experiments were ran, which used randomly generated time series together with the ecological variables. During these experiments, the artificial variables were never selected as relevant, thus confirming the reliability of the proposed statistical analysis.

While the GP/cross validation plus relevance analysis methodology can identify significant relevant variables (in the case of C. finmarchicus, nine out of the 86 initially selected possible drivers), a present limitation of this methodology is that it cannot provide information on the role they play for explaining the dynamics of the target variable. This is a limitation common to all multivariate or correlation analyses. In the C. finmarchicus application, we approached the question of the causal relationship among relevant and target variables with the gradient analysis, by analysing the partial derivatives of the modelling functions, i.e., the members of the population-pool involving only relevant variables. The Gradient Analysis provided information on the relationship between the target variable and the relevant variables incorporated in the modelling functions.

The results achieved by applying this methodology to Calanus finmarchicus in the North Sea are coherent with what can be seen in the literature (see next section), which validates this approach.

2. Application to Calanus finmarchicus in the N. Sea

In this section, we combine the results of GP/CV, relevant, gradient and cumulative sums analyses with a literature review, and discuss their relationship with C. finmarchicus to explain its decline over a 40 year period. The set of variables selected by the GP analysis is justified and evaluated in order to demonstrate the potential use of this method in identifying drivers within a marine system.

Relevant variables

The 9 variables shown in Table 2 and Fig 3 encompass both physical (temperature and circulation) and biological variables (both top-down (predators) and bottom-up (prey) proxies). All the relevant variables are associated to modelling functions (listed in S1 File), which can be used to approximate or forecast the target variable’s abundance.

Relevant physical variables

Out of the large selection of potential physical drivers (Fig 1 and Table 1), only SST and circulation proxies (the flow in/out of the North Sea) were identified as relevant. It is quite interesting that none of the large-scale climate patterns (e.g., NAO, NHT…) were identified by the GP model as relevant variables. These findings contradict a large body of research cited earlier suggesting that large-scale climate anomalies are the main drivers of marine ecosystems. However, our results may have been the outcome of having both large scale and local drivers as model inputs. In other words, when there is availability of both data types, local drivers seem to be more important factors to population dynamics, even if they are controlled by larger scale factors.

Winter and spring SST (Fig 3H & 3I) increased steadily over the time series, and the gradient analysis identified an inverse relationship with C. finmarchicus (Table 2), meaning that a temperature increase corresponded to a decrease in this species. It has been suggested that the increased SST in the North Sea has affected the physiology and development rates of C. finmarchicus individuals [58], since this species, whose thermal critical boundary is around 10° C, is at the edge of its niche in the North Sea [81, 127]. The temperature increase is proposed to be the main cause of a northwards range shift of this and other species, and of its decrease in the study area by several authors [2, 5, 8, 14, 19, 81, 128, 129], although other studies have found C. finmarchicus is still abundant below the thermocline [118], indeed that under higher SST the relative abundance of C. finmarchicus in deeper water increases [57].This suggests that the North Sea population might have moved to a cooler, deeper zone as SST increased, reducing its abundance in the CPR time-series. While increasing spring SST may affect the development and growth of C. finmarchicus [58, 130, 131], changes in winter SST may also be associated to circulation changes, which also can affect plankton abundance in the sea [15, 132]. SST in fact acts as an overarching driver, not only directly affecting C. finmarchicus, but also its predators, its prey [105, 133, 134] and the surrounding habitat. For example, an increase in winter SST is suspected to affect the overwintering habitat of C. finmarchicus in areas such as the Norwegian Sea deepwater. An increase in SST and deeper ocean temperatures reduces the available area for C. finmarchicus to overwinter causing a decrease in the North Sea population by affecting reproductive capabilities and leading to less individuals being advected into the North Sea in spring [79, 96].

N. Atlantic is a circulation proxy for the net amount of flow leaving-entering the North Sea via a transect running from Scotland to Norway. The winter N. Atlantic net flow was northward (positive) throughout the time series, as the flow leaving was always greater than that entering the N. Sea. The flow intensity varied over time, increasing until 1990, and then declining (Fig 3G). The C. finmarchicus time series shows an opposite trend, steadily declining since the late 1970s, then staying low before a year of higher abundance at the end (Fig 3), and in fact the relationship between these variables is overall inverse (Table 2). This inverse relationship between C. finmarchicus and net flow can be related to the transport of zooplankton into/out of the North Sea. The increasing overall flow to the north in the first half of the series could in fact have significantly contributed to reducing the number of C. finmarchicus that were re-entering the North Sea [11, 12, 22]. An increase in winter northward flow could also transport warmer water further northwards which could contribute to the northward shift of this cold-water species and its reduction in the North Sea as well as impacting the size of cold water overwintering areas such as the Norwegian Sea deepwater [8, 96]. It has also been suggested that transport of warmer water into the North Sea has led to increased competition for food from warm water species such Calanus helgolandicus [77].

The English Channel flow represents the flow entering/exiting the southern North Sea via the Strait of Dover (eastward = entering). The summer English Channel eastward flow (Fig 3J) has been eliminated as a likely driver due to the evaluation of its magnitude and potential influence on the North Sea (see Results). It remains, however, a relevant variable for approximating C. finmarchicus abundance.

Relevant biological variables: bottom–up

The bottom-up drivers included chlorophyll, PCI, and nutrients. Of these, only PCI appeared as relevant, but an unlikely driver of C. finmarchicus.

The Phytoplankton colour index (PCI) provides an index based on the colour that the accumulation of green chlorophyll pigments gives to the CPR filtering silk. It is considered to be a semi-quantitative measurement of phytoplankton abundance [120], which is a main food source for C. finmarchicus [90, 91]. PCI saw an abrupt increase during the 1980s (Fig 3D) despite concurrent decreasing nutrient concentrations [135], possibly because of changes in the phytoplankton predator community over the past 40 years [7]. With an increase in phytoplankton one would expect an increase in the abundance of their predator C. finmarchicus. However, our gradient analysis shows the opposite, i.e. a direct relationship (Table 2). This confirms the findings by several studies in the North Sea on increased phytoplankton abundance while C. finmarchicus decreased [19, 70, 136–138]. This discordance suggests that while the feeding environment in the North Sea may have changed to one seemingly more favourable for C. finmarchicus, changes in other drivers might be more important for this particular species. Alternatively, it can be considered that PCI might not reflect C. finmarchicus diet, being an aggregate index which provides no information on species (for example, the abundance of several diatom species decreased in this period [139]). Finally, PCI might have increased because of reduced pressure by C. finmarchicus. Both our analyses and the literature review indicate that it is unlikely that PCI is a direct driver of C. finmarchicus (if anything, the direct relationship with C. finmarchicus suggests that PCI is driven by it), for which reason it was discarded as a driver. However, because of its possible links with overarching physical drivers, such as SST, whose increase could boost phytoplankton biomass but adversely affect C. finmarchicus abundance [7, 136], PCI remains very relevant for approximating and possibly predicting C. finmarchicus abundance, as indicated by its frequency in the GP output.

Relevant biological variables: top-down

Only herring and cod where selected as relevant by the GP approach, while other predators, such as chaetognaths or fish larvae were not selected at all. These two species were the most frequently selected variables, which suggest they are very relevant for C. finmarchicus abundance; however, the overall patterns are very different, and only herring can be considered as a driver.

Herrings are planktivorous throughout its life, and C. finmarchicus can be an important prey [84, 140]. Herring biomass (TSB) is the most frequently selected variable by the GP program, occurring 40 times, and also the selection of herring abundance (TAE) is very high, occurring 15 times (Table 2). The gradient analysis shows that both herring parameters are inversely related to C. finmarchicus (Table 2), which is consistent with the hypothesis that herring can drive C. finmarchicus abundance through predation. This inverse relationship is likely related to the large herring population increase over the time series whilst C. finmarchicus declined at a similar rate (Fig 3B & 3E). The persistent selection of herring by the GP model suggests that herring predation pressure could be one of the most important drivers of C. finmarchicus.

Cod spawning stock biomass (SSB) and cod age 1 abundance (cod1) were also selected as relevant by the algorithm with high frequencies (38 and 13 times respectively, Table 2). The cod variables in fact showed a similar rate of decline as C. finmarchicus over the time series (Fig 3C & 3F). The cause of the decline in cod and cod larvae has been attributed to overfishing and the miss timing of larvae cod and their food, as C. finmarchicus declined and was replaced by other warmer water species of zooplankton which have different life cycles and nutrition quality [36, 82, 133, 134], and to more complex dynamics involving competition with herring, internal feedbacks in the ecosystems, and the reversals of predator-prey roles [43, 74, 141]. Moreover, in this study we used cod SSB and cod age-1 as proxies for cod larvae, for which there were no data available. This is however a gross simplification; in fact, while C. finmarchicus is a major prey for cod larvae, juvenile and adult cods prey on a wide variety of epifaunal species or species associated with hard substrates, and also on the planktivorous herring and so do not directly drive C. finmarchicus [133, 142]. The gradient analysis harmonises with these studies by showing a direct relationship between this copepod and cod aged 1, and a positive but variable (not constant through the series) relationship with cod biomass (Table 2). Both the analyses and the literature review propose that cod is not driving C. finmarchicus. If anything, it is driven by it. Cod SSB and cod-1 still remain very relevant variables for approximating C. finmarchicus abundance, and can possibly be used for predicting it.

Overall, the literature review supports and gives an ecological explanation to the selection of the environmental variables relevant for C. finmarchicus abundance

Temporal changes in the North Sea (the 1980s regime shift)

Abrupt shifts in the means of the relevant variables were analysed with the CuSUM analysis, in order to determine the temporal sequence of changes in the identified variables and thus propose an explaination of C. finmarchicus’s decline (Fig 4). All variables presented a shift in the 1980s. This is a well-known period in time and a large number of studies have addressed the ecological regime shift that involved the North Sea during this period [67, 70] as well as other European basins [15, 68, 69, 124].

It has been proposed that the overfishing of cod caused a reversal in the trophic roles between herring and cod [43, 74]. Herring was released from predation pressure when cod numbers declined and the herring population was then able to increase, shifting around 1982–5 (Fig 4). Herring predation on larval cod subsequently increased, stopping the trophic roles from switching back, even if there was reduced fishing pressure on cod [74]. The increased herring population also predated on zooplankton at a greater magnitude than cod, as herring are planktivorous throughout their lives [63, 84]. This was then intensified by a positive shift in the N. Atlantic net (overall northward) flow in 1985 (Fig 4) which might have contributed to reduced advection of C. finmarchicus from the Norwegian Sea into the North Sea, as well as an increase in central North Sea SST by transporting more heat northwards, possibly impacting the formation of overwintering habitats. The synergistic impacts of these drivers might have caused the abundance of C. finmarchicus to shift to a lower state around 1986 (Fig 4), before an increase (1987) in phytoplankton (PCI), possibly released from the copepod predation or responding positively to increased SST, and a decrease in cod biomass and age 1 larvae (1987), possibly related to their prey reduction and overfishing, occurred (Fig 4). The shift in C. finmarchicus preceded the increase in winter and spring SST in 1988/89, both of which might have contributed to the continuing decline in the abundance of C. finmarchicus in the subsequent years by affecting their geographical range, overwintering habitat space, physiology, food availability and competition with other zooplankton species.

In this proposed scenario, the decline in C.finmarchicus may have been initiated by overfishing and a change in top-down drivers, whose impact has then been magnified by a change in the Atlantic net flow in the North Sea (physical drivers). As the population of C. finmarchicus declined, other drivers, such as the increasing SST, may have caused its abundance to decrease further. This switch in drivers might explain why, while herring have experienced recruitment failures since 2000 [37], the C. finmarchicus population still declined further. The importance of synergistic and time-delayed roles in driving regime shifts is examined in Conversi et al [2], and this work suggests that such a scenario could be identified more frequently than usually thought, provided that the statistical analyses include diverse types of drivers.

Conclusions

This study has shown how the Genetic Programming based methodology can select a small number of variables from a large initial pool. In the application on C. finmarchicus in the North Sea, the ecological meaning of the relevant variables (SST, North Atlantic net flux, English Channel flux, PCI, cod and herring) is backed up by a review of the literature. The analysis of the temporal shifts in these variables can explain C. finmarchicus decline within the 41 years study period.

Moreover, the capability of the proposed analysis approach to produce multiple mathematical models, which use a few variables to approximate the same target, can also be very useful for future robust forecasting applications. These applications could include how the target population respond to climate change using predicted conditions, for example in a global warming scenario. This application could be very important for ecosystem studies and management and we are currently investigating its feasibility. Additional studies of potential interest include comparing this data-driven methodology with knowledge-based mechanistic models.

In the C. finmarchicus application, the list of relevant variables selected by the GP-based analysis, fine-tuned with the relevance and gradient analyses provided 3 likely drivers of this species abundance (Herring, Natl flow, SST) out of the 86 initially selected. This result highlights the importance of both physical and biological drivers on the abundance of this species in the North Sea, and increases our understanding of how climate, circulation, and predation all play a part. The ability to identify relevant drivers indicates future applications for using this GP-based analysis within marine ecosystems, and shows that it can shed light on ecosystem events such as regime shifts.

The application of GP to marine ecology is very novel. With this work we propose that this method holds promises in the near future for the highly non-linear field of ecology, in the same way as the consolidated GP-based methodologies discussed in the introduction have obtained relevant results in many scientific fields, like for example robotics, physics, stock market analysis, and medicine.

Detailed information on the Genetic Programming parameters used in this study.

(PDF)

Click here for additional data file.

Example of an approximating function produced by the Genetic Programming based approach, and cumulative sums analysis results.

(PDF)

Click here for additional data file.

Data availability information.

(DOCX)

Click here for additional data file.