Relating Ctenophore Population to Water Mass Indices in the Northeast U.S. Continental Shelf Ecosystem

Ctenophores exist throughout the Northeast U.S. Continental Shelf Ecosystem, but the underlying mechanisms that control ctenophore populations at this scale are not clear. Ctenophore population data over the last 30 years coincides with changes in several water masses on the shelf, but discovering which water mass was most influential was problematic without mechanistic


Introduction
Large increases in gelatinous zooplankton (e.g., medusae, ctenophores, salps, and siphonophores) have been reported in systems such as the Black Sea and the Mediterranean [1,2], the Caspian Sea [3], and the Bering Sea [4].Increases like these can drive shifts in ecosystem conditions [5].
A recent analysis of observations by NOAA's National Marine Fisheries Service (United States Government) on the Northeast U.S. Continental Shelf between 1980 and 2012 describes spikes in ctenophore abundances [6][7][8].The percent frequency of occurrence of ctenophores, an index developed by Link and Ford [6], spiked from 6 to 42 in the period from 1990 to 1996.
Similar increases were reported for other gelatinous groups including salps and siphonophores [8].Frank [9] described the implications of a spike in the abundance of the ctenophore Pleurobrachia pileus off Southwestern Nova Scotia in 1983.At times when the levels of ctenophores are high, there can be significant ecosystem effects in terms of the consumption of crustacean zooplankton also required by pelagic fish [10].In 1983, Frank reported when concentrations of P. pileus were high, they consumed 26% of the standing stock of copepods in April and 10% in May.This reduction in copepod stock was associated with extremely low levels of haddock (Melanogrammus aeglefinus) larvae produced during the spring of 1983.We feel it is reasonable to assume a similar level of competition between ctenophores and pelagic fish during peaks in our observations.The potential threat to the growth and success of fish populations is an impetus for this study.
We expect the sudden increases and the growth of the ctenophore population on the Northeast Continental Shelf to be associated with changes in oceanographic conditions.The linkage between oceanographic conditions and zooplankton such as copepods is well established and we assume this linkage holds for ctenophores.Inter alia, Greene and Pershing [11], Jossi and Goulet Jr. [12], and Hare and Kane [13] describe these relationships.The oceanographic conditions of the Northeast Continental Shelf can be traced to fluxes in Labrador Subarctic Slope Water (LSSW), Atlantic Temperate Slope Water (ATSW), the position of the north wall of the Gulf Stream, and the phase of the North Atlantic Oscillation.
Potential connections between oceanographic conditions and plankton populations like ctenophores are described in the Ecosystem Status Report [14].Here, changes in LSSW, the position of the north wall of the Gulf Stream, and 2 International Journal of Oceanography the NAO were correlated with a change in the composition of the copepod community on the Northeast Continental Shelf.There was a shift from large copepod species to small copepod species.Since copepods are important prey for ctenophores and ctenophores are not trophically distant from copepods, this correlation supported our expectations that a linkage existed between oceanographic conditions and the ctenophore population.
We attempt to use a new numerical modeling approach to make a first-pass review of the connections that exist between this ctenophore population and the oceanographic conditions of the Northeast Continental Shelf.We will introduce a quantitative technique that will test the fit of thousands of possible models for ctenophore population growth as a function of physical conditions on the shelf.We expect to uncover which oceanographic parameters are influential.By this determination, we expect to make progress toward the development of an indicator of conditions conducive to sudden ctenophore population increase.A leading oceanographic indicator of this type might inform fisheries managers of a shift in the ecosystem.

Methods
We constructed a model of the ctenophore population over time in terms of the physical oceanographic influences of the Northeast Continental Shelf from 1986 to 2008.This unique approach involved three parts.First, the form of the logistic growth model served as the basis for finding functions that relate ctenophore population () and physical conditions.A laboratory study of Pleurobrachia bachei reported that the relationship of body mass over time was sigmoidal [15].This suggested ctenophore growth might be curtailed by finite prey resources.While other investigations have uncovered details of the bioenergetics and growth of ctenophores [16,17], the simpler logistic model appeared to best represent a system where the prey resources of ctenophores are most likely finite.Our model was to serve as a means to evaluate physical oceanographic influences on this population.With the logistic differential equation population model in place, the second step employed a scientific data mining software package, Eureqa [18].Eureqa searches for mathematical patterns within data and has the ability to identify a relationship among parameters within a specified model.Once a formulation for the differential equation was found (and thus the relevant parameters were identified), the third step utilized the capabilities of Mathematica [19] to create numerical and graphical versions of solutions to the differential equation of the model.Completion of this last step allowed us to illustrate how changes in the various indices would impact the ctenophore population as it was modeled.
An estimate of the abundance of ctenophores (percent frequency of occurrence) was obtained through the quantification of stomach contents of the Squalus acanthias, the spiny dogfish (full methodology in [6]).Sampling was conducted across the Northeast Continental Shelf Ecosystem as part of the Food Web Dynamics Program at the NMFS Northeast Fisheries Science Center (NEFSC) in Woods Hole, Massachusetts, USA (http://www.nefsc.noaa.gov/femad/pbb/fwdp/).With this technique, the exact species of ctenophore is not identified and could include species such as Mnemiopsis leidyi, Pleurobrachia pileus, Bolinopsis sp., and Beroe sp.
The percent of Labrador Subarctic Slope Water () and the percent Atlantic Temperate Slope Water () are water mass indices based on oceanographic profile observations near the Northeast Channel at the entrance of the Gulf of Maine from NEFSC surveys.Also, we used another formulation of the percent Labrador Subarctic Slope Water (ℓ) from Mountain [20] based on the same survey data, but including profiles from the Bedford Institute of Oceanography and the World Ocean Database.Since this new formulation of the Labrador Subarctic Slope Water (ℓ) utilized additional inputs, we included it in our analysis alongside the Labrador Subarctic Slope Water from the NEFSC ().The annual anomaly in the position of the north wall of the Gulf Stream () was taken from time series data by Arnold Taylor's Gulf Stream North Wall index (http://www.pml-gulfstream.org.uk/).The annual anomaly in the North Atlantic Oscillation () was the Hurrell Station-Based index taken from differences in atmospheric pressure between Lisbon, Portugal, and Reykjavik, Iceland, for the months of December, January, February, and March (https://climatedataguide.ucar.edu/climate-data/hurrellnorth-atlantic-oscillation-nao-index-station-based).
and its solutions approach /, the so-called carrying capacity for the population [21].While the values ,  are sometimes thought of as constants, this is not a necessary restriction.If  and  depend on other parameters of interest in a particular application, the underlying dynamics of the equation, in which the rate of change of  is either positive or negative depending on whether the value of  is, respectively, below or above a certain level, are unchanged.Herein, we seek a logistic model in which ,  could possibly depend on several parameters of potential relevance in this setting.

Finding a Model.
Models for ctenophore population in terms of these water masses were found using Eureqa [22].This is a scientific data mining software package that searches for mathematical patterns within data [18].We entered into Eureqa the time series data for the years 1981 through 2010 for the ctenophore population, North Atlantic Oscillation index (), and for each of the water mass indices Labrador Subarctic Slope Water (ℓ and ), Atlantic Temperate Slope Water (), and Gulf Stream North Wall ().To fill in missing data points within the collected data, we used Eureqa's interpolation option in which a missing data point is filled in by the mean of two adjacent data points [22].Eureqa will search for any model format but also allows the user to look for models of a specific type.Since we are interested in how the ctenophore population may be affected by NAO index and the water mass indices specified above, we chose to force the software to search for a logistical model of ctenophore population ().That is, we required the software to search for functions  1 and  2 such that where  1 and  2 are functions of , ℓ, , , and .That is, Initial software runs used only the raw data, but output models from Eureqa were highly sensitive to small perturbations of data.As a result, we acknowledged the existence of noise in collected data, and before running Eureqa again, we smoothed the data using Eureqa's smoothing process, which is determined by generalized cross validation among cubic bsplines [18].

Sensitivity Analysis.
Once we select an appropriate model output from Eureqa, the remaining question is, which of the parameters , ℓ, , ,  is the more influential on the ctenophore population ()?We used a standard local sensitivity analysis [23] in which we use small (local) perturbations of each parameter in (2) and examine how each perturbation affects the model output.
For ease of notation let (, ℓ, , , , ) denote the right hand side of (2).That is, let (, ℓ, , , , ) = ( 1 −  2 ⋅ ) ⋅ , where  1 =  1 (, ℓ, , , ),  2 =  2 (, ℓ, , , ), and let   ,  ℓ ,   ,   , and   represent the standard deviations of the time series data sets of , ℓ, , , and , respectively.For example, to illustrate the sensitivity of the ctenophore population to a change in , we solve the differential equation ( 2) with the known values of  modified by adding or subtracting a fraction of   .Since this is a local sensitivity analysis, we chose small varying degrees of change in  by increasing (and decreasing) the values of  by a small fraction of   .In our case, we chose fractional changes of 10% of   and 40% of   .That is, we replaced  with  + 0.1  and found a numerical solution to / = (, ℓ,  + 0.1  , , , ) to see how this solution deviates from the solution to / = (, ℓ, , , , ).We also obtained numerical solutions to / = (, ℓ, −0.1 , , , ), / = (, +0.4,, , ), and / = (, ℓ,  − 0.4, , , ).This process would be repeated for each parameter that appears in Eureqa's model output.Plots of the percent frequency of occurrence of ctenophores () generated from the solution to the differential equation of the model with the perturbed data as a function of year can be compared to the plots of the original, unadjusted model.By changing each parameter one-by-one, we are allowed to see which (if any) parameter has a stronger effect on the model output.

Results
Eureqa generated 7.7029 × 10 7 possible functions, resulting in 18 proposed formulations to the differential equation.These formulations are summarized in Table 1.
In order to select which formulation is the best available, two criteria are kept in mind.The first is that, in the logistic equation ( 2), the values for  1 and  2 should always be positive: this is necessary for the usual dynamics of the logistic population model.The subsequent criterion is to obtain the highest available correlation between the model's predictions and the measured population data.
The three formulations which meet the first criterion are indicated with ( * ) in Table 1.Since the North Atlantic Oscillation index () can take on negative values, the other solutions do not guarantee that both  1 and  2 are always positive.Among the three feasible formulations in Table 1, the following model achieved the best correlation with the (smoothed) input data, with a correlation coefficient of 0.96 and  2 best fit of 0.91: That is, Eureqa found  1 and  2 from our setup equation ( 2) such that  1 = 0.009163 () and  2 = 1/(1 +  −0.2263(−5) ).
Once the modeled formulation for the differential equation was selected from Eureqa, Mathematica [19] was used to find a numerical version of the solution () to the differential equation.We graphed the solution to (3) with a line plot of the actual collected data (Figure 1).For illustrative purposes, we also provide a graph of the solution to (3) with a line plot of a simple smoothed version (moving average with box size 3) of the original data (Figure 2).
Observe that none of the three candidates for our differential equation solution (see Table 1) used any parameter other than  or .Thus, the model output from Eureqa demonstrates that ctenophore population can be effectively modeled using only the North Atlantic Oscillation index () 0.791 0.896 0.918 0.960 0.922 0.961 0.922 0.961 and percent of Labrador Subarctic Slope Water ().Whatever effect the other parameters may have on the population is dwarfed by the importance of these two parameters.Since the model from Eureqa uses only the parameters  and , a sensitivity analysis was necessary only for those parameters.(, ℓ, , , , ) as described in Section 2.3 is merely (, , ), where (, , ) = (0.009163 () − (1/(1 +  −0.2263(−5) ))())().Thus graphs of numerical solutions to / = (,  ± 0.1  , ) are presented in Figure 3, and graphs of / = (,  ± 0.4, ) are in Figure 4. Numerical solutions to / = ( ± 0.1  , , ) (Figure 5) and / = ( ± 0.4  , , ) (Figure 6) are also presented.These illustrations suggest ctenophore population is marginally affected by small perturbations in each index.Larger perturbations show significant departures from the solution to (3).

Discussion
Our results suggest ctenophore population is most closely associated with and sensitive to changes in  (Labrador Subarctic Slope Water as formulated by NEFSC) and North Atlantic Oscillation ().In the case of relatively small changes, the population is more sensitive to Labrador Subarctic Slope Water.Larger changes in either Labrador Subarctic Slope Water or North Atlantic Oscillation correspond to dramatic population changes.This result is consistent with previous analyses of dynamics of copepod populations on the Northeast Continental Shelf [13] where the increase in Labrador Subarctic Slope Water acts to cut off warmer slope water from the Gulf of Maine and instead injects cold fresh water.
The relationship uncovered by our results raises a question for future research about the dynamics that affect ctenophore population on the shelf.Is the change in ctenophore abundance a result of enhanced water conditions (e.g., nutrients, temperature, and stratification) or are the water mass changes we associate with the occurrence of ctenophores transporting ctenophores onto the shelf and inside the range of our surveys?
Our results prescribe a 5-year lag between Labrador Subarctic Slope Water and ctenophore population.Hare and Kane [13] present a timeline from climate to biology on the order of 4 years leading to changes in a copepod population.Therefore, while a detailed understanding of the relationship between ctenophores, Labrador Subarctic Slope Water, and North Atlantic Oscillation does not exist and is beyond the scope of this study, we feel our 5-year lag is not an unreasonable result.Further study is needed to properly attribute the lag we find to ctenophore ecology.

Conclusion
This analytical technique suggests Labrador Subarctic Slope Water and the North Atlantic Oscillation have significant impact on ctenophores on the Northeast U.S. Continental Shelf Ecosystem and can potentially be used as a leading indicator of change in ctenophore population size.With a lack of mechanistic understanding of ctenophore population dynamics at this time and space scale, this technique provided a good first-pass analysis to identify which oceanic forcings are most relevant to ctenophore populations.The nature of the three-step approach seems conducive to many similar problems in biological oceanography.

2. 1 .
The Logistic Model.The changing size of a population over time is often modeled by a differential equation.A variety of standard approaches are available, including the logistic equation.Letting  = () denote the size of a population at time , this equation takes the form   = ( − ) , where  > 0,  > 0,

Figure 1 :
Figure 1: Ctenophore population () as a function of year.The black line plot is raw data from the Link and Ford [6] ctenophore index from 1986 to 2008.The red curve describes the ctenophore population as numerical solution (via Mathematica) to the differential equation provided by Eureqa.

Figure 2 :Figure 3 :
Figure 2: Ctenophore population () as a function of year.The black line plot is a three-year moving average of the raw data from the Link and Ford [6] ctenophore index from 1986 to 2008.The red curve describes the ctenophore population as a numerical solution (via Mathematica) to the differential equation provided by Eureqa.

Figure 4 :Figure 5 :
Figure 4: Ctenophore population () as a function of year.The black line plot is raw data from the Link and Ford [6] ctenophore index from 1986 to 2008.The red curve describes the ctenophore population as a numerical solution to the differential equation provided by Eureqa.The blue line is a solution to the differential equation with  replaced by  + 0.4  .The green line is a solution to the differential equation with  replaced by  − 0.4  .

Figure 6 :
Figure 6: Ctenophore population () as a function of year.The black line plot is raw data from the Link and Ford [6] ctenophore index from 1986 to 2008.The red curve describes the ctenophore population as a numerical solution to the differential equation provided by Eureqa.The blue line is a solution to the differential equation with  replaced by  + 0.4  .The green line is a solution to the differential equation with  replaced by  − 0.4  .

Table 1 :
The proposed solutions to the differential equation generated by Eureqa.