Study site, subjects and dates of recordings

Calls of hinds and calves were recorded in June 2011 and 2012 at the farm of the University of Castilla-La Mancha (Albacete, Spain, 38°57’10”N, 1°47’00”W, 690 m a.s.l.). The population originated in 1994 from 15 male and 50 female pure Iberian deer from a nearby Las Dehesas public game reserve in Alpera (Albacete) and from Cabaneros National Park (Toledo). The animals used in this study were born and kept in four 10,000 m² enclosures on an irrigated pasture. They were fed ad libitum with a diet of barley straw and meal from barley, alfalfa, oats and sugar beets [60].

All mothers and calves were captive-born and kept together in permanent groups (4 groups in 2011 and 3 groups in 2012) separately from adult stags and yearlings. The groups ranged in size from 6 hinds and 2 calves to 30 hinds and 24 calves (mean ± SD = 17.4 ± 7.6 hinds and 15.0 ± 7.7 calves per group). The entire population of animals from which we collected calls counted 61 hinds and 52 calves in 2011 and 61 hinds and 53 calves in 2012 (45 hinds were the same in both years). The age of mothers was 2–19 years (mean ± SD = 10.0 ± 4.6 years). The calves were born from 6 May to 23 June in 2011 and from 13 May to 14 June in 2012. During data collection, the age of calves varied of 1–52 days. All study calves were singletons, excluding two sibling pairs, one in 2011 and one in 2012.

All hinds and calves were individually labeled with Allflex (Palmerston North, New Zealand) plastic ear tags and all hinds were additionally marked with the Allflex color collars with numbers. Three times a day, the farm staff inspected each enclosure with hinds and calves for newborns and labeled them individually with an ear tag, which allowed us to establish a mother for each calf.

Call and body mass collection

For acoustic recordings of hind and calf (48 kHz, 16 bit), we used a solid state recorder Marantz PMD-660 (D&M Professional, Kanagawa, Japan) with a Sennheiser K6-ME66 cardioid electret condenser microphone (Sennheiser electronic, Wedemark, Germany). The distance from the hand-held microphone to the animals was 5–35 m, the level of recording was adjusted during the recordings accordingly to the intensity of the produced calls.

We recorded calls daily, from 6:00–7:00 to 12:00–13:00, often with synchronous video, using a digital camcorder Panasonic HDC-HS100 (Panasonic Corp., Kadoma, Japan). During recordings, individual identities of callers producing calls through the mouth and through the nose were labeled by voice. Recordings have been conducted both inside and outside the outdoor enclosures during different contexts: everyday routine activity, when mothers searched for their offspring which were hidden in the enclosures; at translocations to small paddocks, from where the animals were taken for weightings; during temporal separations of hinds and calves after the weighing and at short separations by presence of researchers in enclosure on the way between a mother and her calf. In total, in 2011 and 2012, we collected 30 hours of audio recordings (16 hours in 2011 and 14 hours in 2012) from 28 individual hinds (9 hinds were recorded in both 2011 and 2012) and from 31 individual calves. All animals were weighted one time with Mettler-Toledo ID1 scales (Mettler-Toledo S.A.E., Barcelona, Spain) as the part of routine farm management [61] during the periods of acoustic recordings. The age of calves at weightings ranged from 1 to 40 days.

Acoustic analyses

For acoustic analyses, we took only calls of good quality with high signal-to-noise ratios sufficient for analysis of all acoustic variables, measured in this study that were not disrupted by wind or the calls of other animals or overloaded during the recordings. We analysed only individually identified calls of known call type (nasal or oral). Calls were classified to nasal and oral call types based on voice comments of researchers made during recording, by video clips, made synchronously with the recordings, by the obvious nasal quality of sound within a recording (Additional files 1 and 2) and by the difference in call energy distribution, shifted towards higher frequencies in oral calls due to the shortening the vocal tract at opening of the mouth. These methods of classification to nasal and oral call types were previously applied for the Iberian red deer [37], for goitred gazelles [31,32] and for saiga antelopes [36]. Two researchers (OS and IV) independently classified all calls, and we took for analysis only calls where both researchers were concordant in their judgments concerning their type. To reduce pseudoreplication, we took calls from different recording sessions per animal and from different parts within session, because calls from the same sequence are commonly more similar in their acoustic structure than calls from different sequences [62]. The mean ± SD number of sessions per animal was 5.3 ± 4.4, and we took up to 15 high-quality calls per individual per call type (nasal and oral) from 28 mothers and 31 young (16 males, 15 females). We analysed 801 calls of mothers (354 oral and 447 nasal) and 469 calls of calves (281 oral and 188 nasal), for 1270 calls in total.

Acoustic analyses were conducted in the same way for hinds and calves and for the oral and nasal calls. For each nasal and oral call, we measured the same set of 14 acoustic variables: 2 temporal, 6 variables of fundamental frequency (f0) and 6 power variables, as they proved their use for estimating vocal individual identity in red deer [28,47,48]. Before analysis, the calls were downsampled to 11025 Hz and high-pass filtered at 50 Hz, to increase frequency resolution and to reduce the low-frequency background noise. We measured the duration of each call and the duration from call onset to the point of maximum f0 (dur-to-max) manually on the screen with the reticule cursor in the spectrogram window (Hamming window, FFT 1024 points, frame 50% and overlap 96.87%) by using Avisoft SASLab Pro software (Avisoft Bioacoustics, Berlin, Germany). Then we performed manual measurements on the screen with the standard marker cursor of the start (f0beg), maximum (f0max) and end (f0end) fundamental frequencies of each call (Figure 1). Measurements were exported automatically to Microsoft Excel (Microsoft Corp., Redmond, WA, USA).

In a 50 ms call fragment symmetrical about f0 maximum, we created the power spectrum, from which we automatically measured fpeak, representing the value of the frequency of maximum amplitude and the q25, q50 and q75, representing the lower, medium and upper quartiles, covering 25%, 50% and 75% of the energy of the call spectrum respectively (Figure 1). On the same spectrum, we estimated (in dB) the power-f0, representing the relative power of the f0 band compared to the peak harmonic, on the screen using two harmonic cursors (Figure 1). The power-f0 was equal to 0 when the f0 band coincided with the fpeak band. In addition, we recorded the peak-harm, representing the order number of the harmonic with the maximum energy.

In addition, we measured the f0 variables following Reby and McComb [63] by using the Praat DSP package [64]. The f0 contour was extracted by using a cross-correlation algorithm (to Pitch (cc) command in Praat). The time steps in the analysis were 0.01 s for hinds and 0.005 s for calves; the lower and upper limits of the f0 range were 50–400 Hz for hinds and 100–1200 Hz for calves. A preliminary visual analysis of the spectrograms in Avisoft showed that the lower limit was lower than the minimum f0 for calls of either hinds or calves. Spurious values and octave jumps in the f0 contour were corrected manually on the basis of the spectrograms. Values of f0min, f0max, the depth of frequency modulation f0 (∆f0 = f0max - f0min) and average f0 of a call (f0mean) were taken automatically by using by using the Pitch info command in the Pitch edit window.

Two different methods of measuring f0max (one using Avisoft and another using Praat) applied to the same calls, resulted in very similar values. Coefficients of correlation, calculated separately for the oral and for the nasal calls of hinds and calves, ranged between 0.994 and 0.996 (0.988 < R² < 0.992). Thus, for subsequent acoustic analyses we could select between these methods and we used the f0 values measured with Praat.

We did not measure formants, as they can only be measured either in low-frequency calls with closely spaced harmonics or in noisy calls (e.g. stag harsh roars), where the sound energy is dispersed over the call spectrum [26,37,55]. In most calls of hinds and in all calls of calves, formant frequencies could not be measured because they fell either between the harmonics and were thus invisible, or they were indistinguishable from the harmonics because they coincided with them in frequency.

Call samples and statistical analyses

Statistical analyses were conducted using STATISTICA v. 6.0 (StatSoft, Tulsa, OK, USA) and R v.3.0.1 [65]. Means are given as mean ± SD, all tests were two-tailed, and differences were considered significant whenever p < 0.05. Distributions of most measured parameter values (excepting fpeak and peak-harm) did not depart from normality and distributions of all mean parameter values did not depart from normality (Kolmogorov-Smirnov test, p > 0.05). As parametric ANOVA and discriminant function analysis (DFA) are relatively robust to departures from normality [66], this was not an obstacle to the application of these tests.

We used a repeated measures ANOVA controlled for individuality, to compare the parameter values between oral and nasal calls. We used 1–15 (9.3 ± 5.2) calls per animal per call type, from 28 hinds (300 oral and 325 nasal calls) and from 31 calves (16 males, 15 females; 281 oral and 188 nasal calls). From hinds which provided calls in both years, we took calls recorded only in one of the years (2011 or 2012). For oral and nasal calls of these 28 hinds and 31 calves, we calculated average values of acoustic variables.

We used DFA to calculate the probability of the assignment of calls to the correct individual for three call samples in both hinds and calves (of nasal calls, oral calls and the pooled sample of oral and nasal calls). We took 5–10 (9.0 ± 1.5) calls per animal per call type from 22 hinds (209 oral and 212 nasal calls in total) and from 17 calves (10 males, 7 females; 149 oral and 134 nasal calls in total). DFA requires balanced sample sizes per group, and thus we excluded from the analyses all the animals with less than 5 oral/nasal calls. We took all the calls from all the animals with from 5 to 10 oral/nasal calls and randomly selected 10 calls per type from animals with more than 10 measured calls of each type. We included 11 of the 14 measured call variables in the DFA, excluding fpeak and peak-harm (for not meeting the criterion of normality), and f0min (because it was used for calculating another variable).

Then we investigated the stability of acoustic individuality of hind oral and nasal calls between years for hinds that provided calls in both years. We classified hind calls from 2012 year with DFA functions derived from 2011, considering the value of the correct cross-validation as a measure of the retention of individuality over time [39,43,44,67]. We used 7–15 (12.2 ± 3.4) oral calls per animal per year from 5 hinds, which provided a sufficient number of oral calls in both years (in total 65 oral calls in 2011 and 57 oral calls in 2012). Also, we used 10–15 (13.7 ± 1.8) nasal calls per animal from 9 hinds which provided sufficient numbers of nasal calls in both years (in total 125 nasal calls in 2011 and 121 nasal calls in 2012).

We used Wilks’ Lambda values to estimate how strongly acoustic variables of calls contribute to the discrimination of individuals. With a 2×2 Yates’ chi-squared test, we compared the values of correct assignment of nasal and oral calls to the correct individual. We used a repeated measures ANOVA to compare the percentages of correct assignment of oral and nasal calls to particular individuals. To validate our DFA results, we calculated the random values of correct assignment of calls to individual by applying randomization procedure with macros, created in R. The random values were averaged from DFAs performed on 1000 randomized permutations on the data sets as described by [68]. For example, to calculate the random value of classifying oral calls to individual hinds, each permutation procedure included the random permutation of 209 calls among 22 randomization groups, respectively to 22 individual hinds which were examined, and followed by DFA standard procedure built-in in STATISTICA. All other permutation procedures were made similarly. Using a distribution obtained by the permutations, we noted whether the observed value exceeded 95%, 99% or 99.9% of the values within the distribution [68]. If the observed value exceeded 95%, 99% or 99.9% of values within this distribution, we established that the observed value did differ significantly from the random one with a probability p < 0.05, p < 0.01 or p < 0.001 respectively [44,68,69].

Because body mass should theoretically be proportional to the cube of a linear dimension like body size, for hinds we used log₃ body mass to calculate Pearson’s correlation between body mass and acoustic variables of oral and nasal calls. For calves we used ANCOVA (with sex as a fixed categorical factor and body mass as continuous factor) for estimating effects of sex, body mass and the conjoint effect of sex*body mass on acoustic variables of oral and nasal calls. We used one-way ANOVA to compare body mass values between male and female calves.