Individual differences in heart rate reveal a broad range of autonomic phenotypes in a free-living seabird population

Animals in the same population consistently differ from each other in their physiological and behavioural responses to the environment (Carere and Maestripieri, 2013), but the proximate causes of these differences, particularly the underlying neural mechanisms, remain poorly understood (Snell-Rood, 2013; Duckworth, 2015). The autonomic nervous system regulates wide-ranging physiological functions such as heartbeat, blood pressure, gastrointestinal activity, immunity, metabolism and reproductive functions to support the physical demands of behaviour (e.g. locomotion, eating, sex) or other internal changes (e.g. haemorrhage, infection), while maintaining homeostasis (Jänig, 2006; Kuenzel, 2015). Much of the phenotypic variation observed in wild populations may therefore be generated by individual differences in autonomic activity, but this has not yet been studied. Moreover, autonomic activity, because of its direct synaptic connections to the central nervous system, provides a window into an animal's brain while it activates behavioural responses (Thayer et al., 2012; Beissner et al., 2013) and can give new insight into the stability versus plasticity of the neural processes underlying physiology and behaviour.

The autonomic nervous system is composed of two independently regulated neural branches extending from the brain to the body, where they have largely opposing effects. The sympathetic branch drives the ‘fight-or-flight’ response (including an increase in heart rate), which helps prepare an animal for danger, and the parasympathetic branch promotes ‘rest-and-digest’ processes and self-maintenance (including a decrease in heart rate; Jänig, 2006; Kuenzel, 2015). Heart rate therefore reflects the balance between sympathetic and parasympathetic activity. In addition, the sympathetic and parasympathetic branches generate oscillations in heart rate at different frequencies, so heart rate variability can be analysed to separately measure the activity of each of the two autonomic branches (Carravieri et al., 2016; Müller et al., 2017).

We investigated, for the first time in a free-living animal population (the streaked shearwater, Calonectris leucomelas Temminck 1835, a long-lived pelagic seabird), whether individuals consistently differ in autonomic activity, over time and across contexts. We repeatedly recorded electrocardiograms from individual shearwaters, and quantified individual repeatability of heart rate and heart rate variability indexes that reflect parasympathetic and sympathetic activity, within and across years, and across the different contexts of stress, recovery from stress and a non-stress baseline, to assess the stability versus plasticity of individual autonomic responses. We found a wide range of autonomic phenotypes in the shearwater population that persisted across contexts and remained stable even across years.

We performed fieldwork on breeding adult streaked shearwaters in a large colony (84,000 breeding pairs, M.Y., unpublished data) on Awashima Island (38°18′N, 139°13′E) in the Sea of Japan during the chick-rearing season (August to October) in 2014 and 2015. Shearwaters build nests inside narrow burrows excavated in the soil on a steep coastal slope facing the sea. Burrows are typically 10–20 cm wide and ca. 0.5–1 m deep. Streaked shearwaters lay only one egg per season and the two parents contribute equally to parental care (Ogawa et al., 2015). During the chick-rearing period, adults spend the entire day at sea foraging for fish. They return to the colony after sunset to feed their chicks (spending most of this time inside their nest burrows), and depart for the sea again just before sunrise.

Fieldwork was performed at night, between 20:00 h and 04:00 h, when adults were present in the colony. Adults were captured from their nest burrows, identified by their permanent metal rings (unringed birds were given a ring), equipped with externally attached miniaturized electrocardiogram (ECG) data loggers and then returned to their nest burrows, so we could measure heart rate and parasympathetic (rest-and-digest) and sympathetic (fight-or-flight) indexes of heart rate variability during handling stress (just after return to the burrow, 0 min post-handling), during recovery from stress (after 20 and 90 min of resting in the burrow post-handling) and at baseline (120 min post-handling). We performed 55 tests in 2014 (49 individuals) and 140 tests in 2015 (69 individuals). Twenty-three birds were captured at least once in both years.

We used Little Leonardo ECG loggers (model W400-ECG, 21×109 mm cylindrical logger, 1 ms sampling interval, voltage range ±5.9 mV, 60 g, 2 GB memory) and Neurologger 2A ECG loggers (0.625 ms sampling interval, voltage range ±3 mV, 20 g, 1 GB memory, Evolocus LLC; for more details, see Müller et al., 2017). Although the two logger types differ in sampling interval, they both produce an ECG trace that shows a clear signal for each heartbeat; measurements of heart rate and heart rate variability did not differ between records from the two types of logger. Three lead wires extended from the ECG logger with small safety pins (electrodes) soldered to the ends that we subcutaneously attached to the breast skin of the birds (Yamamoto et al., 2009). We wrapped the wires around one side of the bird and secured the logger to the dorsal feathers with Tesa tape (Yamamoto et al., 2009). Compared with gluing, the use of subcutaneous pins has several advantages: it requires no feather removal, it results in quicker logger attachment (and therefore reduced handling time) and it causes no lasting damage (birds that are repeatedly tested and recaptured within a few days of a previous test exhibit no skin wounds or irritation from previous ECG logger attachments; Müller et al., 2017). It has therefore become the standard method for seabirds (e.g. Ropert-Coudert et al., 2006; Carravieri et al., 2016; Müller et al., 2017). We cleaned the pins and skin with alcohol wipes before attaching the loggers and replaced the pins several times during each season. After logger attachment (after a total handling time of 7–12 min), birds were placed back into their burrows for 2 h. At the end of each test, we retrieved the bird from the burrow, removed the logger and measured bill length with callipers, and with a Pesola spring scale (5 g accuracy). All fieldwork was authorized by the Japanese Ministry of the Environment. All procedures were approved by the Animal Experimental Committee of Nagoya University.

We analysed ECG data using Igor Pro version 6.37 (Wavemetrics, USA) in 5 min intervals, based on von Borell et al. (2007). In the PQRS complex, which is the cluster of graphical deflections that comprise a single heartbeat in an ECG wave, the R peak (occurring with the depolarization of the right and left ventricles of the heart), in particular, is very prominent in this species (Müller et al., 2017). We identified R peaks in ECG recordings primarily using the software Ethographer (Sakamoto et al., 2009), which permits smoothing of the wave and enhances the length of R peaks to facilitate peak detection. We manually identified R peaks when necessary. We created a data frame of the timing of each heartbeat (in milliseconds).

Using the RHRV package (http://CRAN.R-project.org/package=RHRV) in R (version 3.2.1), we filtered the beat positions to eliminate spurious beats from other prominent (non-R) peaks in the wave caused by muscle noise. The filtered dataset was then used to calculate inter-beat intervals (IBIs). Heartbeat positions plotted over time reveal oscillations in heart rate caused by the autonomic nervous system, which generates most of the heart rate variability (Müller et al., 2017). Oscillations occurring at a high frequency (between 0.3 and 2 Hz in this species, or every 0.5–3.3 s; Müller et al., 2017) reflect variability in heart rate generated by the parasympathetic nervous system and correspond to respiration: during inhalation, heart rate accelerates; during exhalation, heart rate slows, making oxygen delivery more efficient (respiratory sinus arrhythmia; Stauss, 2003; Taylor et al., 2014; Carravieri et al., 2016). The strength, or amplitude, of these oscillations in heart rate is therefore an index of parasympathetic activity. Oscillations occurring at low frequency (0.04–0.3 Hz in this species, or every 3.3–25 s; Müller et al., 2017) are generated by both the sympathetic and parasympathetic nervous system, and the amplitude is an index of combined sympathetic and parasympathetic activity (Malik et al., 1996; von Borell et al., 2007; Yamamoto et al., 2009; Carravieri et al., 2016).

We calculated these indexes from IBI data, using RHRV. We calculated the standard deviation of the differences between successive IBIs (‘rMSSD’), which reflects the amplitude of high-frequency oscillations and therefore parasympathetic activity, the standard deviation of all IBIs (‘SDNN’), which reflects the amplitude of low-frequency oscillations and therefore the combined sympathetic and parasympathetic activity (hereafter sympathetic+parasympathetic index), and the ratio between the SDNN and rMSSD (‘SDNN:rMSSD’) which therefore is approximately the sympathetic:parasympathetic ratio (Malik et al., 1996; von Borell et al., 2007; Kjaer and Jørgensen, 2011; Shaffer et al., 2014; Carravieri et al., 2016; see Müller et al., 2017, for more details about heart rate variability analysis in this species). We also used RHRV to compute average heart rate (which reflects the balance between sympathetic and parasympathetic activity, increasing or decreasing, respectively, when their level of activity increases) over the course of each 5 min interval.

Statistical analyses were performed using R (version 3.2.1). All indexes were log transformed to achieve normality except for heart rate, which was already normally distributed. To test how autonomic activity changed between acute stress just after handling (0 min post-handling), during recovery in the nest burrow (20 and 90 min post-handling) and at baseline (120 min post-handling), we performed mixed models on heart rate and heart rate variability indexes using the lmer function from the lme4 package (http://CRAN.R-project.org/package=lme4) and lmerTest (http://CRAN.R-project.org/package=lmerTest) to determine statistical significance. Time interval (0, 20, 90 and 120 min post-handling) was included in the model as a continuous predictor. Fifty-five birds were tested more than once, so we included individual ID as a random factor in all models (n=780 observations from 197 tests from 97 different individuals).

Repeatability is a standardized index that reflects the proportion of the variation in a phenotypic trait that comes from between-individual variation (Lessells and Boag, 1987). Thus, high repeatability values (closer to 1) indicate large and consistent differences between individuals in a trait, due to large between-individual differences relative to within-individual variability, whereas values closer to zero indicate that differences between individuals are small and intra-individual variability is high.

We calculated repeatability for each autonomic parameter (heart rate, rMSSD, SDNN, SDNN:rMSSD), for each time interval (0, 20, 90, 120 min), within and across years. Repeatability was calculated as the between-individual variance component divided by the sum of the within-individual and between-individual variance components, which were derived from linear mixed models (LMMs, with restricted maximum likelihood; Nakagawa and Schielzeth, 2010). We also included additional potentially confounding variables (timing in season, year) in the LMMs that could incorrectly inflate within- or between-individual variance estimates and so could bias our repeatability values; therefore, we calculated ‘adjusted repeatability’ sensu Nakagawa and Schielzeth (2010). The LMMs were constructed in the following way: they contained an autonomic index (e.g. heart rate) as a dependent variable, and the random factor ‘individual ID’, the fixed covariate ‘calendar date’ (timing in the season) and the random factor ‘year’. The variance of the random factor ‘individual ID’ represents the between-individual variance component, and the ‘residual variance’ component represents the within-individual variance. Calendar date corrected for changes in autonomic activity in all birds across the season that could bias results, either by artificially increasing or reducing consistency estimates for individuals if they were repeatedly tested at a similar time or at very different times, respectively, or artificially increasing between-individual differences if all tests for one bird were performed at a very different time in the season than all tests for another bird were performed. Year was included in case autonomic activity differed between years, as some birds were tested only in 2014 and others were tested only in 2015 (effects of calendar date and year on autonomic activity are reported elsewhere). The repeatability estimate was further corrected based on the recommendations of Nakagawa and Schielzeth (2010). As we used mean values of heart rate or heart rate variability over a continuous 5 min interval, and heartbeats from a 5 min interval are not independent, it was more appropriate to use a repeatability estimate of measurement means. Whether a trait was significantly repeatable or not did not differ depending on whether we used uncorrected repeatability or corrected repeatability for measurement means and therefore the type of repeatability estimate we use does not qualitatively change the results or interpretation.

We performed two sets of analyses: within-year repeatability and between-year repeatability. In our models producing variance components for within-year repeatability (n=45 individuals, 120–123 observations, see sample size details below), for birds that were tested in both years, we included only data from the year with the most tests for that individual. We also only included data from individuals that were tested two or more times within a year. The average time between consecutive tests from the same individual within a year was 6.72±6.13 days (range 1–30 days). The average time span between the first and last test from the same individual within the same year was 11.67±8.59 days (range 2–41 days).

For our analyses of between-year repeatability (n=23 individuals, 46 observations), if birds were tested more than once in one or more of the years, we selected data points from each year that were collected on the most similar calendar date. The difference between calendar dates of tests from the two years was 9.21±8.68 days (range 1–37 days between dates of tests).

Repeatability estimates were adjusted using the equation for unequal sample sizes from different individuals (Lessells and Boag, 1987). We estimated 95% confidence intervals (CIs) directly from a simulated distribution of repeatability generated by parametric bootstrapping (1000 iterations, described in detail by Faraway, 2006; as recommended by Nakagawa and Schielzeth, 2010). We performed likelihood ratio tests to test for statistical significance of variance of the random effect of individual ID (Bolker et al., 2009).

To visually compare the size of within- versus between-individual variance among autonomic indexes, we performed the same linear mixed models on the data after it had been standardized (x−mean)/s.d. and extracted variance components (Figs 1B and 2B,D). Data were standardized after removing outliers. Outliers were individuals that showed heart rates higher than 335 beats min⁻¹ during the recovery (20 or 90 min post-handling) or non-stress (120 min post-handling) phase, as such high heart rates indicate a stress response. Only one outlier was removed from the recovery phase at 20 min, and three outliers were removed from the non-stress phase at 120 min. The 24 h ECG recordings of birds at rest (incubating inside their nests) revealed no detectable circadian rhythm in heart rate or heart rate variability (Müller et al., 2017) so we did not correct for time of night in our analyses.

Shearwaters are very flexible in their autonomic responses, as evident in the large changes in sympathetic (fight-or-flight) and parasympathetic (rest-and-digest) activity across the contexts of stress, recovery from stress and non-stress: heart rate (which reflects the balance between sympathetic and parasympathetic activity) decreased from circa 300 beats min⁻¹ during the period of stress to circa 180 beats min⁻¹ at baseline (slope b=−0.8466, s.e.=0.0267, P<0.001; Fig. 1A). Heart rate variability indexes also revealed large changes in the activity of individual autonomic branches: parasympathetic activity increased sharply between the contexts of stress and baseline (rMSSD, the standard deviation of the differences between successive inter-beat intervals; Fig. 2A, b=0.0073, s.e.=0.0003, P<0.001), sympathetic+parasympathetic activity (SDNN, the standard deviation of all inter-beat intervals; Fig. 2C, b=0.0038, s.e.=0.0003, P<0.001) also increased to baseline, and sympathetic:parasympathetic balance (SDNN:rMSSD, b=−0.0035, s.e.=0.0003, P<0.001) decreased towards baseline, in line with the expectation that rest-and-digest activity increases and fight-or-flight activity decreases as an animal goes from a state of stress to a baseline resting state.

Despite this flexibility in autonomic activity, individual autonomic responses in a given context consistently differed from each other, revealing a wide range of autonomic phenotypes in this free-living population. Repeated recordings within the same year from 45 individuals showed that their heart rates consistently differed from each other and these differences persisted across different contexts (Fig. 1B, Table 1 f_H): during stress (0 min post-handling, Fig. 1C), during recovery from stress (after 20 and 90 min in the nest post-handling) and at baseline (after 120 min in the nest, Fig. 1D). Heart rate variability analysis showed that the wide range of differing autonomic phenotypes in this population, measured from heart rate, were driven by individual differences in parasympathetic rest-and-digest activity (rMSSD), which was also highly repeatable across the different contexts (Table 1, Fig. 2B), and showed repeatability even across years (Table 1). Sympathetic activity, in contrast, is elevated only during stress in this species (Müller et al., 2017) and did not show consistent individual differences mainly due to high within-individual flexibility – the SDNN index, which reflects combined sympathetic+parasympathetic activity, produced significant repeatability only after birds had fully recovered from stress at a time when sympathetic activity is negligible (120 min; Table 1, Fig. 2D).

The sympathetic and parasympathetic branches are regulated – through direct synaptic connections – by two, mostly separate, and inversely activated, brain networks (Thayer et al., 2012; Beissner et al., 2013) that are well studied in humans (Fox et al., 2005; Buckner et al., 2008), are present in other mammals (Rilling et al., 2007; Vincent et al., 2007; Lu et al., 2012) and have functionally and anatomically homologous structures in birds (Shanahan et al., 2013). These brain networks regulate not only physiology (via the autonomic nervous system) but also behaviour (via the somatic nervous system). The activity of the two autonomic branches, measured from heart rate and heart rate variability, therefore provides a real-time window into the activation of these brain networks and the mental state of the animal while it interacts with its environment (Jänig, 2006; Beissner et al., 2013). Examining consistency versus plasticity of individual autonomic responses over time and across contexts thus provides insight into the stability versus flexibility of the neural circuitry that regulates physiology and behaviour. Although we found the birds to be very flexible in their autonomic responses between contexts (as evident in the large changes in sympathetic and parasympathetic activity between stress and non-stress phases; Figs 1A and 2A,C), individual birds consistently differed from each other in their responses in the same context (Figs 1B and 2B), demonstrating significant stability in the activation of one or more brain networks in a given situation. This stability was not evident for the network regulating the sympathetic branch, which is activated during focused attention on specific tasks/events including threats (Thayer et al., 2012; Beissner et al., 2013): sympathetic activity exhibited high within-individual variability in sympathetic-active contexts (stress and recovery) even within years (Fig. 2D, Table 1 SDNN). In contrast, activation of the ‘default-mode’ brain network, which regulates the parasympathetic branch and a mental state of broadly tuned outward watchfulness and monitoring of the external environment (Thayer et al., 2012; Beissner et al., 2013), appeared to be very stable within individuals even across years (Fig. 2B, Table 1 rMSSD).

Our demonstration of consistent individual autonomic phenotypes in a wild free-living animal points to autonomic activity as a key neural mechanism driving broad phenotypic variation in natural populations, which has important ecological and evolutionary implications. Many eco-physiological and life history traits and trade-offs hinge on the allocation of limited resources. Because the autonomic nervous system regulates metabolism (reflected in heart rate; Romero and Wingfield, 2016), distinct autonomic phenotypes can mediate different solutions to resource allocation trade-offs. Divergent phenotypes in the same population are often favoured in competitive environments and can make way for adaptive individual niche specialization (Bergmüller and Taborsky, 2010; Dall et al., 2012). Continued success of such phenotypic variants is a step towards evolutionary change (Fusco and Minelli, 2010). Autonomic phenotypes in humans are 50% heritable (Neijts et al., 2015) and show a partial genetic basis in laboratory animals as well (Koolhaas et al., 1999; Korte et al., 1999; Kjaer and Jørgensen, 2011). Our study identifies the autonomic nervous system as a potentially important mediator of life history evolution.

We thank Masaki Shirai, Sakiko Matsumoto and Giacomo Dell'Omo for assistance with fieldwork.