Using skin temperature and activity profiles to assign chronotype in birds

Abstract

Chronotypes describe consistent differences between individuals in biological time-keeping. They have been linked both with underlying variation in the circadian system and fitness. Quantification of chronotypes is usually by time of onset, midpoint, or offset of a rhythmic behaviour or physiological process. However, diel activity patterns respond flexibly to many short-term environmental influences, which can make chronotypes hard to identify. In contrast, rhythmic patterns in physiological processes, such as body temperature, may provide more robust insights into the circadian basis of chronotypes. These can be telemetrically recorded from skin-mounted, temperature-sensitive transmitters, offering minimally invasive opportunities for working on free-ranging animals in the wild. Currently, computational methods for deriving chronotype from skin temperature require further development, as time series are often noisy and incomplete. Here, we investigate such methods using simultaneous radio telemetry recordings of activity and skin temperature in a wild songbird model (Great Tit Parus major) temporarily kept in outdoor aviaries. Our aims were to first develop standardised selection criteria to filter noisy time series of skin temperature and activity, to second assign chronotype based on the filtered recordings, and to third compare chronotype as assigned based on each of the two rhythms. After the selection of rhythmic data using periodicity and autocorrelation parameters, chronotype estimates (onset and offset) were extracted using four different changepoint approaches for skin temperature and one approach for activity records. The estimates based on skin temperature varied between different approaches but were correlated to each other (onset: correlation coefficient r = 0.099–0.841, offset: r = 0.131–0.906). In contrast, chronotype estimates from skin temperature were more weakly correlated to those from activity (onset: r = −0.131–0.612, offset: r = −0.040– −0.681). Overall, chronotype estimates were less variable and timed later in the day for activity than for skin temperature. The distinctions between physiological and behavioural chronotypes in this study might reflect differences in underlying mechanisms and in responsiveness to external and internal cues. Thus, studying each of these rhythms has specific strengths, while parallel studies of both could inform broadly on natural variation in biological time-keeping, and may allow assessment of how biological rhythms relate to changes in the environment.

Background

Diel activity patterns are widespread across all living organisms. For example, plant species rhythmically open and close their leaves [1] and zooplankton migrate vertically through the water column between the photo- and scotophases [2]. These phenomena are broadly synchronous to the daylight patterns of the 24-h day (but see [3]). Animals show a wide range of diel activity patterns, including crepuscularity, nocturnality or diurnality, as well as ultradian rhythmicity spread across day and night [4, 5]. Diel rhythms not only occur in behaviour, but also in physiological processes, such as metabolic rate, immune defence and body temperature. For instance, in endotherms body temperature shows a diel sinusoidal pattern of rise and fall, with higher temperatures during the active phase, e.g., during the day for diurnal organisms (reviewed in [6, 7], but see also [8,9,10]). Even when temporal environmental cues are experimentally removed, most organisms retain an approximately diel (i.e., circadian) rhythmicity of behaviour and physiology due to the internal circadian clock [11,12,13,14]. The circadian clock system, in turn, coordinates animals’ responses to the environment, and is itself based on cellular cycles of clock proteins and several interacting physiological pathways [15, 16]. In natural environments, rhythms of most species entrain to the natural diel light–dark cycle. However, many further factors of both the external and internal (e.g., state) environment can affect timing [17].

The timing of activity and physiological processes is often essential for the success, health and well-being of organisms, including humans [18]. Therefore, biomedical sciences explore in particular the circadian basis of diel rhythms, whereas ecologists are increasingly interested in environmental context and fitness implications of rhythmicity [19, 20]. Timing can vary considerably between individuals within populations (e.g., [21,22,23]). When diel timing of trait expression is consistent within individuals, this can be defined as a chronotype. Chronotypes are scaled on a continuum from early to late activity relative to their conspecifics, based on characteristics, such as onset, offset or midpoint of activity (e.g., [24, 25]). These individual differences are genetically and environmentally determined, but may also be obscured by additional variation [17]: first, chronotype is shaped by genetic differences in the circadian system and in linked physiological systems, including those that regulate environmental responses (e.g., light sensing). Second, due to the environmental responsiveness of the circadian system, local factors (e.g., light exposure) and behavioural context (e.g., reproductive stage) can contribute to individual differences. Third, additional individual variation, i.e., deviations from the chronotype, can arise from differences in the individual’s internal state (e.g., energy expenditure, health, age) and from short-term adjustments (temporal flexibility or “masking”), which can allow individuals to cope with sudden, unpredictable changes in the environment [17].

It is increasingly recognised that as temporal environments are modified by indirect or direct anthropogenic disturbance, e.g., via artificial light or hunting, chronotype may shift, and that fitness implications can be major [26, 27]. However, the relative benefits and costs for the different chronotypes are mostly unclear, and it is also unclear whether the alignment between behaviour and physiology remains intact when diel rhythms shift [28]. Therefore, there is a need to better understand both, behavioural and physiological, aspects of chronotype. So far, chronotypes have been assigned to wild individuals using different types of remotely collected activity measurements: overall activity [29, 30] or sleep [31, 32], but also diel rhythms of feeding [33], parental care [34, 35] or incubation [24]. More recently, fast technological developments improved remote recording not only by making data loggers and radio transmitters lighter, but also smarter. In particular, physiological measurements of heart rate, EEG or body temperature by the tags can also be transmitted over larger detection ranges making techniques developed in the 1980s for laboratory use [36, 37] suitable for small free-ranging animals in the wild [38]. Body temperature in particular is widely assessed in the study of biological rhythms. While time series on activity and physiology are thus now readily available, a limiting factor for understanding chronotype are methods for extracting meaningful descriptors of temporal behaviour. Existing analytical methods have been mostly developed for studies under laboratory conditions, which typically yield high-quality data series (e.g., [10, 37, 39]). However, remotely collected data on chronotype exposed to multiple abiotic and biotic factors [40] are usually of much poorer quality. Therefore, the overarching aim of this study is to develop computational models that facilitate the use and analysis of information on chronotype from remotely collected data, in particular on body temperature.

Analyses of physiological and behavioural cycles are based on the extraction of different curve characteristics [18, 41, 42]. As first steps, rhythms are sometimes distinguished from noise by testing for autocorrelation between subsequent data points [43]. Autocorrelation of rhythmic data is particularly strong between the data points with specific time lags of one full cycle of behaviour or physiology (i.e., period length). Thus, for diel rhythms, autocorrelation should peak at 24 h. Overall levels of the traits and the strength of the rhythms can be assessed using mean (mesor) and amplitude (difference between maximum or minimum with mean level) of the curves, respectively. This information can be used for exploration of individual variation and flexibility of the rhythmic traits. Phase markers (points repeated in each cycle of the rhythm) can be derived as reference points to characterise a rhythm, in particular onset, mid-point, offset, minimum, maximum, or inflection points of increase or decrease. These phase markers are then related to environmental reference points, for example sunrise or midnight, and can be used to assign chronotype.

Various customised time series applications have been developed to extract curve characteristics described above, of which many are open source and based on the statistical software environment R [44]. For instance, the lomb (2.0 Ruf 2021 based on [45]) and cosinor (Sachs 2014 based on [46]) packages are useful for basic analyses directly inside R, whereas more customised applications such as Chronoshop [47] and RhythmicAlly [48] provide user-friendly interfaces showing informative diagrams and can analyse multiple data sets simultaneously. Even though many settings and features are available in the latter, the programs have their limitations in the flexibility of adjustment to different data types, in the curve characteristics to be extracted, and in the required engagement with circular data. Furthermore, curves can deviate considerably from standard models, for example for activity by sudden shifts, or for body temperature by a spectrum from sinusoidal curves to time series with plateaus [6, 10, 49, 50]. Therefore, many studies use generic procedures, such as visual data extraction [24] or employ changepoint analyses to extract phase markers of interest in their specific type of data [29, 51, 52]. Depending on the research question and the data available, changepoints can be assigned based on moving averages or on changes in variation between subsequent data points. These various methods likely differ in the chronotype descriptors they produce, and in their suitability for assigning chronotype to time series on behaviour and especially physiology of wild animals. To our knowledge, this has so far not been investigated.

Here, we make use of temperature-sensitive radio telemetry that continuously and simultaneously records activity and skin temperature, applied to a small songbird, the Great Tit (Parus major, Linnaeus 1758). Telemetry in the wild often results in time series with gaps, depending on the distance of animals to the receivers [45]. As in many studies of wild animals, we examined peripheral body temperature from the skin, rather than core body temperature, to reduce invasiveness and to simplify remote data collection without a need for recapture [53, 54]. Data on skin temperature carry particular challenges. Although peripheral body temperature correlates with core temperature, these measurements are confounded by effects of ambient temperature, wind and solar radiation [54, 55]. Moreover, absolute temperatures from external skin measurements can be confounded by differences in the attachment of the transmitters (i.e., tightness and insulation). This can make mean and amplitude estimates unreliable [54, 56], whereas changepoints may be more robust, because they refer to intra-individual data. These challenges mandate cautious data selection and analyses to avoid an arbitrary choice of suitable time series and subjective assignment of chronotype.

Methods

Aims

Using skin temperature and activity records from the transmitters, we thus aim to develop a standardised protocol to filter and process skin temperature data to assign chronotype from a physiological rhythm. To develop these analyses, we use data collected under semi-natural conditions in outdoor aviaries, where long recordings are available for some individuals. Our protocol includes (1) the proposal of an automatic and generalised protocol to select for rhythmic data based on periodicity and autocorrelation analyses and (2) the exploration of four different changepoint approaches to extract chronotype estimates from skin temperature patterns, i.e., for the diel onset and offset (i.e., increase and decrease) in temperature. In addition, ambient temperatures will be investigated to account for potentially confounding effects on skin temperature measurements. (3) We will then compare chronotype estimates from physiology to those from behaviour, using activity onset and offset derived from a well-developed Behavioural Changepoint Analyses (BCPA, [29]). Through these protocols, we hope to stimulate use of the rapidly increasing time series data for investigating chronotype of behaviour and physiology, and for exploring potential discrepancies between them that may inform on responses to environmental change.

Data collection

20 Great Tits (5 females and 15 males) were caught during roost checks on 1 November 2012 in southern Sweden 25 km east of Lund, individually ringed, and brought to the field station Stensoffa (55°42ʹN, 13°27ʹE). All birds were kept in four outdoor aviaries exposed to natural light–dark cycles and outside temperature fluctuations. Ambient temperature was measured every 10 min in the inside ceiling of the aviary roof using a small temperature logger (iButton DS1922-L, Maxim Integrated Products, CA, USA; accuracy: ± 0.5 °C). Each individual bird was provided with a nestbox for roosting and had ad libitum access to food and water. Housing conditions are further detailed in Nord et al. [54]. As part of another experiment, individuals received a temperature sensitive PIT tag implanted subcutaneously at the right flank (for details see [54]).

Telemetry

On 29 November 2012, 19 individuals were tagged with a temperature-sensitive radio transmitter (PicoPip, Biotrack, Dorset, UK; temperature sensing: ± 2 °C accuracy and 0.1 °C resolution; 0.55 g, 3.1% of mean body mass) after which data on skin temperature and activity were recorded for up to 19 days. In temperature-sensitive transmitters, the duration of the interval between subsequent emitted signals changes with recorded temperature, such that shorter intervals indicate higher temperatures. The transmitters were individually calibrated from 22.9 to 48.1 °C by the manufacturer shortly before deployment. The externally recorded temperatures from the transmitters are further referred to as skin temperatures. Variation in signal strength, generated by movements of the animal relative to the receiver, was used to assess the individual’s activity (e.g., [29, 30, 55, 57]). To enlarge the attachment surface, the transmitters were sewed onto a 10 × 20 mm piece of cloth. For attachment, feathers on the back of the bird were brushed to the sides. After clipping feathers of a patch equalling the cloth size, a transmitter was attached to the skin and the clipped-down feathers, using a thin layer of eyelash glue (DUO Eyelash Adhesive, American International Industries, CA, USA) and small drops of superglue at the cloth edges (Sekunden Alleskleber, UHU, Buehl, Germany). Following attachment of the transmitter (Additional file 1: Figure S1), feathers from adjacent plumage parts mostly covered the transmitter (for details see [54]).

The transmitter signals were recorded using a telemetry receiver (SRX400A, Lotek Wireless, Newmarket, ON, Canada) with a Yagi antenna (Lintec flexible 3-element Yagi antenna, Biotrack, Dorset, UK) on the rooftop of the building. The receiver was set to record each frequency (one frequency per transmitter: n = 19) for 8 s, resulting in data intervals of about 2.5 min per bird. The receiver filtered automatically for noise signals of > 20 ms length and recorded the calculated temperature from 3 consecutive signals.

Data selection and diel rhythmicity

R (version 4.2.1) and RStudio (version 2022.07.01) were used for data processing and further analyses (documentation of the R Script used for data selection and changepoint analyses are reported in Additional file 2). The first days (including 3 December) have been excluded from the data set due to disturbance for other studies (see [54, 58]).

Ambient temperature was checked for periodicities using the Lomb–Scargle periodogram (lsp function of lomb package) and for autocorrelation (acf function) across the complete time series. Because no clear diel rhythmicity in ambient temperature could be found in our data, it was no longer considered in further skin temperature analyses.

We then used the complete time series on signal strength and skin temperature of each individual for data exploration. The skin temperature data were filtered to retain only records higher than 23 °C and lower than 40 °C to account for the calibrated temperature range and exclude unlikely high skin temperatures in winter [54, 59]. This excluded a malfunctioning transmitter (frequency 173.951; 4.5% of the data set) and an additional 19.5% (10,561 of 54,111) of records from the remaining transmitters. The remaining data were assigned to 10 min bins by averaging values for skin temperature and signal strength. To indicate whether an animal was active, changes in signal strength were calculated as absolute deviation between two consecutive 10 min bins (e.g., [55]). We then used these deviations for further analyses, rather than the raw data, to account for between day differences in mean signal strength. The skin temperature data of each individual were examined for rhythmicity using a 3-day sliding window for calculation of Lomb–Scargle periodograms (lsp function of lomb package). The days when sliding windows showed diel rhythmicity (defined by significant period length between 23 and 25 h; further referred to as diel days) were retained (for details see Additional file 2: Selection of Rhythmic Days). When sliding windows stopped showing diel rhythmicity, the day that caused deviation was excluded (Additional file 1: Figure S2). Furthermore, remaining data were analysed for positive autocorrelation (acf function) at a time lag of 24 h (Fig. 1D) and periodicity (Additional file 1: Figure S3A, Additional file 2: Autocorrelation and Periodicity). Periodicity and autocorrelation were also scrutinised for the activity data (Fig. 1E, Additional file 1: Figure S3B), which resulted in no further exclusion of data points. For creating bird-specific diel profiles of skin temperature and activity, we averaged the binned data per daytime (24 h = 144 times of day) over all included diel days and calculated standard errors for these mean skin temperature estimates as measurement for their accuracy (Additional file 2: Average Diel Skin Temperature Profiles). Subsequently, highly inaccurate skin temperature estimates (i.e., SEMs of more than three standard deviations from the mean accuracy: 21 of 1292 data points) were excluded from the diel profiles as well as estimates with only one observation per daytime (i.e., accuracy/SEM = NA: 27 of 1292 data points). All highly inaccurate estimates belonged to the time series of one particular bird (individual 173.258 M) because of low numbers of diel days available and some extreme noise at night, resulting in its exclusion due to large data gaps (> 25% of missing data, for details see Additional file 2: Average Diel Skin Temperature Profiles). The final data set consisted of eight individuals with 4–16 days of recording (Table 1).

Fig. 1

Time profile of ambient temperature, and rhythmicity of skin temperature and activity. Upper three panels show time profiles of ambient temperature (A), skin temperature (B) and activity (C) (10-min binned) against date and time for an exemplary Great Tit (frequency 173.204, male) kept in an outdoor aviary during winter. Bottom two panels show corresponding autocorrelation plots for skin temperature (D) and activity (E), including a red vertical line to indicate 24 h

Table 1 Final selection steps for inclusion of birds’ records into chronotype analysis. Sample size and applied criteria (i.e., diel periodicity of 23–25 h and presence of an autocorrelation peak at 24 h) for skin temperature and activity for the ten individual Great Tit [Female (F) and Male (M)] retained after initial filtering steps. Diel days available refers to the number of rhythmic days included while filtering skin temperature in the 3-day sliding window approach (incl. the total number of recorded days in brackets). The remaining nine filtered-out individuals with zero diel days are not shown. Two additional individuals (marked in bold font) were excluded from the chronotype analyses as they did not fulfil all selection criteria (173.391 M had no positive autocorrelation at 24-h peak and 173.258 M had > 25% of data missing in the diel skin temperature profile (173.258 M)

Overall diel profiles (means of all individuals) indicated that skin temperature showed a broadly sinusoidal pattern with a tendency of plateau phases around minimum and maximum temperatures followed by an increase in the early morning before sunrise and a decrease in the afternoon, respectively (Fig. 2A). In contrast activity data showed a sudden increase around sunrise and decrease around sunset (Fig. 2B).

Fig. 2

Diel profiles for skin temperature (A) and activity (B) for each individual Great Tit (unsmoothed data in different colours) and mean (black solid curves) ± standard error (black dashed curves) for the full data set. Blue curve shows the fitted sinusoidal curve (with 3 harmonics) for skin temperature. Horizontal red lines show overall mean (levels). Daylight phases are indicated in yellow (sunrise–sunset)

Before assigning chronotype (i.e., onset and offset) based on skin temperature and activity, further preparations were necessary. A first step was to delineate regions of interest (here called “windows”) during which we searched for onset and offset. This was done to maximise accuracy and functionality of the search algorithms for detecting chronotype. We used the overall diel profiles (Fig. 2) to define the centre of the windows, i.e., one for onset and one for offset, as applied previously for seasonal timing [60]. For skin temperature, identifying the centre first required smoothing of the data series for noise reduction, for which we fitted a three-harmonic sinusoidal curve (Fig. 2A: blue curve). The choice of a three-harmonic curve was based on a visual comparison of the relative advantages of different numbers of harmonics (for one versus three harmonics, see Fig. 3), whereby fewer harmonics diminish noise but also diminish temporal resolution [61].

Fig. 3

Bird-specific skin temperature profiles, centred around minima and maxima to estimate onset and offset. Shown skin temperature curves for onset (A) and offset (B) are smoothed with 1 (red, dashed) and 3 (blue, solid) harmonics. Graphs also show means (black dots) ± standard error (grey bars) per 10-min bin. Daylight phases are indicated in yellow (sunrise–sunset). Sex: female (F) and male (M))

From the smoothed skin temperature time series, we then selected the times of minimal (03:10 h) and maximal (12:10 h) temperatures for centring two 24-h analysis windows for onset and offset, respectively. 24-h windows were used to avoid a priori limitation of onset and offset times, given the relatively large variation in the timing of diel skin temperature patterns observed in our data. For activity chronotype, we used unsmoothed data. As centre points of the windows, we used overall activity onset and offset, defined by the intersections of the overall diel activity curve with the overall mean activity level (onset at 07:30 h and offset at 16:20 h; Fig. 2B: red line). We set window size to 6 h in total, and thus narrower than for skin temperature, to avoid artefacts from quick but intense activity level changes during night and day. Using these settings, the next step was to derive bird-specific times of onset and offset.

Chronotype estimates from skin temperature

For a comparison of methods to assign chronotype from skin temperature, we surveyed previous applications of identifying changepoints (here, onset and offset) in various types of time series that showed similar curve characteristics (i.e., increases and/or decreases) as seen in temporal patterns of diel body temperature [52] but also of torpor [51, 62], seasonal and diel activity [29, 30, 60] and moulting [63]. We identified three approaches useful to our study. These included two approaches based on segmented linear regressions (as in [52, 60]), whereby one (Broken Stick Regression) fits a single changepoint of slopes [51], whereas the second (Structural Changes) generates multiple changepoints of slopes, from which an optimal changepoint can be selected. A third approach estimates changepoints from distributions (Behavioural Changepoint Analyses, BCPA, by [29]), whereby multiple changepoints are generated, and the one that results in the best model fit is selected. Since these methods sometimes failed to pick up changepoints that visually aligned with the onset (i.e., start of increase in skin temperature) and offset (i.e., start of decrease), we also developed a fourth one specifically for our data type (named Changepoint here: see below).

Various R packages are available to extract these changepoints. For the Broken Stick Regression we chose the lm.br function (from lm.br package by Adams 2021 based on [64]) used in Lewden et al. [51]. Different type settings (i.e., line–line, threshold–line, line–threshold) can be adjusted and applied in combination with a selected analysis window depending on the customer’s needs. Confidence intervals can be extracted directly from the model. Here, we fitted two lines (type line–line) to assign one changepoint close to minimum or maximum for onset and offset, respectively.

For Structural Change the function breakpoints (from the strucchange package by Zeileis et al. 2019 based on [65, 66]) was used to fit multiple models with several lines to the data. The function, then, selected the model with the optimal number of changepoints (with the lowest Bayesian Information Criterion (BIC)). A maximum number of fitted changepoints could be enforced (using breaks setting) resulting in poorer fits (higher BIC) and slightly different changepoint estimates. To avoid poorer fits, we left breaks undefined, and instead used the best fitting model to select one changepoint that was closest to the temperature minimum or maximum for onset and offset, respectively. Confidence intervals for the changepoints can be calculated using the confint function.

In the BCPA, models with two normal distributions, and thus, one changepoint, are fitted (fitdistr function) for each potential timepoint in the observed window [29, P. Capilla pers. comm.]: one distribution for the data before a potential changepoint and another one for afterwards. The Log likelihood is calculated (logLik function) for each distribution and the Akaike Information Criterion (AIC function) is obtained for each whole model. The model with the lowest AIC is selected and the associated changepoint is used as a chronotype estimate. As AIC showed two minima in a 24-h window, i.e., one before and one after the skin temperature peak, the chosen changepoint (with lowest AIC) was ambivalent before or after minimum or maximum. To detect increase (i.e., onset) or decrease (i.e., offset) in temperature, we reduced the window to 12 h using the second part of the analysis window only (i.e., 12 h from minimum or maximum). BCPA applied to skin temperature data is further referred to as BCPA_Tskin.

Our new method, Changepoint, was based on the slope of change in the temperature data. This allowed us to extract the start of increasing temperature in the morning, which might reflect the rewarming of the body in preparation for the day. Likewise, we extracted the start of decreasing in the afternoon, which could be interpreted as the end of the physiological day, where thermoregulation is (partly) turned off allowing for passive cooling. In particular, we looked at the deviation of temperature in subsequent bins over time, i.e., relative to the previous 10-min bin. We calculated the slopes of the bins following the temperature minimum or maximum for onset and offset, respectively. One can apply different thresholds for the slopes to estimate the changepoint, ranging from zero to maximum slope, and thereby adapt the algorithm to account for temperature plateaus with low slopes. Here, we used the first increase or decrease (slope > 0 or < 0) for estimating skin temperature onset and offset (for details on all changepoint analyses see the R script in Additional file 2 and 3).

Chronotype estimates from activity

For assigning chronotype parameters, onset and offset of activity were determined from the distinct changes of variation in signal strength using BCPA [29, P. Capilla pers. comm.]. In our data, changes in variation were extracted from the previously calculated activity levels (i.e., absolute deviation of signal strength, see above and the R script in Additional file 2 and 3). Changepoints extracted from the activity data are further referred to as activity chronotypes to distinguish them from the BCPA_Tskin output.

Comparison of chronotype estimates of skin temperature and activity

For comparison of methods and chronotype estimates from skin temperature and activity, we computed means, standard errors and a Pearson correlation matrix (cor function) for onset and offset separately. The correlation matrix was fed with estimates from Broken Stick Regression, Structural Change, Changepoint and BCPA_Tskin for skin temperature as well as activity chronotypes of all eight individuals.

Results

Diel rhythmicity of ambient temperature, and of activity and skin temperature

Ambient temperature was explored for diel rhythmicity to assess confounding effects on rhythmic patterns of skin temperature. The diel temperature profile showed no clear 24-h pattern, but between day variation in profiles (Additional file 1: Figure S4). The data points were decreasingly autocorrelated with increasing time lags up to 48 h (acf = 0.37 for 24 h), instead of periodic autocorrelation peaks at 24 h, 48 h, 72 h etc. as seen in the diel rhythms of the birds (compare Fig. 1D to Additional file 1: Figure S5A). The Lomb–Scargle analysis showed weak, but significant, periodicity at a time lag of 1 day and additional significant periodicities, at 2 and 3.4 days, with increasing normalised power reaching infinity for periods longer than 4 days (Additional file 1: Figure S5B).

Skin temperature and activity data were filtered using inclusion criteria described above resulting in eight individual data sets with a mean period length of 24.04 ± 0.05 h (mean ± SEM) for skin temperature and 23.94 ± 0.08 h for activity (Table 1).

Chronotype estimates from skin temperature

Skin temperature onsets and offsets were estimated from diel temperature profiles using the four described approaches. In most cases onset and offset estimates by Broken Stick Regression, Structural Change and Changepoint analyses fell between the time of relatively constant temperature plateaus (close to the minimum or maximum, respectively) and the first change (increase or decrease) in slope. In contrast, BCPA_Tskin estimates were usually between steepest slope and the timepoint where temperature started flattening again within the chosen 12-h window (Fig. 4, Additional file 1: Table S1).

Fig. 4

Determination of onset and offset of diel change in skin temperature. Bird-specific onset (A) and offset (B) estimates are calculated by four methods (different colours) for eight Great Tits (female (F) and male (M)) kept in outdoor aviaries during winter. Mean (± SEM) chronotype estimates (in numeric daytime hours) and information on colour coding are given for each method in the inlay. Each panel represents skin temperatures of one Great Tit with mean (black dots) ± standard errors (grey bars) per 10-min bin, the fitted sinusoidal curve (blue, based on 3-harmonic smoothing), and the skin temperature chronotype estimates as vertical lines. Daylight phases are indicated in yellow (sunrise–sunset)

All methods differed in their estimates, with on average earliest onset times (mean ± SEM) of 01:38 ± 0:23 h from the Broken Stick Regression and latest onsets of 09:07 ± 0:19 h from the BCPA_Tskin. Onset estimates from the other two methods were in between and similar to each other, with Changepoint (03:06 ± 0:28 h) being overall 9.6 min later than Structural Change (02:56 ± 0:34 h). For the offset, BCPA_Tskin had again the latest time estimates of on average 19:30 ± 0:30 h, clearly distinct from the other three methods which were between 12:19 ± 0:47 h and 12:52 ± 0:21 h. Of those, the offsets of Structural Change were earliest, followed by the estimates from Changepoint (29 min later) and Broken Stick Regression (34 min later; for details see Additional file 1: Table S1 and Figure S6). The standard errors for onset and offset were similar to each other and ranged from 19 to 47 min.

Chronotype estimates from activity

Using BCPA, individual activity onsets and offsets were estimated from diel activity profiles. From visual inspection the estimates and the AIC values fitted the change of variation in the data reliably without arbitrary assignments (Additional file 1: Figure S7). Individuals became active on average at 07:30 h (SEM: 4.2 min) with low between individual variation (range: 07:10 to 07:50 h) and 56 min before mean sunrise (08:26 h). Similar, but reverse, patterns were found for activity offset which occurred on average at 16:15 h (SEM: 6.0 min), with a range from 16:00 to 16:40 h and 41 min after mean sunset (15:34 h) (Additional file 1: Table S2).

Comparison of chronotype estimates from skin temperature and activity

The relationship between chronotype from skin temperature and behavioural activity depended on the estimation method from skin temperature. Most methods showed a much earlier increase, and also earlier decrease, of skin temperature than of activity. Behavioural activity started on average 5.86 h (Broken Stick Regression), 4.56 h (Structural Change) and 4.4 h (Changepoint) after the skin temperature estimates for onset, but it was 1.62 h earlier than the BCPA_Tskin onset estimates. Similarly, activity ended on average 3.38 h (Broken Stick Regression), 3.94 h (Structural Change) and 3.46 h (Changepoint) after the skin temperature offset estimates, but 3.25 h before BCPA_Tskin offsets (Fig. 5). Variation in activity chronotype estimates was about four- to eightfold lower than found in estimates from skin temperature for onsets (SEM: 4.2 vs. 19.2–34.2 min) as well as for offsets (SEM: 6.0 vs. 21.6–46.8 min).

Fig. 5

Individual chronotype estimates based on individually averaged diel patterns of skin temperature and activity. Colours indicate different analytic methods for skin temperature, estimates for activity chronotype are shown as black diamonds. Great Tits (female (F) and male (M)) are ranked by activity onset. Daylight phases are indicated in yellow (sunrise–sunset)

Skin temperature estimates from Broken Stick Regression, Structural Change and Changepoint were strongly positively correlated with each other with r = 0.646–0.841 for onset and r = 0.719–0.906 for offset. In contrast, they were more weakly correlated with activity chronotypes (onset: r = − 0.131–0.138 and offset: r = − 0.394– −0.681). In contrast, BCPA_Tskin estimates showed weaker correlations with the other skin temperature estimates (onset: r = 0.099–0.195 and offset: r = 0.131–0.676) and different correlations with the activity onset (r = 0.612) and offset (r = − 0.040) than the other methods (Table 2).

Table 2 Correlation table between estimates of chronotype from different methods, for onset (below diagonal) and offset (above diagonal, and italicised)

Discussion

Minimally invasive, remotely collected data can provide great opportunities to study diel patterns of behaviour and physiology in free-ranging animals. However, these data also bring new challenges for data processing and analyses. Here, we processed and analysed activity and skin temperature records from externally attached radio transmitters to extract and compare estimates of physiological and behavioural chronotype.

Evaluation of methodology

Although in some cases relatively clean rhythmic records can be achieved [57, 67], data on skin temperature are prone to confounding effects. In particular, ambient temperature, wind [54, 56, 68] and solar radiation [55] can affect skin temperatures, and such confounds can be accentuated depending on the tightness between skin and transmitter (e.g., if the attachment loosened). Strong effects of solar radiation, depending on exposure, in particular can make skin temperature measurements from daytime harder to interpret than those from the night [55]. Therefore, for estimating chronotype, measures at night and in the early morning are particularly suitable, when body temperature increases in diurnal animals. In addition to confounding the skin measurements, ambient temperature can also physiologically alter the daily profiles of body temperature and activity. Experimental studies in small mammals have reported altered amplitudes and absolute levels, but divergent effects on the shape and phase (i.e., timing) of rhythmic profiles. For example, with colder ambient temperatures, (core) body temperature decreased and levels of activity increased, whereas the shape and phase of both rhythms did not change in two studies [69, 70]. In contrast, if ambient temperatures reached the individual’s physiological constraints (e.g., energetical depletion or overheating) and its internal state was affected, activity patterns were shifted in phase or altered in shape either asynchronised [71] or synchronised with the individual’s physiological rhythms [72].

In our study, skin temperature records showed high variability in noise across the data set. Attachment of transmitters was sometimes relatively loose, because we did not pluck feathers and also attached the transmitters to a piece of cloth to increase area of adhesion leading to partly noisy data and high variability in data quality of the skin temperature measurements (also discussed in [54]). Thus, prior to analyses we first assessed ambient temperature fluctuations. We found high variation between days and weak diel periodicity, which suggests that effects of ambient temperature on rhythmicity were small. Furthermore, all individuals were exposed to the same ambient temperatures, and diel profiles of skin temperature and activity were averaged across days accounting for between day variation of ambient temperature. Therefore, we focussed on skin temperature analyses without further consideration of ambient temperature. However, ambient temperature should always be recorded for further investigation, especially when individuals are exposed to different environmental temperatures.

Subsequently, we filtered our noisy time series to high quality levels, which greatly diminished the availability of data for analyses of chronotype. For these remaining data, however, diel fluctuations in skin temperature were well suited for analysis. Our time series were successfully filtered using rhythmicity parameters from periodicity and autocorrelation, and the derived diel skin temperature profiles showed the expected, roughly sinusoidal curves, with lower temperature during the inactive phase (e.g., [6, 8, 50]).

Extracted chronotype estimates from skin temperature (onset and offset) were correlated but differed between the four changepoint approaches because of the ways the changepoints were assigned, and because of the windows selected for analyses. In our study, the Broken Stick Regression and the BCPA_Tskin assigned only one changepoint. However, for different applications, settings can be adjusted and analysis windows reduced and/or shifted, to extract different changepoints of interest. In our study, there was high variation in skin temperature chronotypes so that one general analysis window did not fit all individuals equally well, especially for the BCPA_Tskin. In contrast, a standard window is suitable for traits with low variation within or between individuals, such as chronotypes from locomotion activity, where BCPA_Tskin selects distinct changes in variation and mean [29, 30]. Further important considerations for the selection of windows also include variation in environmental conditions, such as changes in light–dark cycle across the season. Issues relating to window selection can partly be solved by fitting models with multiple changepoints, and then selecting the one in the region of interest, as in our Structural Change approach.

One criterion for method selection, given our noisy time series, was relative independence from between individual differences in recorded absolute skin temperatures. Even though three of our methods (Broken Stick Regression, Structural Change and BCPA_TSkin) used the absolute temperature records for analyses, the algorithms estimated changepoints based on the input data set and, thus, for each individual independently. Our Changepoint method additionally accounted for this by focussing on temperature changes at the regions of interest (close to minimum or maximum) using a slope-based threshold rather than one absolute temperature threshold for all individuals which seems suitable for defining deep torpor [62]. Similarly, individual-based temperature thresholds have been proposed [39, 62, 73] and discussed [57] to account for individual variation in absolute temperature recordings in the face of methodological differences, e.g., in transmitter attachment.

In our study, variation in skin temperature profiles between individuals was the focal interest. Thus, high temporal resolution was crucial to detect patterns, but at the same time, data needed to be of high quality. If differences between individuals are large, as is the case for our chronotype estimates from skin temperature, they are easy to detect despite binning and smoothing. However, if temporal resolution is too coarse, and filtering is too strict, it can hamper the detection of low amounts of variation and obscure small differences [54, 70]. Thus, sampling interval and degree of filtering need to be carefully chosen to minimise potential loss of information.

Comparison of chronotype estimates of skin temperature and activity

Our data showed an increase in skin temperature anticipating the activity onset as shown in previous studies for birds and mammals, supporting the hypothesis of physiological preparation for diel activity [6, 52, 61, 74]. This rise occurred several hours before sunrise and activity onset according to most of the methods we used. As in humans and other species, skin temperatures of the birds we studied started to fall before the end of the active day. Thereby, most estimates of skin temperature offset occurred well before the end of activity, shortly after mid-day, and decreased thereafter [18]. The only exception was chronotype estimates based on BCPA_Tskin, for which onset occurred around sunrise and offset well after sunset. The reason was that BCPA_Tskin was fitting one distribution to the steep slope and another one to the plateau thereafter, whereas the three remaining methods aligned with skin temperature changes close to minimum or maximum.

In general, birds can show different magnitudes of temperature decrease during inactive phases ranging from shallow rest-phase hypothermia (few degrees) up to deep torpor (occasionally up to several tens degrees) depending on species, environmental conditions and the individuals’ states [75, 76]. Small passerines, such as the Great Tit, exhibit mainly the rest-phase hypothermia with decreases of max. 10 °C only [8, 37, 77]. The timing of hypothermia onset and offset depends on internal mechanisms and external cues associated with the light–dark cycle, ambient temperature and food availability [10, 78, 79]. In addition, locomotion activity itself can generate heat resulting in accelerated rewarming, and higher body temperatures during the active phase [7, 13, 52, 57]. In our study, the timing of temperature decrease must be carefully interpreted due to the potentially confounded skin temperature measurements during the active phase and effects of smoothing. Overall, skin temperature onsets (= offset of hypothermia, but onset of the physiological day) might be a more reliable tool to assign chronotypes than skin temperature offsets, for the reasons explained above.

Overall, we found large differences between individuals in skin temperature onset and offset, but the magnitude of these differences varied between analytic methods used. Furthermore, the chronotype estimates from most methods based on skin temperature differed considerably from estimates based on activity. In addition to differences in phase (timing) between skin temperature and activity explained above, estimates also differed in variation. Whereas activity chronotypes showed small variation as in previous studies [23, 29, 30], standard errors of the mean chronotypes (here, indicating variation between individuals) of skin temperature onset and offset were up to eightfold larger. This might partly be due to the chosen methodology, i.e., the chosen changepoint approach and the external measurement of skin temperature, but can also indicate large differences in the regulation of skin temperature and activity. One possibility could be differences in the adjustment of the two rhythms to the light–dark cycle. Especially activity is highly and directly affected by light cues, such as sunrise and sunset, light pulses and artificial light at night, either in the short- or long-term [10, 12, 47, 80, 81]. For example, various species first “light sample”, i.e., assessing the light levels after awakening but before leaving roosting sites [82]. In contrast, body temperature rhythms might rather reflect physiological and internal time, with reduced immediate responsiveness to the light–dark cycle and other environmental effects [17].

Closer investigation of day-to-day variation of skin temperature chronotypes could inform on such temporal flexibility within individuals and on their responsiveness to external cues. The methodology we have proposed can also be applied for day-to-day analyses. Thereby, onsets and offsets could be separately estimated for a series of days. These daily data could be used to disentangle consistent chronotypes from individual flexibility (i.e., calculate consistency and repeatability), and to quantify influences of specific environmental factors on skin temperature rhythmicity. However, this requires high-quality data and will bring additional challenges, such as the interpolation of data gaps, choice of selecting the analysis windows, and statistical analysis of potentially autocorrelated day-to-day estimates of onset and offset [83].

In our study, for the eight individuals with sufficient data, onset and offset of skin temperature and activity were only weakly and partly negatively correlated. Several other studies, usually under strictly controlled environmental conditions, have described both synchronisation and desynchronisation of body temperature and activity rhythms (reviewed in [6]). Often studies focus on changes in amplitudes and mean levels of diel rhythms, i.e., suppression or amplification (e.g., [39, 69, 84]). However, shifts in timing are important for understanding variation in chronotype, as well as responses of wildlife to anthropogenic changes [85]. External and internal factors can alter both rhythms, but effects and dynamics of change can differ for skin temperature versus activity, and results are heterogenous. Some studies found no change or synchronised shifts in body temperature and activity rhythms for stressors. Whereas both rhythms were unaffected by social stress situations in rats [39], synchronised shifts have been observed in free-running rats under constant dim light conditions [86] as well as in pigeons that were subject to manipulation of melatonin levels [87]. In contrast, food supplementation of resveratrol, a natural phenol from grapes, affected timing of body temperature and activity differently depending on age and the treatment duration in lemurs [88]. In cattle the two rhythms have been desynchronised with increased intraspecific competition shifting them in opposite directions [89]. Desynchronisation can also be forced by manipulation of the light–dark cycle either into constant free-running conditions [90] or into extended days of 28 h [91] resulting in shorter period lengths of body temperature than activity. Discrepant findings of different studies depend on the stressor and study model used, highlighting the complexity of the underlying mechanisms and functions of the two rhythms. Much of this still needs to be discovered in wildlife.

Conclusions

In this study we proposed methods to assess chronotype from skin temperature data. We developed standardised protocols to filter skin temperature data for rhythmicity, using periodicity and autocorrelation parameters. The protocols worked well and helped avoid arbitrary data selection but came at a cost of greatly reduced sample sizes. We then tested four different changepoint analyses for extracting chronotype estimates, i.e., the morning onset and evening offset of sinusoidal skin temperature curves. We found that three analyses (Broken Stick Regression, Structural Change and Changepoint) performed similarly well and yielded closely correlated, but somewhat different, estimates of chronotype. The fourth method, BCPA_Tskin, gave weakly correlated estimates of chronotype that were starkly differently phased. These discrepancies suggest that the methods are picking up different shape characteristics of skin temperature profiles. The choice of suitable methods might, therefore, depend on recorded time series curves, and further adjustments can be made to account for different species and study questions. Compared to chronotype derived from activity, estimates based on skin temperature were timed far earlier, were more variable, and magnified between individual differences. Estimates of chronotype based on activity and skin temperature were only weakly correlated. Through our protocols, we hope to stimulate use of the rapidly increasing time series data for investigating chronotype of behaviour and physiology, and for exploring potential discrepancies between the two processes that may indicate responses to environmental change. Future studies might include the comparison of day-to-day variation within individuals and overall between individual variation to disentangle temporal flexibility from consistent chronotypes. A main observation in our study are the distinct differences between skin temperature and activity in time series patterns as well as in chronotype and its variation. The activity chronotype is a useful tool to study behavioural and ecological timing, whereas body temperature chronotype might inform about individual timing of physiological processes broadly related to thermoregulation and metabolic rate. However, simultaneous investigation of behavioural and physiological rhythms has the advantage of exploring similarities and differences in proximate and ultimate processes in wildlife.