An affordable and automated imaging approach to acquire highly resolved individual data—an example of copepod growth in response to multiple stressors

Individual trait variation is essential for populations to cope with multiple stressors and continuously changing environments. The immense number of possible stressor combinations and the influence of phenotypic variation makes experimental testing for effects on organisms challenging. The acquisition of such data requires many replicates and is notoriously laborious. It is further complicated when responses occur over short time periods. To overcome such challenges, we developed an automated imaging platform to acquire temporally highly resolved individual data. We tested this platform by exposing copepods to a combination of a biotic stressor (predator cues) and a toxicant (copper) and measured the growth response of individual copepods. We tested the automatically acquired data against published manually acquired data with much lower temporal resolution. We find the same general potentiating effects of predator cues on the adverse effects of copper, and the influence of an individual’s clutch identity on its ability to resist stress, between the data obtained from low and high temporal resolution. However, when using the high temporal resolution, we also uncovered effects of clutch ID on the timing and duration of stage transitions, which highlights the importance of considering phenotypic variation in ecotoxicological testing. Phenotypic variation is usually not acknowledged in ecotoxicological testing. Our approach is scalable, affordable, and adjustable to accommodate both aquatic and terrestrial organisms, and a wide range of visually detectable endpoints. We discuss future extensions that would further widen its applicability.


INTRODUCTION
The life of any organism is a continuous struggle with different stressors, be it from other organisms or the physical environment. Since the last century organisms are also exposed to novel artificial substances of anthropogenic origin, such as chemical toxicants, as well as rapid changes in the environment due to human activities. In nature, stressors never act independently of one another. In addition to anthropogenic stressors, organisms must cope with natural biotic stress from parasites, competition for resources, and predation risk. While direct consumption is detrimental for prey, non-consumptive effects are complex (Heuschele et al., 2014), and include behavioural, morphological, and physiological changes (Verity & Smetacek, 1996). Most studies on mixture toxicity and multiple stressors focus either on interactive effects of two toxicants, or one pollutant in interaction with changes in the physical environment (Gunderson, Armstrong & Stillman, 2016). The effects of combined stressor exposures range from synergistic to antagonistic when compared to single stressor exposure (Rose, Warne & Lim, 2001;Fischer, Roffler & Eggen, 2012;Holmstrup et al., 2010), and the scales and timing of response differ widely, rendering it challenging to predict the outcome of additional stressors (e.g. Segner, Schmitt-Jansen & Sabater, 2014).
However, such interactions between biotic stress and toxicants might be the rule rather than the exception and complicated indirect effects on predator-prey relationships seem to be common in aquatic communities (Rohr & Crumrine, 2005;Langer-Jaesrich et al., 2010;Trekels, Van de Meutter & Stoks, 2013).
Phenotypic differences within a population are another source of variation that complicates predicting multiple stressor effects. Recognizing this trait variation has profound and practical consequences for ecotoxicology, but also human medicine (Evans & Relling, 2004). Nevertheless, effects of stressors are often tested separately using laboratory populations of limited genetic diversity (OECD, 2012;Macken, Lillicrap & Langford, 2015). The use of unique strains and clones inherently misses the ecotoxicological target 'to predict effects in real populations' (Lam & Gray, 2001).
The tremendous number of possible stressor combinations and the influence of phenotypic variation on the biota response poses a grand challenge for ecotoxicology, which is seemingly impossible to deal with, and yet unavoidable. Gathering individual data is more laborious compared to gathering pooled or group data, especially when it comes to following the development of individuals. The sensitivity of individuals is likely highest at specific life stages and stage transitions. Therefore, our focus should be to identify these states rather than only addressing the larger time scale responses related to life span, fecundity, and adult survival. While the development rate of an animal can be estimated from daily observations, to actually measure growth parameters e.g. length or examine rapidly occurring ontogenetic events e.g. timing and duration of stage transitions, requires the frequent observer presence and potentially the repeated handling of the experimental animals. Even the mere observer presence may involuntarily alter the behaviour and development of individuals (Mallet et al., 1987;Baker & McGuffin, 2007). Hence using traditional methods to gather individual data, such as manually extracting, measuring, and placing them back into their holding container, can lead to biased results from greater effects of handling than treatment.
This challenge calls for the application of efficient and automated testing platforms to execute well-designed and manageable experiments that compare stress responses of organisms. In recent years, automated experimental systems have been developed to monitor water quality by quantifying motility and other endpoints in indicator species ranging from single-celled Euglenas (Tahedl & Häder, 2001;Lee, Zheng & Yang, 2012), over Daphnia (Häder & Erzinger, 2017) to fish (Cunha et al., 2008). There is a range of automated systems which find use as in real-time monitoring of behavioural responses of aquatic organisms (reviewed in Bae & Park, 2014).
Among automated systems microfluidic Lab-on-a-chip systems (Campana & Wlodkowic, 2018) allow for the precise dosing in biochemical assays treatments. Due to their small size and cost efficiency they could be used in high-throughput screening of new chemicals. Measured endpoints range from bioluminescence production (Zhao & Dong, 2013), to viability (Gammoudi et al., 2014), and motility (Huang, Campana & Wlodkowic, 2017). One drawback is however that small volumes can limit the size of the testable animals (Campana & Wlodkowic, 2018) and might hinder a 'natural' response to the treatment, especially when measured throughout the complete development period of an organism. The aim of this study was thus to develop and test an affordable automated imaging system that allows for the continuous observation of a large number of separately kept individuals. We validated our approach by comparing our detailed data and results to the ones from Lode et al. (2018) which are based on temporally less resolved data from the same experiment.
We used copepods as model organisms, as they are key players in marine pelagic food webs and the most abundant metazoans on the planet (Humes, 1994;Naganuma, 1996). In recent years, copepods are also increasingly used as models in ecotoxicology (Macken, Lillicrap & Langford, 2015;Raisuddin et al., 2007).
The experiment followed the development of individual copepods under the influence of a biocide (copper) and a natural stressor in the form of chemical cues of a fish predator. In our case these include both kairomones, chemical cues emitted by the fish, which benefits the receiver and potentially harms the emitter, and Schreckstoff, cues from eaten copepods that can warn the other individuals. For simplicity we refer to these as predator cues in the remainder of the manuscript. We used a sublethal concentration of copper. We thus expected small but accumulative effects that would likely affect specific developmental stages. We anticipated that the highly resolved data would allow us to determine the most affected development stages, and also uncover subtle changes in growth trajectories and in the duration of stage transitions.

Automated platform for image acquisition
We used an automatized imaging setup to follow individual growth at an hourly resolution, with a self-made experimental system that is capable of filming single culture plate wells repeatedly over the course of a copepod's development time. We used a DIY plotter kit (Makeblock Co., Ltd., Shenzhen, China: XY-plotter Robot Kit V2, see Fig. 1) as a basis for the system. On the movable platform, we mounted an upward facing infrared-capable camera with an image resolution of 2,592 Â 1,944 pixels (Raspberry Pi NOIR with C mount and a six mm adjustable-focus lens). The images were saved as jpgs with a moderate compression of 85 to reduce file size while maintaining details, with 100 representing the maximum possible quality. Above the camera system, we installed a platform made of transparent Plexiglas. On this, we placed four 24-well plates containing the animals. Two stepper motors move the camera position. Instructions are taken from a python script on a microcomputer (Raspberry Pi), which also controls the camera. We programmed the system to sequentially take one image (2,592 Â 1,944 pixels) of each well approximately every hour, for a period of 13 days. Two LED lamps (Camlink CL-Studio10) provided constant illumination to the setup from opposing sides, providing dark field imaging conditions. After each well plate the system reset its positioning system using two contact switches at each axis. This prevents a continuous systematic error in cases when the system gets misaligned. We later discovered that such errors occurred sporadically to be due to non-optimal baud rate settings of the serial port. The source code for the imaging system is included in Supplemental Material S1. The material costs of the system summed up to less than 500 EUROs.
The build of the setup is for the largest part easy as the plotter kit is targeted at juveniles. The execution of the script and adding changes to the number of plates, timing of recordings can be done by any person with basic python knowledge. Userfriendliness could however easily be improved by adding a simple graphical user interface to the script.

Experiment and animals
We tested the platform with an experiment on combined effects of predator cues and copper exposure on copepod age and size at maturity. We used the harpacticoid copepod Tigriopus brevicornis as the model organism. The laboratory stock cultures originated from a splash water pool in Drøbak (Norway) and one from Tjärnö (Sweden). We kept the stock cultures at 30 psu, 18 C and a 12/12 light dark cycle for more than six months before the experiment. Stock cultures were fed ad libitum three times a week an equal mix of Rhodomonas salina, Isocrysis galbana, and Dunaniella tertiolecta.
Prior to the experiment, we picked single females with egg sacs from the stock culture. We placed them individually into 24 well plates. Every 30 min we checked manually if nauplii had hatched. If nauplii were detected, we registered the time, assigned a clutch ID, and removed the female. The nauplii were then placed individually in wells containing 2.5 ml of water with the respective treatment.
We exposed copepods to one of four treatments: predation risk (predator cues), copper (20 mg L -1 ), combined predation risk and copper (20 mg L -1 ), or control (pure seawater). All seawater was taken from the outer Oslofjord and filtered at 1.2 mm prior to use. We prepared seawater with predator cues by incubating three-spined stickleback (Gasterosteus aculeatus) for 48 h in filtered seawater at 18 C, 30 psu, two fish l -1 . The fish were fed with T. brevicornis first at initiation, and once more after 24 h. Following incubation and removal of fish, the water was filtered (GF/C, 1.2 mm) and frozen (-18 C). Water without predator cues was prepared similarly but without addition of fish and copepods. These frozen bottles were thawed daily to prepare the exposure solutions. For Cu and the combined treatment, we then added Cu through a 2-step dilution process of a 0.1 M CuSO 4 stock solution. Instead of Cu, we similarly added distilled water for the control and the predation risk treatment.
We replenished 72% of the exposure solution daily for each individual. We inspected the individuals daily using a binocular microscope and recorded survival, took a photo for subsequent length measurements, and most importantly assessed the development stage. We used the numbers of exuviae to determine the stage of the copepods. For a more detailed description of treatment preparation and the general procedure see Lode et al. (2018).

Image acquisition
We incubated 72 individual copepods from nine clutches of different mothers. One was lost during the setup of the experiment, and two could not be followed due to a misalignment of the robot. From the 69 other ones, we managed to extract on average 168.9 ± 31.1 measurements for each individual during development. The upper number of images and thus the maximum amount of measurements was 438.

Length measurements
In total, we acquired 11,657 images that were suitable for measurements. Some of the other pictures were misaligned due to occasional glitches in the computer drives platform, or because the system took them during the daily water change, when the wells were removed from the platform (Lode et al., 2018) These images had to be removed from the analysis. We imported the images into the imaging software Fiji (Schindelin et al., 2012). The high magnification reduced the depth of field and the copepods could only be measured when they were in close proximity of the bottom (Fig. 1C). If the animal was distinctly visible and in focus, we manually measured its length. Due to the benthic lifestyle of the nauplii compared to the demersal behaviour of adults, we were able to obtain more measurements from younger stages. To compare the length measurement from machine images and the one obtained from microscope image by Lode et al. (2018), we calculated a daily average of the machine length measurements.

Moult from nauplii to copepodite
We determined the time of moult to the first copepodite stage (i.e. moment) by manually screening the pictures for the first occurrence of the copepodite stage. Earlier moults from nauplii to nauplii or moults from copepodite to copepodite were not as obvious and could not be easily determined directly from the images. Sometimes a water change, a temporary misalignment of the camera, or the temporary removal during the manual screening led to unusable pictures. If such unusable pictures preceded the first appearance of a copepodite stage, we were not able to accurately determine the time point. To be able to account for the varying degrees of uncertainty we noted down the number of 'uncertain' pictures, and included this information in the statistical analysis.
From the images we saw that most nauplii remained motionless for several frames during their transition to the first copepodite stage. To be able to determine whether the treatments influenced the moult duration to copepodite, we recorded the number of 'motion free' pictures before the first appearance of the copepodite. We included the first picture in which they were detected in the same position.

Statistical analysis
All statistical analyses and data manipulation were conducted using the statistical software R (Version 3.5.0) (R Core Team, 2018). We used all acquired length measurements for the analyses. We converted the measurement time to the actual age of the individual using the time of birth (±30 min uncertainty). To meet normality assumptions, copepodite transition time and length data were log-transformed before analysis, and data of the duration of the nauplii to copepodite transition was square root transformed. If not otherwise stated, we used an information theoretic approach to select the best model for the linear models based on the corrected Akaike information criterion (AICc) for all measured response variables. If there were several models within a DAICc < 2 of the best one, we averaged model estimates. We used the package MuMIn (Barto n, 2018) for model selection and averaging.
To test whether length measurements from the images where comparable to the manual measurements from Lode et al. (2018), we used a linear model with machine as dependent, microscope measurements as fixed factor, and individual ID as random factor to control for repeated measurements of the same individual. We also included treatment as fixed factor in the initial model, to test for differences in measurement 'accuracy' between treatments.
We tested the influence of treatment and clutch ID on transition timing from the last nauplii to the first copepodite stage with a two-way ANOVA. We allowed for an interaction between both independent factors in the initial model. We controlled for the uncertainty in transition time, by including the number of 'uncertain frames' as weights in the analysis.
To test whether males and females were affected differently by the treatments, we censored the data to include only matured individuals. We then examined the influence of treatment and gender in a separate model with these two factors as fixed factors, and copepodite transition time as the dependent variable. As before we allowed for all interactions between the independent factors and the number of uncertain frames as the weighting factor.
We tested for differences between clutch IDs and treatment in the moult duration of the nauplii to copepodite transition using a linear model with duration as dependent variable, and clutch ID and treatment as fixed factors. We allowed for interactions between the fixed factors in the initial model, and always included the number of uncertain frames as the weighting factor.
We tested for the influence of copper, predator cue, and clutch ID on length development using general additive models (GAM) with thin plate regression splines, using the mgcv-package (Version 1.8-23) for general additive modelling (Wood, 2003(Wood, , 2011. Our models included individual ID as random factor, allowing for a random smoothing over time for each individual. In this case, to choose the best model, we started out with a model allowing for full interactions between the fixed factors and took advantage of the inbuilt model selection tool of the mgcv-package where the smoothing parameter estimation allows for model terms to approach zero. This procedure results in a final best fitting model. The final model parameters and smoothing functions were then evaluated using the function gam.check() and based on k`, estimated degrees of freedom (edf) and p-values. We further visually inspected the distribution of the residuals, quantile-quantile plots, and residual vs linear predictions.
Except for the moult from nauplii to copepodite, it was challenging to determine the exact time point from one nauplii stage to the next by manual (visual) inspection of the acquired images. Therefore, we tested for differences in transition timing between nauplii stages by analyzing the predicted individual growth increments over time. For this, we first ran a GAM that included random smoothers for each individual only. We then calculated isochronal body length predictions on a 72-minute resolution, similar to the one of the raw data. To not use the model beyond the data range, individual level length predictions were limited to the respective time the individuals spent in the experiment.
From these, we derived the individual growth increments, which showed distinct peaks representing the stage transition phases.
We then used a new GAM with growth increments as the dependent factor to test for interactive effects of Clutch ID, copper, and predator cues on the transition dynamics. We also included negative predictions of growth to keep the normally distributed nature of the data and to emphasize the growth spurts during moults, although they are biologically impossible in this species.
To test whether copper, predator cues, or clutch ID influenced the final size of the copepods at the end of the experiment we restricted the data to measurements taken during the last day of the experiment and averaged them for each individual. We then analyzed this relationship using a linear model with copper, predator cues, and clutch ID as independent factors, and the log-transformed (averaged) body length as the dependent factor. In this case the endpoint size encompasses both potential effects of developmental delay and on size. To test whether treatment influenced only the final size of males and females differently we further restricted the averaged length data to include only lengths of fully matured individuals. We then used a similar model structure to the previous one but added gender as a fixed factor in the initial model.

RESULTS
The daily manually measured length data by Lode et al. (2018) and data acquired using the robot  Figure 3 Treatment effects on the timing of N6 to C1 metamorphosis. The moment in time when nauplii N6 turned copepodites C1 depending on the copper treatment, with shapes indicating the sex as identified by the end of the experiment. Sex "Copepodites" are individuals that did not reach maturity by the end. The red dot and line indicate the mean and the SE of the raw data. although automated length estimates were generally larger than the manual measurements, especially during the early copepodite stages (Fig. 2). Treatments did not influence the relationship between machine and manual measurements. Copper delayed the time of the moult from nauplii N6 to copepodite C1, while the transition was independent of an individual's clutch ID ( Fig. 3; Table 1). Model estimates for the N6 to C1 transition timing from nauplii to copepodite in males and females were similar, and the effect was driven by individuals that did not reach maturity by the end of the experiment (Fig. 3).
In contrast to the time of moult, the moult duration from nauplii N6 to copepodite C1 was mainly influenced by the individual's clutch ID in interaction with copper ( Fig. 4; Table 1). Individuals of six clutches showed a reduced moult duration, while in the other clutches the duration was prolonged compared the control (Fig. 4). The duration of the moult ranges from the time between two subsequent recordings (∼72 min) to more than 400 min.
The development of individuals, measured as individual growth, is influenced by a complex interactive effect of copper, predator cues, and clutch ID ( Fig. 5; Fig. S1; Table 2). The treatments left two clutches unaffected (Fig. S1), while individuals in all other families responded with delayed growth. When we focus on the treatment effects, copper alone delayed development while predator cues did not have an impact. The combination of both led to a stronger delay in the late copepodite stages (Fig. 5).
The conversion of length measurements to growth increments showed distinct peaks which revealed the moults of the copepods, with variation in growth increments explained by an interaction between Clutch ID, copper, and predator cues ( Fig. 6; Table 3). Overall, the most significant delay in development occurs in individuals exposed to copper and predator cues combined when they metamorphose from nauplii to first copepodite (Fig. 6), which confirms the results of the manually screened timing of this major moult ( Fig. 3; Lode et al., 2018). However already during the third transition the averaged peak height is reduced in the combined stressor treatment, which indicates a larger variability between exposed individuals. With an increased age of the individuals at the specific developmental stages, the moult cycles became less synchronized and measurement error became larger, and an overall trend between treatments was harder to detect, but also detecting individual peaks or moults becomes harder (Fig. 6). Therefore, we refrained from analyzing intermoult durations  Figure 4 Copper and clutch ID effects on moult duration. Boxplot of the differences in N6 to C1 moult duration depending on clutch ID of the animals and the presence and absence of copper, the box shows median, quantiles, and the 1.5-time interquartile range is indicated by vertical lines.
Full-size  DOI: 10.7717/peerj.6776/ fig-4 based on detected growth maxima. However, the peak heights of the averaged growth increments are highest in the control and predator cues treatment compared to the other two treatments, which means that there was less variation in transition timing in control and predator cues treatment. Or in other words, the effects of copper seem to be strongly affected by individual and clutch variation. The final body size at the end of the experiment depends on additive effects of treatment and Clutch ID. Both competing models (DAICc < 2) included additive effects by copper and clutch ID, while only one had a negative effect of predator cues (Table 4). Some clutches were unaffected by treatment while for most others the exposure by both treatments led to reduced length (Fig. S1). However, the 23 unmatured copepodites at the end of the experiment in the copper and combined treatment biased these results. Therefore, we restricted the analyses to include mature individuals only (n ¼ 49). The analysis showed several competing models (Table 4), which revealed that sex is the most important factor and to a lesser degree Clutch ID and copper. In general, males were slightly larger than females (Table 4).

DISCUSSION
The intensity of adverse responses to toxicants in the natural environment is challenging to predict due to the almost infinite number of possible interactions with biotic and abiotic factors (Segner, Schmitt-Jansen & Sabater, 2014), and begs for efficient methods to handle many replicates. In this study, we used an automated imaging approach to measure the combined impact of a toxicant (copper) and a biotic stressor (predator cues) on copepod development. We validated and evaluated the added benefits of our approach with the findings obtained using traditional manual methods at a much lower time resolution (Lode et al., 2018). We find the same complex interactions between the copper and predator cues treatment, and individual's clutch ID, in determining the growth trajectory of an individual facing multiple stressors. Compared to the daily measurements of Lode et al. (2018) our highly resolved data allowed us to zoom in on individual moult events, which is not possible for large sample sizes using traditional methods. Especially during the naupliar stages, we detected clear peaks of moult events. These revealed that the treatment effects first affected the N3 transition and got more pronounced from then onwards. The biggest effect is visible during the naupliar to copepodite (N6-C1) metamorphosis. The significant differences between individuals and different clutches in their response to the treatments led to a wider distribution of the moult timings (Fig. 6). The strong influence of an individual's clutch ID on the major intermoult duration also suggests a genetic role in the resistance to multiple stressors, which is a major challenge in ecotoxicology (Evenden & Depledge, 1997;Wirgin & Waldman, 2004).
The concurrent results show the potential of our semi-automated system to tackle large sample sizes and detect small developmental differences in individual organisms, while still reducing the workload and the handling of the animals. Our setup can thus ease the collection of individual trait data and be used to answer questions in both toxicology and ecology. A focus on trait-based responses is especially helpful in studying the responses to multiple stressors. While we used it for small aquatic invertebrates, the imaging system is customizable and adjustable to accommodate different container-and species sizes. It thus is in line with the successful use of automatic monitoring systems in ecotoxicology like for example the Multispecies Freshwater BiomonitorTM (Gerhardt & Schmidt, 2002), LeDaphNet (Lagergren et al., 2017), and DaphniaTox (Häder & Erzinger, 2017) and lab-on-a-chip systems (Zhu et al., 2015;Campana & Wlodkowic, 2018).
In its current state, the system can reliably capture the size and movement of animals with a primarily benthic lifestyle. Examples include surface cruising animals such as benthic copepods, snails, trematodes, and nematodes. Using a bottom mounted camera works best when individuals are close to the bottom of the holding container, or in our case the well plate. Potential research questions include testing the influence of abiotic and biotic factors on the settlement of planktonic larvae of benthic animals, growth development of invertebrates at different nutrient concentrations, egg hatching times, and in this context also the factors which drive the emergence of resting eggs. Our system can even quantify the behavioural variability in populations from recorded movie sequences. As individual variability is the foundation and currency of personality research, i.e. the study focusing on repeatable and correlated behaviours, the imaging system can easily be used to capture consistent differences between individuals and in consequence 'personalities' of animals. In recent years it has become clear that such characteristics that are not only present in 'higher' organisms but also in invertebrates (Kralj-Fišer & Schuett, 2014;Sih, Bell & Johnson, 2004).
The depth of field of the camera is one of the apparent drawbacks of the current system. With increasing age, the copepods became increasingly active and explored the whole water column. Thus, the number of images from which we could reliably determine body length decreased with age. This problem could be solved by using a smaller aperture or by adding a servomotor to the camera. Moving the plane of focus while recording a movie, would increase the chances of acquiring a sharp and complete image of the animal.
While the recording of the animals is automatized, our approach currently still relies on the manual screening and measuring of the images. This step is necessary due to the imperfect image quality (variable light conditions, occasional blurriness). An even background illumination using electroluminescent sheets or diffuse LED light sources does however increase the image quality to a point where it would be possible to implement image analyses based on neural networks. For example, using the tensorflow library (Abadi et al., 2016) we could then automatically classify stage data, measure size data, and other traits (gut content). The type of images could also be analyzed using crowd-based annotation services such as Quanti.us (Hughes et al., 2018).

CONCLUSIONS
Our results illustrate the need to study the interactive effects of natural and anthropogenic stressors, and they underscore the necessity to consider the phenotypic and genetic variation in stress response if we want to use ecotoxicological studies to predict the consequences of toxicants for natural populations. Our system takes the idea of autosamplers, lab-on-a chip, and other high-throughput ideas, and applies it to questions related to the development and potentially the behaviour of small invertebrates. It uncovered differences in moult duration and the timing of copepod metamorphosis which would be difficult to detect using manual approaches. Given that it is easy to build, affordable, and runs with open source imaging and analysis software, it can be scaled to accommodate for high-throughput testing of multiple treatment combinations and gradients. When data are gained at both individual and population levels, they can be combined conceptually in adverse outcome pathways and increase the value of risk assessment in ecotoxicology (Kramer et al., 2011).