Reproductive state alters vocal characteristics of female North American red squirrels (Tamiasciurus hudsonicus)

Female advertisement of reproductive state and receptivity has the potential to play a large role in the mating systems of many taxa, but investigations of this phenomenon are underrepresented in the literature. North American red squirrels (Tamiasciurus hudsonicus) are highly territorial and engage in scramble competition mating, with males converging from spatially disparate territories to engage in mating chases. Given the narrow estrus window exhibited in this species, the ubiquitous use of vocalizations to advertise territory ownership, and the high synchronicity of males arriving from distant territories, we hypothesized that female vocalizations contain cues relating to their estrous state.  To test this hypothesis, we examined the spectral and temporal properties of female territorial rattle vocalizations collected from females of known reproductive condition over 3 years. While we found no distinct changes associated with estrus specifically, we did identify significant changes in the spectral characteristics of rattles relating to both female body mass and reproductive state relative to parturition. To the best of our knowledge, this is the first evidence of changes in vocal characteristics associated with late pregnancy in a non-human mammal.

Methods
Methods listed below are from:

Hare, McAdam, Dantzer, Lane, Boutin, and Newman (2024). Reproductive state alters vocal characteristics of female North American Red Squirrels (Tamiasciurus hudsonicus). Journal of Mammalogy (In Press).

"Study system
We studied a wild population of T. hudsonicus associated with the Kluane Red Squirrel Project (KRSP) from 2018 – 2020 in southwestern Yukon (61˚N, 138˚W) within the Champagne and Aishihik First Nations’ traditional territory (Dantzer et al. 2020). Audio recordings were collected from 108 individual female squirrels at 6 separate study sites spaced 5 km apart along the Alaska Highway. Each individual in the study was tagged with unique alphanumeric ear tags in each ear (Monel #1 National Band & Tag Co., Newport, KY, USA). We also threaded colored wires through ear tags for visual identification from a distance (McAdam et al. 2007). As part of the KRSP, all individuals were routinely trapped and handled to establish territory ownership, assess body mass and reproductive condition, and determine matrilineal relationships (McAdam et al. 2007). All data were collected following guidelines in Sikes (2016). 

Audio recording
Continuous audio recordings were taken in 3-day bouts using Zoom H2N digital audio recorders (Zoom North America, NY, USA) placed at the center of each female midden. In brief, recorders were fixed to trees at the center of each midden at the beginning of each recording session using bungee cables at ~1.5m height prior to squirrel emergence at sunrise (or 5 am when close to the solstice). Batteries and SD cards were changed every 24 hours. Stereo recordings were taken at a 44.1 kHz sampling rate in a 16-bit WAV format then mixed down to mono prior to analysis.  

Audio analysis
Audio recordings were selected for inclusion if they took place within ± 10 days of a live trapping event in which individual mass was recorded in order to control for the effects of body size on the pitch of vocalizations (Bradbury and Vehrencamp 1998). This range was selected to include only recordings from squirrels for which we had relatively high certainty of reproductive state through palpation and maintain the most proximate mass measurements given the 2-week periodicity of regular trapping that is part of the KRSP trapping protocol. Selected recordings were analyzed with Kaleidoscope Pro bioacoustics software (Wildlife Acoustics, Maynard, MA, USA) using the cluster analysis function, which differentiates red squirrel rattles from other signals based on an existing template, as used previously (Siracusa et al. 2019). This process identified 2.23 million acoustic signals that could be red squirrel rattles from 9.85 TB of uncompressed continuous recordings constituting 8,214 total days of sampling. Signals were then run through the ‘Seewave’ package (Sueur et al. 2008) in R version 4.1.0 (R Core Team 2022) to extract a range of spectral and temporal values from each respective clip (sensu Siracusa et al. 2019).

To differentiate calls made by territory owners from those produced by neighbouring red squirrels, all clips identified through Kaleidoscope were run through a Linear Discriminant Function analysis, built using previously published acoustic parameters (Siracusa et al. 2019) through the ‘caret’ package in R (Kuhn 2008) using a ground tested dataset of 10,650 manually identified acoustic signals (Siracusa et al. 2019). In short, an observer recorded continuous audio while remaining stationary on the territory of 49 focal squirrels, recording the time whenever a rattle vocalization was produced by the focal individual. Acoustic signals extracted from these recordings were then labelled as either an owner rattle, neighbour rattle, or not a rattle vocalization (Supplementary Data SD4). Applying this linear discriminant function reduced the sample to 11,488 potential owner rattles at a certainty of 95% for manual validation.

Potential owner rattles were reviewed manually through both visual examination of the spectrograph and aural examination of the audio clip. Audio clips that had been detected as owner rattles were selected for inclusion if they were actually a red squirrel rattle, had a minimum of a 3:1 signal to noise ratio, had no spectral or temporal clipping of the vocalization, and there were no competing signals within the vocalizations from conspecific, heterospecific, or anthropogenic sources. This manual validation produced 897 good quality owner rattles for use in our analysis for which we had high certainty of individual identity. It is worth noting, however, that these data were collected autonomously—we did not witness the squirrel vocalization at the time the recording were made.

Estimating estrus date
As we were interested in the timing of acoustic changes relative to estrus dates, owner rattles were filtered further to include clips from individuals from which we could accurately estimate estrus date.  This was done using known parturition dates from live-trapping and monitoring associated with ongoing research at the KRSP (McAdam et al. 2007), then backtracking 35 days, which is within the average reported gestation period for a North American Red Squirrel (Dolbeer 1973; Lair 1985; Steele 1998).  From this point, a possible window in which a mating chase could occur was defined as ± 5 days, allowing for individual variation in gestation length, which has been documented as occurring between 31 to 35 days (Lair 1985; Steele 1998; Dantzer et al. 2011). Calls that occurred within ± 100 days of parturition were included to capture the full extent of changes that would occur across a reproductive cycle—a range that was selected to encompass meaningful reproductive changes across the reproductive season, omitting samples collected in the winter from a single individual. This filtering criterion reduced our sample size from 897 rattles to a final dataset of 120 owner rattles from 22 individuals in which we were able to capture their estrus window while being confident of signaller identity.

Using known parturition dates also allowed us to assign different reproductive states to callers with fairly high certainty. In brief, any calls included that were over 35 days to parturition were considered as coming from non-reproductive squirrels, while any calls occurring within 35 days were considered as coming from pregnant squirrels (Steele 1998; McAdam et al. 2007). Finally, any calls that were collected within 71 days following parturition were classified as coming from lactating squirrels that had not yet weaned their pups (McAdam et al. 2007; Dantzer et al. 2011). Each of the 120 owner rattles included were assigned to one of these 3 reproductive states for further analysis of changes across reproductive condition (Fig. 1).

Statistical analyses
Each acoustic parameter of interest was first visually examined for normality and homoscedasticity before being analyzed via a generalized additive model (GAM) to test for differences across days to parturition using the ‘mgcv’ package in R (Wood, 2017). We selected GAMs to test our hypothesis regarding estrus because the additive nature of these models would allow for any pronounced changes occurring within our defined estrus window to be detected, while linear models would likely miss these discrete changes. Smoothing parameters were used on continuous variables using the restricted maximum likelihood (REML) method on days to parturition at recording and scaled individual mass at most proximate trapping. Julian date of recording was omitted because it displayed high concurvity with days to parturition across all variables of interest. Year of recording and field site were coded as factorial covariates. Individual age in years and age in years squared were included as continuous integer coefficients to capture any linear and quadratic effects of age. Finally, individual squirrel identification was included as a random intercept to control for individual variation. This was particularly important for our models because there was uneven sampling of individuals across reproductive states, and the inclusion of a random intercept for individual identity accommodated for the effects of repeated sampling of certain squirrels (Fig. 1). Response variables were grouped into either spectral or temporal variables, with spectral variables including mean frequency (Fmean), fundamental frequency (F0), and Shannon entropy (entropy)—while temporal variables included call duration (duration), rattle trill rate (rate), and rattle duty cycle (duty cycle; for definitions see Supplementary Data SD3). Trends for each response variable were examined visually using a GAM spline with nine knots along days to parturition, specifically looking for distinct changes within the ± 5-day estrus window. Nine knots were selected for visual inspection as this was the number of smoothing terms produced through the REML smoothing method. 

Terms that exhibited statistically significant effects in our GAMs were further examined through post-hoc generalized linear mixed models (GLMM) to determine whether effects were statistically significant. Days to parturition was replaced with the factorial covariate reproductive condition as calculated relative to parturition in order to test for significant changes across reproductive state. The effective degrees of freedom (edf) calculated through GAMs were used to inform GLMM construction, including quadratic terms for mass when it was indicated to have a non-linear relationship with response variables. All statistical tests were run on R version 4.2.1 (R Core Team 2022) and for all tests, statistically significant trends were inferred at a level of significance of  = 0.05 following visual examination of confidence intervals. Due to the number of tests run, a Holm-Bonferroni correction for multiple comparisons was applied to any obtained P-values, then compared at our level of significance. For any factorial variables, parameter estimates ± standard errors are reported along with relevant test statistics and associated P-values are reported in-text when significant and otherwise can be found in the supplementary data (see Supplementary Data SD6, SD8, SD10, SD12, SD14, SD16 – SD20)."