Genetic integration of behavioural and endocrine components of the stress response

The vertebrate stress response comprises a suite of behavioural and physiological traits that must be functionally integrated to ensure organisms cope adaptively with acute stressors. Natural selection should favour functional integration, leading to a prediction of genetic integration of these traits. Despite the implications of such genetic integration for our understanding of human and animal health, as well as evolutionary responses to natural and anthropogenic stressors, formal quantitative genetic tests of this prediction are lacking. Here we demonstrate that acute stress response components in Trinidadian guppies are both heritable and integrated on the major axis of genetic covariation. This integration could either facilitate or constrain evolutionary responses to selection, depending upon the alignment of selection with this axis. Such integration also suggests artificial selection on the genetically correlated behavioural responses to stress could offer a viable non-invasive route to the improvement of health and welfare in captive animal populations.


19
The vertebrate stress response comprises a suite of behavioural and physiological traits 20 that must be functionally integrated to ensure organisms cope adaptively with acute  (Stearns 1989;Roff 1992), behavioural ecology (Sih & Bell 2008)), but has not 82 been tested for explicitly using quantitative genetic approaches.

83
The most compelling evidence for genetic integration of behavioural and 84 physiological stress response traits to date comes from artificial selection experiments 85 on domestic animal populations (e.g., rainbow trout Oncorhynchus mykiss (Pottinger and 86 Carrick, 1999), Japanese quail Coturnix coturnix japonicus (Jones, Satterlee & Ryder   87 1994), house mice Mus musculus domesticus (Veenema et al. 2003b)). For example, lines 88 of rainbow trout selected for stress-induced plasma cortisol levels (Pottinger and 89 Carrick, 1999) experienced correlated evolutionary changes in behaviour (Øverli,

104
While selection experiments illustrate that genetic integration of behaviour and 105 physiology can occur, estimation of the genetic variance-covariance matrix (G) through 106 quantitative genetic modelling provides a complementary strategy that also allows 107 investigation of exactly how multivariate genetic variation is structured within 108 populations. In the context of the stress response, this should provide insights into both most commonly for suites of morphological traits (following, e.g., Cheverud 1982), but 114 that is equally applicable to any aspect of the phenotype (see, e.g., Hine

124
Here we estimate G for behavioural and physiological components of the acute 125 stress response in Trinidadian guppies (Poecilia reticulata). This enables us to 126 determine not only (i) whether these components are genetically integrated on the 127 major axis of genetic (co)variation (i.e., the first eigen vector of G, gmax), but also (ii) 128 whether the structure and orientation of this axis suggests variation in overall stress 129 responsiveness and/or 'coping style' (explained further below; Koolhaas et al. 2010; 130 Boulton et al. 2015 these stress response traits, and then determine the additive genetic contribution to this 153 (G; the genetic variance-covariance matrix for this suite of stress response traits). We 154 predicted that individual traits will be heritable and that G will contain strong genetic 155 correlation structure between behavioural and physiological components of the stress 156 response consistent with genetic integration. We also predicted that both behavioural 157 and endocrine components of the stress response would load on the major axis of 158 genetic variation in multivariate trait space, gmax. The 'stress coping style' model 159 (Koolhaas et al. 1999

216
The absence of a strong positive correlation between track length and area 217 covered ( Fig 1A) is notable and potentially biologically informative; if fish moved 218 randomly with respect to direction in the arena then area covered would increase 219 monotonically (to an asymptote at 100%) with track length. A possible explanation is 220 that a long track length arises sometimes from a (putatively) less stressed fish exploring 221 the arena (fish 1 in Fig 1B) and sometimes from a (putatively) more stressed fish 222 exhibiting a typical 'flight' response (fish 4 in Fig 1B). These two types of response can 223 be discriminated based on whether, in a given trial, higher track length is associated

318
There is strong evidence for phenotypic integration of Cortisol with behaviour at the 319 among-individual levels (see Table S5). To test for and characterise the hypothesised

336
The resultant estimate of G (Table 2)

501
Were one to accept this view, the trait loadings on gmax would suggest that thigmotaxis and freezing 502 behaviour are strong indicators of stress coping style. Notably, here we find continuous variation 503 along this major axis of genetic variation (see Figure 4-figure supplement 2), rather than the bimodal 504 distribution suggested by some of the coping styles literature (Koolhaas et al. 1999(Koolhaas et al. , 2007(Koolhaas et al. , 2010). An 505 alternative interpretation, however, is that gmax represents variance in stress response 'magnitude' 506 rather than coping style per se. This is because the putatively 'reactive' end of gmax might simply reflect 507 a low magnitude stress response rather than a particular style of 'coping' with a stressor. Such an 508 interpretation would help explain why fish (or genotypes) that produce lower levels of cortisol 509 following handling and confinement also have more 'exploratory' movement (higher relative area 510 covered) and reduced thigmotaxis (i.e., increased time in the middle) in the OFT. We note that the 511 distinction between 'style' and 'responsiveness' may be rather moot if, for example, fish with proactive 512 styles are also more responsive.

615
Heavily gravid females were then isolated in 2.8L brood tanks to give birth (and were returned to the 616 breeding tanks either after producing a brood or two weeks of isolation). Any offspring produced in 617 the breeding tanks were excluded from the experiment as maternal identity could not be positively 618 identified. For the following generations, after 3 months of isolation from males we moved females 619 into individual 2.8L tanks, with 1 male then circulated among 3 females. Males were moved between 620 females every 5-8 days. In this way, females did not have to be moved to brood tanks, and any 621 offspring could be assigned to mothers definitively. In this setup, offspring were moved to a separate 622 brood tank on the day of birth. Note that as the gestation period for guppies is approximately 1 month, 623 any brood produced by a female less than one month after exposure to their designated male was 624 recorded in the pedigree as having unknown paternity.

625
Within 24h of a female producing a brood we recorded her weight (g) and brood size. We kept 626 juvenile fish in full-sib family groups in 2.8L tanks before moving them to 15L 'growth' tanks at an behavioural trial, the individual tested was weighed and then moved to a temporary 'holding tank'.

642
Once a full group (as described above) had been tested, all were moved from the holding tank back to 643 their home tank. We replaced the water in the testing and holding tanks between groups to reduce the

729
Beakers were placed within cardboard 'chambers' to prevent fish from seeing each other or 730 experiencing outside disturbance. One fish was transferred every 30s, alternating across holding 731 tanks, such that all fish were in their beakers within 10min of the initial netting. After 60 mins in the beaker, each fish was removed by pouring its sample through a clean net into a second beaker, with 733 the fish then quickly checked to confirm ID and returned to the holding tank until the entire group 734 could be returned to its home tank.

735
We immediately filtered each water sample using Grade 1 filter paper ( All data handling and analysis was performed in R version 3.6.3 (Team 2020). We used the 'tidyverse' 763 set of packages for data handling and visualisation (Wickham 2017), and ASreml-R v4 (Butler 2021) 764 for fitting linear mixed effects models (as described in full below). We also used 'nadiv' for pedigree 765 preparation (Wolak 2012

771
When testing a single random effect (variance component), we assumed the difference to be 772 asymptotically distributed as an equal mix of 2 0 and 2 1 (denoted 2 0,1; Self and Liang, 1987; 773 Visscher, 2006            Eigen vector Proportion of variance