Observation of animal movements on small spatial scales provides a means to understand how large-scale species distributions are established from individual behavioral decisions. Small-scale vertical movements of 14 Summer Flounder Paralichthys dentatus residing in Chesapeake Bay were observed by using depth data collected with archival tags. A generalized linear mixed model was employed to examine the relationship between these vertical movements and environmental covariates such as tidal state, time of day, lunar phase, and temperature. Vertical movements increased with warming water temperatures, and this pattern was most apparent at night and during rising and falling tides. Fish generally exhibited greater vertical movements at night, but the difference between vertical movements in the day and those at night decreased as fish increased in size. Results from this study fill a void in understanding the small-scale movements of Summer Flounder and could be incorporated into individual-based models to investigate how species distributions develop in response to environmental conditions.
The observation of animal movements on small spatial scales may provide a means to understand how large-scale species distributions are established from individual behavioral decisions (Humston et al. 2004; Roshier et al. 2008). In response to environmental heterogeneity, individuals within a population will alter their behavior to minimize physiological stress, minimize predation risk, and maximize foraging success (Wannamaker and Rice 2000; Lowe and Bray 2006). These individual behaviors subsequently influence the distribution and structure of populations (Humston et al. 2004). As fisheries management moves toward ecosystem-based approaches, it is important to understand how individual behaviors lead to large-scale species distributions and community development. The Summer Flounder Paralichthys dentatus is one example of an upper-trophic-level predator (Latour et al. 2008) that has the potential to influence community dynamics in coastal ecosystems and for which small-scale movement patterns are not well understood (Overholtz et al. 2000; Link et al. 2002).
Large-scale, seasonal migration patterns of Summer Flounder have been described from conventional mark-recapture studies conducted throughout most of the species' range (Poole 1962; Murawski 1970; Lux and Nichy 1981; Monaghan 1992; Desfosse 1995; Kraus and Musick 2001). During fall and early winter, Summer Flounder migrate to the continental shelf to spawn. After spawning is complete, individuals return to coastal habitats, where they reside during the spring and summer. These fish tend to return to the same inshore locations in subsequent years, but fish that emigrate tend to be recaptured in coastal habitats north of their release locations (Lux and Nichy 1981; Desfosse 1995).
In contrast to the well-documented seasonal migration patterns, relatively few studies have examined small-scale movements of Summer Flounder during the period of inshore residency. Studies using acoustic telemetry have found that Summer Flounder move in response to changes in tidal state (Szedlmayer and Able 1993; Henderson et al. 2014), low dissolved oxygen concentrations (Miller 2010), oncoming storms (Sackett et al. 2007; Henderson et al. 2014), lunar phase (Henderson et al. 2014), and decreased light levels at night (Miller 2010; Capossela et al. 2013; Henderson et al. 2014). Due to the limitations of acoustic telemetry, these studies could observe only meso-scale horizontal movements (hundreds to thousands of meters) but not the small-scale (1–100-m) vertical behaviors of these fish. Archival tags continuously record environmental information (i.e., depth and temperature) over long durations, making them an ideal technology for use in observing the small-scale vertical movements of fish (Block et al. 2001; Wilson et al. 2005). An understanding of small-scale behaviors is a critical component in developing individual-based models to examine how species distributions change in response to environmental conditions (Humston et al. 2004; Roshier et al. 2008).
We used temperature and depth data recorded with archival tags to observe the small-scale vertical movements of Summer Flounder during their residency within and dispersal from Chesapeake Bay. We focused our research on fish originating from Chesapeake Bay because it is the largest estuary in the Summer Flounder's range and is an important seasonal habitat for both juveniles and adults (Packer et al. 1999).
Tagging.—During 16 d in August and September 2009, we released 262 Summer Flounder with archival tags that were configured to record water temperature (range = -1°C to 40°C, resolution = ±0.03°C) every 60 min and depth (range = 1– 250 m, resolution = ±0.08 m) every 20 min. Fishing effort was distributed throughout the lower portion of Chesapeake Bay, but the majority (71%) of the fish were captured near the Chesapeake Bay Bridge Tunnel and at a few other sites near the bay mouth (Figure 1). Smaller percentages of fish were captured at other sites inside the bay (25%) or just outside the bay (4%). Summer Flounder were captured with hook-and-line gear, with the exception of four fish (1.5%) that were captured in a 13.7-m bottom trawl that was towed for 30 min (Bonzek et al. 2010).
Archival tags (Star-Oddi DST milli-L) were attached externally by using two nickel pins that pierced the dorsal musculature (Figure 2). The external attachment method was modified from a procedure used to attach similar archival tags to another flounder species (Cadrin and Moser 2006). On the nonpigmented side of the fish, small plastic discs were used to secure the pins. We allowed about 4 mm of space between the plastics discs and the epidermis of the fish to permit growth. A rubber earring backing was used to secure the plastic discs, and nickel pins were clipped and crimped around the earring backing to secure the archival tag to the fish. A T-bar anchor tag (Hallprint) was also inserted into the dorsal musculature as a secondary identification tool.
To maximize survival of fish after tag attachment and to avoid abnormal behaviors associated with application of a tag that was too large, only fish that exceeded 290 mm TL were tagged (range = 295–714 mm; mean = 413 mm; Figure 3). Archival tags measured 12.5 mm in diameter × 38.4 mm in length, weighed 9.2 g in air, and weighed 5 g in water. A length-weight regression from a previous study indicated that a 290-mm Summer Flounder in Chesapeake Bay weighs approximately 250 g (M.J.H., unpublished data). Therefore, we estimated that the archival tag weight was less than 4% of the weight of each fish.
Because it was necessary to recover the archival tags to retrieve the recorded temperature and depth data, we offered a $200 reward for returned tags. We also instituted an extensive advertising campaign at ports and fish processing houses throughout the mid-Atlantic coast to increase the probability of recovering tags. Data from recovered tags were downloaded and processed to remove measurements that were recorded prior to the tag's deployment date and after the tag's retrieval date.
Analysis.—Depth data recorded by the archival tags were used to examine the small-scale vertical movements of Summer Flounder during their residency in Chesapeake Bay. To ensure that the observed movements were not influenced by unusual behaviors due to the tagging procedure, we excluded depths that were recorded within 24 h of release. We also restricted our analysis to dates prior to October 15, 2009, to maintain a sample size of at least five fish. Restricting our analysis to these dates also reduced the probability that the analysis would be influenced by behaviors related to the spawning migration. Small-scale vertical movements were inferred from changes in depth between subsequent measurements from the same fish. While fish were presumed to be resident in Chesapeake Bay, we corrected for tide-related changes in depth by subtracting the predicted change in tidal amplitude from the observed depth change; tidal amplitudes were estimated using the Tides and Currents software program (Nobeltec, Beaverton, Oregon) and were similar throughout the study area. Because we did not know an exact location for each fish, we calculated the mean tidal amplitude for the lower Chesapeake Bay. Based on tidal corrections that were used to predict tides at various locations in the lower bay, we estimated that the mean tidal amplitudes were within 20 min of the actual tides experienced by our tagged fish.
Due to the importance of tidal state in the horizontal movements of Summer Flounder (Szedlmayer and Able 1993; Henderson et al. 2014), we elected to examine vertical movements during 1.5-h time intervals centered around four tidal stages: low, rising (the midpoint between low and high tide), high, and falling (the midpoint between high and low tide). The time of day for each tidal interval was assigned based on the sunrise and sunset times, which were obtained from Tides and Currents. If the midpoint of the tidal time interval occurred prior to sunrise or after sunset, the time of day for that interval was specified as “night”; otherwise, the time of day was specified as “day.” We did not include crepuscular periods in our model, as preliminary analyses indicated that movement behaviors taking place during crepuscular periods were similar to behaviors that occurred during the night (Henderson and Fabrizio 2011). Lunar phase was assigned using the moon phase output from the Tides and Currents program. Based on the percentage of the moon illuminated, the lunar cycle was divided into eight phases: new moon, waxing crescent, first quarter, waxing gibbous, full, waning gibbous, third quarter, and waning crescent.
We developed a generalized linear mixed model (GLMM) to investigate the effect of fish size, tidal state, time of day, lunar phase, and water temperature (°C) on the movements of Summer Flounder that were resident in Chesapeake Bay. Fish size (TL, mm) and temperature data were centered (i.e., the mean was subtracted from each observation) to reduce collinearity (Quinn and Keogh 2002). The response variable (mean depth change observed during each tidal interval) did not satisfy the linear model assumption of normality and was therefore transformed by using a Box-Cox transformation (Box and Cox 1964). The Box-Cox transformation was calculated by
where yi (λ) is the transformed response, yi is the untransformed response, and λ is a power parameter. The most appropriate value for λ was estimated with maximum likelihood using the solver function in Microsoft Excel 2010. We then used the MIXED procedure in the Statistical Analysis System (SAS Institute, Cary, North Carolina) to fit a repeated-measures GLMM to the transformed depth change data. A mixed model was used because our data included both fixed effects (i.e., fish TL, tidal state, time of day, lunar phase, and temperature) and an individual fish random effect (Littell et al. 2006; Bolker et al. 2008). A repeated-measures model was used because observations of the same fish recorded closely in time were serially correlated. The repeated measure in the model was the sequential number of the tidal interval since the start of the study. The serial correlation among the repeated measures was estimated with a specialized variance-covariance structure (Littell et al. 2006; Rogers and White 2007). We developed a three-step model selection procedure to determine the most parsimonious model: (1) we determined the appropriate random effects structure, (2) we built a global model, and (3) we selected the most parsimonious model from all subsets of the global model. A more detailed description of this procedure is provided in the Appendix.
Release and recapture information for 15 Summer Flounder that were archival tagged in Chesapeake Bay. Refer to Figure 1 for release and recapture locations.
Summer Flounder with archival tags were recaptured primarily within Chesapeake Bay shortly after release. Fifteen archival tags were recovered (6% recovery rate) from Summer Flounder that were at large from 1 to 810 d (Table 1). A t-test found no difference in mean TL between tagged and recaptured fish (t = 0.02, P = 0.98; Figure 3). With the exception of three fish, all recaptures occurred within 75 d of release. We chose to exclude one of these fish from the analysis because its unattached tag was recovered on a beach in North Carolina, and we were uncertain when the tag was shed. Of the remaining recaptured fish, seven (50%) were from a single release date (August 25, 2009), although they were recaptured on different dates (Table 1). These fish were released at two popular fishing locations near the bay mouth and were recaptured at five different locations (Figure 1). We could not discern an obvious reason (i.e., mean fish TL, water temperature, release location, or recapture location) for this apparent anomaly; therefore, we assumed that these fish were independent samples in our analyses.
Most of the archival-tagged Summer Flounder appeared to remain on or near the bottom for long periods (>2 consecutive weeks) while residing in Chesapeake Bay. During this time, observed depth changes reflected only tidal variation. Along-bottom intervals were occasionally interspersed with rapid changes in depth, when fish appeared to move either within the water column or to deeper or shallower habitats. In our statistical analysis, we did not differentiate between off-bottom and along-bottom movements, but our observations suggested that along-bottom movements to new locations were characterized by a 2–5-m shift in the mean depth followed by the resumption of a regular tidal pattern (Figure 4). In contrast, what we believe to be off-bottom movements were generally of a higher magnitude (2–10 m) and shorter duration (20–40 min), and the fish would come within 2–4 m of the surface. However, it is important to note that we could not observe movements that occurred between the depth measurements recorded every 20 min by the archival tags.
Summer Flounder movements modeled with a GLMM indicated that behavior varied among individuals and depended on fish size and environmental factors. The first-order autoregressive moving average was identified as the most appropriate variance-covariance structure, and the model also included the random effect associated with individual fish. The top-three models had a cumulative Akaike weight of 0.8 (Table 2), providing considerable evidence that the predictors within these models were the most influential on Summer Flounder movements observed with the archival tags. We present the results from the second-best model because it included all of the potential predictors of interest, which was not true for the top model, and we were primarily interested in qualitative understanding rather than quantitative predictions (Bolker et al. 2008). Furthermore, estimates of the predictors' effects were similar among the three models (Table 3); hence, our conclusions would be similar regardless of which of these models we selected. The model selection process indicated that movements were related to fish TL, tidal stage, time of day, and temperature,Table 2). Support was moderate for the relationship between depth change and fish TL, the TL × time of day interaction, and the time of day × temperature interaction (Table 2). Fish that were smaller than 400 mm TL exhibited larger depth changes at night than during the day, whereas time of day had little effect on the activity of larger fish (Figure 5). Summer Flounder movements tended to increase with increasing temperatures; this was especially evident during night (Figure 6) and during rising and falling tides (Figure 7).
Fixed-effects model selection table showing the predictors (fish TL, tidal stage [tide], time of day [TOD], lunar phase [lun], and temperature [temp]) and interactions included in the top-five models (selected using Akaike's information criterion [AICc]) for Summer Flounder depth changes. Also shown are the AIC difference (▵AIC; difference between the given model's AICc value and the minimum AICc value) and the Akaike weight (w), which represents the probability that the given model is the most parsimonious model. The relative importance of each predictor variable is the sum of the w-values for all candidate models containing that variable.
Restricted maximum likelihood parameter estimates for the top-three generalized linear mixed models that were selected to describe depth changes of individual Summer Flounder (see Table 2; temp = temperature; tide = tidal stage). Tidal stage estimates are relative to the rising tide. Time-of-day (TOD) estimates are relative to night. These models included an individual fish random effect and were fitted with an autoregressive moving average covariance structure. Note that the estimates are based on Box-Cox-transformed data, and the effects of TL and temperature are for centered data.
In addition to the information on small-scale depth changes, observations from the two archival tags that were at liberty for the longest durations provided observations reflecting what we believe to be Summer Flounder migration behavior. Examination of the depth histories from these two tags (Figure 8) indicated that the mean daily depth change was relatively small during the summer months but increased drastically during the fall and winter. The summer months correspond with the period in which Chesapeake Bay Summer Flounder are known to reside in inshore feeding habitats, whereas the fall and winter months are when they participate in the spawning migration (Desfosse 1995; Kraus and Musick 2001). Because this abrupt behavioral change appeared to be correlated with the timing of the spawning migration, we postulate that depth change behavior could be used to discern when these fish migrate. One of the two fish that were at liberty into the fall (tag number 199; 541 mm TL) was recovered on November 19, 2009, approximately 8.05 km (5 mi) east of Virginia Beach. Based on the depth history from the tag, this fish apparently initiated migration in mid-October 2009 (Figure 8a), when the mean water temperature recorded by the National Data Buoy Center ( www.ndbc.noaa.gov) at the mouth of Chesapeake Bay was 21°C and when the mean water temperature recorded by the archival tag was 20°C. The second fish (tag number 241; 454 mm TL) was recovered in November 2011 after 810 d at liberty. Unfortunately, tag memory was exceeded after 452 d, so no data were recorded after November 2010. This fish experienced temperatures as high as 27°C in summer and as low as 6°C during winter. Although we cannot be sure of where the fish was located at any time postrelease, depth change behavior that we considered to be indicative of migration was first exhibited by this individual in late November 2009. Fish 241 then resumed the relatively small depth changes indicative of inshore residence, maintaining this behavior from May through early November 2010 (Figure 8b). In both 2009 and 2010, this individual initiated behavior indicative of migration when temperatures recorded by the archival tag and at the Chesapeake Bay data buoy decreased to approximately 15°C. We emphasize that because archival tags cannot be used for geolocation, we can only infer that these depth change behaviors were associated with the offshore spawning migration.
Movements that were observed with archival tags revealed more intricate behaviors than could be observed with acoustic telemetry or conventional tagging data, providing useful information on the small-scale movements of Summer Flounder. Our analysis revealed that small-scale vertical movements were primarily related to fish TL, time of day, tidal stage, and temperature. The environmental covariates that were related to Summer Flounder vertical behaviors were similar to those related to meso-scale (hundreds of meters) horizontal movements observed with acoustic telemetry (Henderson et al. 2014). In fact, lunar phase was the only environmental covariate related to horizontal movement that was not also related to vertical movement; we suspect that this was due to an insufficient number of archival-tagged fish recoveries because we did observe behaviors during the lunar cycle that were similar to those observed with acoustic telemetry (Henderson 2012). The similarities between the vertical and horizontal behaviors imply that movements on these two axes are related. Both acoustic and archival telemetry observations could be incorporated into individual-based models to determine how small-scale behavioral decisions produce large-scale distributions of Summer Flounder. These individual-based models could then be used to understand the distribution and behavioral responses of Summer Flounder to environmental variability and climate change (Humston et al. 2004).
Based on the behaviors we observed with archival tags as well as the behaviors that were observed with acoustic telemetry (Henderson et al. 2014), we believe that Summer Flounder movements within Chesapeake Bay were primarily related to foraging. Summer Flounder migrate inshore during the spring and summer to feed and increase their energetic reserves for spawning, which occurs during the winter months (Packer et al. 1999). We suspect that the movements observed with archival tags were related to foraging behavior because hungry fish are more likely to be mobile to improve their chances of encountering prey (Gibson 2005). For example, mysids are the primary prey of Summer Flounder smaller than 375 mm (Latour et al. 2008; Buchheister and Latour 2011). Mysids are generally more active at night (Hurlburt 1957) and may elicit increased movements of smaller Summer Flounder during that time of day. In addition, smaller Summer Flounder may be more likely to forage at night to avoid any potential predators. However, Summer Flounder comprise only a minor component of the diets of large predators (Link et al. 2002). Thus, it is unlikely that predator avoidance greatly influenced their movements.
In contrast to the smaller Summer Flounder, the diets of larger Summer Flounder are dominated by fish (Latour et al. 2008), and this may require larger individuals to primarily utilize ambush tactics to capture prey (Staudinger and Juanes 2010). The larger Summer Flounder may be less active than smaller fish at night because fish that utilize ambush predation must remain sedentary for long periods. However, two of the larger tagged individuals (>450 mm) were more active at night than would be expected based on predictions from the GLMM. This result implies that Summer Flounder can employ multiple foraging strategies depending on available prey.
Based on our observations, we hypothesize that temperature influences Summer Flounder foraging behavior, but it is also possible that these fish move into different habitats as temperatures increase. Summer Flounder activity levels in Chesapeake Bay increased with rising water temperatures—a pattern that was most pronounced at night and during the rising and falling tides. In laboratory experiments, feeding rates were observed to increase with increasing water temperatures (Summer Flounder: Malloy and Targett 1991; Pacific Halibut Hippoglossus stenolepis: Stoner et al. 2006). The relationship between feeding rate and temperature is most likely an adaptation to meet increased metabolic requirements at higher temperatures (Malloy and Targett 1991; Fonds 1992; Claireaux and Lagardere 1999; Capossela et al. 2012). Hence, when temperatures increase, Summer Flounder may increase their activity levels at night and during rising and falling tides to feed on prey items that may be more available during these periods (e.g., mysids and zooplanktivorous fish; Taggart et al. 1989; Gómez-Gutiérrez et al. 2007). An alternative hypothesis is that when temperatures increase, Summer Flounder are more likely to move in an attempt to find a more suitable habitat. In this case, individuals may use tidal currents associated with rising and falling tides to move between habitats at night, when they are less likely to encounter predators. Tidal stream transport has previously been suggested as an energy-saving mechanism used by Summer Flounder to move between locations (Szedlmayer and Able 1993; Sackett et al. 2007; Miller 2010). To our knowledge, studies examining the diets of Summer Flounder relative to time of day and tidal stage are lacking; such studies would be useful in understanding the factors that drive the behavior of Summer Flounder in mid-Atlantic estuaries.
Although our sample size was small, observations from two archival-tagged Summer Flounder that were at liberty into the fall indicated that inshore residence behaviors could be differentiated from offshore migration behaviors based on mean daily depth changes. The mean daily depth changes of these two fish increased abruptly in the fall, corresponding with the period when Chesapeake Bay Summer Flounder are known to initiate their offshore migration (Desfosse 1995; Kraus and Musick 2001). One of the fish that exhibited this behavior change (fish 199) was recaptured in November approximately 8.05 km (5 mi) offshore of Virginia Beach. Based on the timing of the behavior change and the offshore recapture of this individual, we hypothesize that the observed change in behavior was related to the offshore migration. In 2009, the behavior of fish 199 changed in mid-October, while the behavior of fish 241 changed in late November—a difference of nearly 1.5 months. The single fish for which we had two consecutive years of data, fish 241, exhibited behavioral changes nearly 1 month earlier in 2010 than in 2009. If the observed behavior change was associated with offshore movement, then our anecdotal observations suggest that individuals do not rely on seasonal cues (i.e., photoperiod) to initiate the offshore migration. Based on emigration timing from a mid-Atlantic coastal lagoon, Capossela et al. (2013) postulated that Summer Flounder responded to decreasing temperatures for initiating offshore movements. In both 2009 and 2010, the archival-tagged fish in our study exhibited behavioral changes when water temperatures at the mouth of Chesapeake Bay declined to approximately 15°C, which is consistent with the hypothesis that these fish use temperature as a cue to initiate migration. An alternative hypothesis is that Summer Flounder do not initiate their offshore migration until they have achieved sufficient energetic reserves to sustain the migration. Unfortunately, our archival tag data were insufficient to discern the true motivation driving the initiation of offshore migration. Future studies on the cues used by Summer Flounder to initiate emigration from inshore habitats are warranted and would increase our understanding of how climate change could impact the migration behavior of these ecologically and economically important fish.
Steve Cadrin (University of Massachusetts Amherst) provided guidance on archival tag specifications and attachment techniques. Captain Chandler Hogg, Captain Craig Paige, and numerous volunteers from the Virginia Institute of Marine Science (VIMS) were instrumental in collecting Summer Flounder for the archival tagging project. We are grateful for comments provided by Rob Latour (VIMS), Tracey Sutton (VIMS), Harry Wang (VIMS), Steve Cadrin, and Susanna Musick (VIMS) on a previous version of this paper. Funding was provided by the National Oceanic and Atmospheric Administration's Sea Grant Population Dynamics Fellowship to M.J.H. and VIMS. This paper is Contribution Number 3349 of the Virginia Institute of Marine Science, College of William and Mary.
APPENDIX: GENERALIZED LINEAR MIXED MODEL SELECTION
We used a three-step process to select the most parsimonious generalized linear mixed model (Figure A.1). The first step was to identify the random effects structure of the model. The random effects structure includes the variance-covariance structure and the individual fish random effect, which represents the between-subject variability. We selected a preliminary random effects structure using models that contained the five main effects (i.e., fish TL, tidal stage, time of day, temperature, and lunar phase) and no interactions. Restricted maximum likelihood (REML) was used to compare models with different random effects structures (Pinheiro and Bates 2000; Zuur et al. 2007). The variance-covariance structures tested were variance components, compound symmetry, first-order autoregressive, and first-order autoregressive moving average. These variance-covariance structures allowed us to model the correlation with the repeated-measures response, which for these data was the depth change recorded by the archival tags. As recommended for repeated-measures models, the Kenward-Roger approximation was used to calculate the df and to adjust the estimated SEs (Littell et al. 2006; Bolker et al. 2008). We selected the preliminary random effects structure that best described the data as the model with the lowest value of Akaike's information criterion corrected for small sample size (AICc; Burnham and Anderson 2002).
After identifying the random effects structure, we developed a global model that included all of the main effects as well as any potential interactions (Zuur et al. 2007). To avoid testing thousands of models with every combination of main effects and interactions, we individually added all possible two-way interactions to the model with only the main effects. In this step, all models were fitted by using maximum likelihood (ML) and the preliminary random effects structure previously discussed. “Important” interactions were identified as those interactions that reduced the AICc value by more than 1 unit. We graphically examined the important interactions to determine whether such interactions were due to small sample sizes and simply reflected noise. Only informative interactions were added to the global model. The procedure used to evaluate the random effects structure (step 1) was repeated using the global model. This was necessary because a change in the mean structure (i.e., the fixed effects included in the model) affects the random effects model selection criterion calculated with REML (Littell et al. 2006). Thus, we validated that the correct random effects structure was used to develop the global model. We repeated this entire process until there was no difference between the random effects structure selected in step 1 and the global model selected in this step (step 2).
Our final step was to identify the fixed effects and interactions that best described the variation in Summer Flounder movements. Here, we used ML to fit models by using all possible combinations of main effects and the interactions identified in step 2 (Littell et al. 2006; Zuur et al. 2007). Once again, the df were estimated using the Kenward-Roger approximation, and AICc was used to select the most parsimonious model with the best fit to the data. The final model predictions shown in the Results were estimated by using REML (Zuur et al. 2007).
Bolker, B.M., M. E. Brooks, C. J. Clark, S.W. Geange, J. R. Poulsen, M. Henry, H. Stevens, and J.-S. S. White. 2008. Generalized linear mixed models: a practical guide for ecology and evolution. Trends in Ecology and Evolution 24:127–135.
Burnham, K. P., and D. R. Anderson. 2002. Model selection and inference: a practical information-theoretic approach. Springer-Verlag, New York.
Littell, R. C., G. Milliken, W. W. Stroup, R. Wolfinger, and O. Schabenberger. 2006. SAS for mixed models, 2nd edition. SAS Institute, Cary, North Carolina.
Pinheiro, J. C., and D. M. Bates. 2000. Mixed-effects models in S and S-Plus. Springer-Verlag, Statistics and Computing Series, New York.
Zuur, A., E. N. Leno, and G. M. Smith. 2007. Analysing ecological data. Springer, New York.