To capture the first-order characteristics of the distribution of eDNA shed by mesopelagic species, we used a mechanistic model with one spatial dimension (the vertical dimension). It treats migrating organisms as a continuous point source of eDNA and simulates temporal evolution of the eDNA concentration vertical profiles in a one-dimensional ocean with climatological conditions. Horizontal variation of the distribution of eDNA shed by mesopelagic species is not considered in this study. The point source in the model represents a number of individuals within a single species. The simulation includes relevant eDNA transport processes using a vertical Price-Weller-Pinkel (PWP) dynamical module42. We use this model to first explore the sensitivity of eDNA concentration profiles to a variety of biological and physical parameters and to assess vertical and temporal variability of the eDNA concentration. Table 1 outlines the parameters, separated by type, and ranges of values for each parameter used in the sensitivity analysis (see “Sensitivity analysis” section below) and in a case study of how the model can be used to determine the percent of individuals that migrate (see “Application to ecological questions” section below). We then discuss how the mechanistic model can be used in conjunction with field sampling to test ecological hypotheses, highlighting how analysis of eDNA concentration profiles can provide information regarding where and when organisms are located in the water column. Finally, we comment on implications of this work and how the model can be used to both design and interpret observations.
Model overview
Governing equations
To maintain simplicity and also provide a more realistic scenario where both settling and neutrally buoyant eDNA particles28 are considered, we assume that the eDNA material at any given time in the system consist of eDNA particles of two size classes. Here, large eDNA particles are on the order of 10 s of µm to 1 mm and represent materials such as fecal pellets, chunks of tissue, or gametes that would be subject to both settling and physical breakdown over time. Small eDNA particles are on the order of 0.1 µm to 10 s of µm and would be any eDNA shed from an organism that has a density such that they are near-neutrally buoyant including extracellular eDNA or forms such as sperm, urine, blood, or single cells. This size class of small particles is so small that we can neglect their settling. In practice, these two size fractions of eDNA would be captured using the common method of filtering water through a < 1 µm pore size filter43. Both small and large particles in the model are subject to the same eDNA decay rate constant and large particles break down over time into small particles.
The equation describing the change of eDNA concentration in the form of large particles is:
$$begin{aligned} frac{partial {C_{LP}}}{partial {t}} = – w_vfrac{partial {C_{LP}}}{partial {z}} + frac{partial {}}{partial {z}}left( kappa _zfrac{partial {C_{LP}}}{partial {z}}right) – kC_{LP} – delta C_{LP} + w_sfrac{partial {C_{LP}}}{partial {z}} + hat{S_{LP}} end{aligned}$$
(1)
where (C_{LP}) is the concentration of large eDNA particles, (w_v) is the vertical velocity [L/T] with a positive value pointing upward, (kappa _z) is the vertical diffusivity [L(^{2})/T], k is the first order decay rate constant [1/T], (delta) is the breakdown rate of particles [1/T], (w_s) is the settling rate of eDNA [L/T], and (hat{S_{LP}}) is the shedding rate of new large eDNA particles.
Similarly, the equation describing the temporal evolution of the concentration of small particles of eDNA, (C_{SP}), is:
$$begin{aligned} frac{partial {C_{SP}}}{partial {t}} = – w_vfrac{partial {C_{SP}}}{partial {z}} + frac{partial {}}{partial {z}}left( kappa _zfrac{partial {C_{SP}}}{partial {z}}right) – kC_{SP} + delta C_{LP} + hat{S_{SP}} end{aligned}$$
(2)
where the parameters are the same as Eq. (1). Note that there is no settling term in Eq. (2) as small eDNA particles are assumed to be neutrally buoyant and that the (delta) term in Eq. (2) has the opposite sign as that in Eq. (1), reflecting eDNA volume conservation during the breaking-down of the large particles into small particles.
Physical model
The 1-D vertical model is implemented in MATLAB. The physical module is based on a previously published 1-D model of physical processes in the mixed layer42 and is modified to include the relevant organism and eDNA parameters. The model spans from the surface down to 1500 m deep with a vertical resolution of 0.5 m and has a time step of 10 s. Each simulation runs for a total of 90 days, representing one of the four seasons: summer (July, August, September – JAS), fall (October, November, December – OND), winter (January, February, March – JFM), and spring (April, May, June – AMJ). The model has a zero-gradient open boundary condition at the bottom. The model is forced on the surface by climatological, 8-hourly meteorological conditions from the NCEP/NCAR Reanalysis at a location in the Northeast Atlantic Ocean (37.047(^circ) N, -71.25(^circ) W). Inputs include air temperature, short wave radiation, long wave radiation, precipitation, air pressure, air humidity, and winds. The wind stress and heat fluxes are calculated using the bulk formulae44.
The initial vertical distribution of the temperature and salinity in the model are seasonal climatological profiles obtained from the NCEI World Ocean Atlas 2018 at a site in the Northwest Atlantic slope sea (39.125(^circ) N, −70.875(^circ) W; north of the Gulf Stream). The model also simulates the influences of vertical mixing and advection on the temperature, salinity and eDNA particles. The mixing influence is incorporated by prescribing synthetic vertical profiles of vertical diffusivity that combines observed seasonal mean mixed layer depth45 and simulated seasonal mean vertical diffusivity profile from an operational model46, both at a site in the Northwest Atlantic slope sea (see Supplemental Text S1 and Supplemental Fig. S1). Vertical advection profiles are set to increase linearly from 0 ms(^{-1}) at the surface to the maximum value of ((10^{-4}) ms(^{-1})) at 200 m and then decrease linearly back to 0 ms(^{-1}) at 400 m (Supplemental Fig. S2). This prescribed advection profile does not change seasonally and represents enhanced vertical motions associated with sub-mesoscale processes in the upper water column (e.g.,47). Note that the result of this study is not sensitive to the prescribed vertical profiles of diffusivity or vertical velocity (see below).
The physical parameters that will be adjusted in this study include: the initial temperature and salinity profiles, the mixed layer depth, the vertical diffusivity profile, and the vertical velocity profile (see Ocean parameters in Table 1).
Organism movement
This study focuses on eDNA shed by vertically migrating organisms living in the mesopelagic ocean. For simplicity, we assume the organisms reside at the constant depth of 500 m during the day and the constant depth of 50 m during the night. These depths fall within the range of where migrating organisms are found at these times41,48. As the objective of this study is to develop a qualitative understanding of the vertical distribution of the eDNA concentration, the exact depth a particular mesopelagic species resides at in the daytime is not crucial. Organisms begin migrating to the surface two hours before sunset, and the upward migration ends one hour after sunset. The time of sunset and sunrise are determined seasonally using NOAA’s ESL solar calculator for the year 2019 at 42.35(^circ) N, −71.05(^circ) W. Both upward and downward migrations are assumed to be in constant speed, and the migration times do not change within each season (Fig. 2).
Representative organism migration curves for each season. Migration times change with sunrise and sunset for each season. This illustration assumes that 50% of individuals migrate. The shading shows the percentage of individuals with black indicating 100%, i.e., all organisms, and grey indicating 50%. These numbers vary among the simulations.
The migration parameters that can be adjusted in the model include: daytime residing depth, nighttime residing depth, the start and end times of migrations (and thus duration), the percent of individuals that migrate, and the layer thickness (i.e., “width” of school). (See Organism parameters in Table 1). In this study, we focus on the start and end times of migrations (and thus duration) and the percent of individuals that migrate.
eDNA parameters
We assume organisms are continuously shedding eDNA, and the shedding rate is fixed in time in each simulation. The eDNA then remains in the water column, where it is subject to transport and decay as described by Eqs. (1) and (2). Although several studies have characterized eDNA shedding rates of different marine organisms29,49,50, shedding rates are highly variable across studies and have high error associated with them, likely due to high temporal variability29. Here, we use a constant shedding rate of 1 mass unit per time step (10 seconds). Note that this study focuses on the relative vertical distribution of the eDNA concentration and its temporal variability, rather than its absolute value. The particle size distribution of eDNA (and the breakdown rate of large to small particles) and the settling rate of eDNA are largely unknown and expected to be time-varying28. Here we use a breakdown rate of fecal pellets to apply to large eDNA particles breaking down into small eDNA particles51. The most well characterized parameter of eDNA is decay rate, which several studies have characterized as a function of water temperature29,49. The upper and lower limits of the decay rate have also been well established29,31. Finally, there are currently no estimates of eDNA settling rates. We use values of marine snow settling rates52 to apply to large eDNA particles that are subject to settling in the model.
The eDNA-related parameters that are adjusted in this mechanistic model include: settling rate of large particles; ratio of large particles to small particles that are shed by an organism; eDNA breakdown rate (large particles to small particles); the eDNA decay rate. (See eDNA parameters in Table 1). Note that shedding rates in the model do not change over the course of each simulation.
Sensitivity analysis
A series of 90-day simulations were carried out with altered values of the aforementioned parameters to examine the sensitivity of the eDNA vertical distribution to the parameters. Table 1 provides a list of the parameters that were adjusted in the sensitivity analysis. Note that the diffusivity profile (both the shape of the profile and the mixed layer depth), the temperature profile (which impacts the eDNA decay rate), and the daytime length (which impacted migration times) all change with the season. Other sensitivity parameters include vertical advection, settling rate of large particles, the ratio of large to small eDNA particles shed by organisms, and the percent of individuals that migrate. A total of 972 simulations were conducted with the following sets of parameters: 4 mixing profiles (one for each season), 3 vertical advection profiles (upwelling, downwelling, none), 3 settling rates of large eDNA particles, 3 decay rate constant scenarios (two constant values representing high and low values in the literature31, one temperature-dependent rate29 for each season, Supplemental Fig. S3), 3 scenarios for the percent of individuals that migrate, and 3 ratios of large to small eDNA particles shed by the organisms.
In order to assess the impact of the parameters on the eDNA profiles, several metrics were defined. First, the water column was divided into three sections: surface layer (0–100 m), mid-depth (100–450 m), and deep water (450–550 m), based on our representative daytime and nighttime depths (500 m and 50 m) used for the simulations. For each simulation, the mean and maximum eDNA concentration in each depth bin was recorded. Also, the cumulative eDNA concentration in each depth bin was normalized to calculate the proportion of eDNA in each depth bin at any given time during the simulation. For each simulation, the mean and standard deviation of the proportion of eDNA in each depth bin were calculated over the course of the 90 day simulation.
A first-order comparison of the vertical length scales of eDNA transport by different processes were conducted. Here, we use a time scale of eDNA decay, T(_{90}), i.e., the time it takes for 90% of the released eDNA to decay, to determine the vertical length scale of each process. In particular, we would like to estimate, the vertical distances eDNA is transported by advection, mixing, and settling before 90% of the released eDNA has decayed, L(_{mix}), L(_{advect}) and L(_{settling}), respectively. These quantities are defined mathematically as,
$$begin{aligned} T_{90} = frac{-ln(0.1)}{k_{avg}} end{aligned}$$
(3)
$$begin{aligned} L_{mix} = sqrt{kappa _v T_{90}} end{aligned}$$
(4)
$$begin{aligned} L_{advect} = w_{vm} T_{90} end{aligned}$$
(5)
$$begin{aligned} L_{settle} = w_s T_{90} end{aligned}$$
(6)
where (k_{avg}) is the average decay rate constant [T(^{-1})] for the whole water column over the 90 day simulation, (kappa _v) is the maximum vertical diffusivity coefficient [L(^{2})T(^{-1})] , (w_{vm}) is the maximum vertical velocity [LT(^{-1})] , and (w_s) is the settling rate [LT(^{-1})]. Note that because k is temperature dependent, the value of the decay rate constant will change both vertically (e.g., deeper water will be colder and have a lower decay rate constant) and temporally (e.g., at a given depth, the water temperature will be warmer in the middle of the day and thus will have a higher decay rate constant).
Application to ecological questions
After obtaining a first-order understanding of the eDNA concentration profile and the impact of the parameters on the eDNA distribution, we ran another set of model simulations to examine if vertical and temporal eDNA concentration variability can shed light on an ecological question regarding what percentage of individuals within the same species migrate on a daily basis. To do this, for each season, all parameters were held constant except for the percent of individuals that migrate, P(_{m}). In each season, the value of P(_{m}) varies from 0 to 100% with an increment of 10%. A total of 44 such simulations (4 seasons x 11 percent migrate scenarios) were carried out. We then compared mean and maximum eDNA concentrations in the surface layer (0–100 m) to those in the deep water (450 to 550 m), and examined the relationship between the surface-to-deep ratios and P(_{m}).
Source: Ecology - nature.com
