Study area and animals
We collected data registrations from 12 captive female moose, ages 2–19 years old, at the Kenai Moose Research Center, on the Kenai Peninsula, Alaska, USA (Fig. 1). The Kenai Moose Research Center is operated by the Alaska Department of Fish and Game and is located on the Kenai National Wildlife Refuge (60°43’ N, 150° 26’ W). The study area is characterized as boreal forest, including dense forests, wetlands, and open meadows. Present tree species are dominated by white spruce (Picea glauca), black spruce (Picea mariana), Alaska birch (Betula neoalaskana), Scouler’s willow (Salix scouleriana), and quaking aspen (Populus tremuloides). In addition to wild moose, the surrounding area is also inhabited by brown bears (Ursus arctos), black bears (Ursus americanus) and gray wolves (Canis lupus). The study animals were kept outdoors in two 2.6 km2 enclosures with free access to water and natural forage, which enabled the desired frequency of individual-specific fecal sample collection.
In May and July 2021, the captive female moose were immobilized with a combination of Thiafentanil oxalate (0.001–0.004 mg/kg estimated body mass; 10mg/mL; Wildlife Pharmaceuticals Inc., Windsor, CO, USA) and Xylazine (0.03–0.05 mg/kg estimated body mass; 100mg/mL; Lloyd Laboratories, Shenandoah, IA, USA) by intramuscular hand-injection and equipped with Vertex Plus Global Positioning System (GPS) collars (Vectronic Aerospace GmbH, Berlin, Germany) and Vaginal Implant Transmitters (VITs; size = big; Vectronic Aerospace GmbH; Berlin, Germany). Immobilization was reversed by a combination of Atipamezole HCl (0.005 mg/kg estimated body mass; ¼ dose intravenous, ¾ dose intramuscular; 5g/mL; Zoetis, Parsippany, NJ, USA) and intramuscularly administered Naltrexone HCl (100mg/mg Thiafentanil oxalate intramuscular; 50mg/mL; ZooPharm LLC, Laramie, WY, USA). GPS collars and VITs were manually removed from all animals without sedation on November 23rd, 2021.
Vaginal temperature and collar activity
The VITs transmitted vaginal temperature (Tv; °C) every 5 minutes to the GPS collar and recorded vaginal temperature every 17.10 minutes on board the VIT. During the aforementioned immobilization procedure, a VIT was inserted into the vaginal canal of each individual by using a lubricated (OB Lube; Jorgensen Laboratories Inc., Loveland; CO, USA), sterilized speculum (Sterile Disposable Vaginal Speculum; Jorgensen Laboratories Inc.) following previously established procedures (51). Tv was successfully downloaded from 11 individuals (from GPS-collar n = 10 and from the VIT n = 1). One VIT stopped recording on the 12th of August 2021 (i.e., before the start of fecal sample collection), and so this individual was excluded from further Tv analysis. Over 150,000 registrations of Tv were made during the dates with corresponding fecal progestagen values (i.e., August 23rd to October 15th, 2021). Daily mean Tv was calculated, resulting in 594 registrations during this period.
The GPS collars included a triaxial accelerometer which registered forward-backward, left-right, and up-down movements, stored as X, Y, and Z, respectively. Movement was recorded in 5-minute intervals as average values of each axis ranging between 0-255 at 6–8 Hz and stored in the collar. Overall activity was calculated by summing the X, Y and Z axes, giving values ranging from 0-765 every 5 minutes, with 0 representing no or low activity and 765 the highest activity level. Activity registrations were successfully stored for all 12 individuals resulting in over 185,000 registrations for each axis (X, Y, and Z) during the period of fecal sampling (i.e., August 23rd to October 15th, 2021). Daily summed activity during the fecal sampling period was calculated, resulting in 234 to 288 registrations a day.
VIT validation
We validated the accuracy and precision of the VIT temperature logger within the range of moose body temperature (i.e., within the range of vaginal temperature in this study and prior research (37)). We placed the VITs in a warm water bath heated to 37.00°C, and we then increased the temperature of the water bath by 0.50°C until 40.50°C. The water bath temperature was measured with a National Institute of Standards and Technology-certified thermometer (Traceable ® Digital Thermometer Model 90080-09, accuracy ± 0.05°C, 0 to 100°C, resolution 0.001°C (Webster, TX, USA)). We also tested the accuracy and precision of the VIT in an ice bath, which has been traditionally used to calibrate temperature logger VITs (52). We used linear mixed model regression, with individual VIT as a random effect, to determine accuracy of the loggers.
Fecal samples
A total of 468 fecal samples were collected from mid-August to mid-October 2021. The study period was selected based on previous literature on the timing of the species’ reproductive season, personal observations of the study population’s estrous behavior in previous years, and a desire to collect samples during the transition between anestrus to estrus (17, 21, 24). Individual samples were collected every other day from the 23rd of August to the 21st of September 2021. After this, samples were collected daily until the 16th of October 2021. All sampling occurred between 07:00–19:00 Alaska Daylight Time. Each sample weighed between 100-200g and was collected in pre-labeled Whirl-Pak bags (Nasco Whirl-Pak®, Fort Atkinson, WI, USA) during or shortly after defecation and stored on ice in a cooler bag, before being transferred to a -18°C propane freezer. The frozen samples were freeze dried (Labconco model 7752020, Kansas City, MO, USA) and homogenized to a powder to further distribute metabolites evenly before randomly subsampling 5g (53). The subsamples were shipped to an endocrine laboratory (Applied Biosciences, TX, USA) for further extraction and analysis. Progestagen concentration was determined by radioimmunoassay (RIA; Catalog #07-270102; ImmuChem Double Antibody, 125I RIA Kit, MP Biomedicals, Costa Mesa, CA, USA), following the manufacturer’s protocol. Fecal progestagen concentration is expressed as nanogram (ng) of immunoreactive fecal progestagen hormone metabolites per gram (g) dry fecal weight.
Environmental variables
Within the Kenai Moose Research Center there is a National Oceanic and Atmospheric Administration (NOAA), U.S. Climate Reference Network weather station which recorded environmental variables including ambient air temperature (Ta, °C), relative humidity (%), solar radiation (W/m2), windspeed (m/s), and precipitation (mm) every 5 minutes (54). Ambient air temperature and relative humidity recorded by the NOAA weather station were used to calculate dew point temperature (°C), which was further used to calculate actual vapor pressure (hPa) (55). The NOAA weather station had 1 day of missing data for ambient temperature, relative humidity, and solar radiation (25th of August, 2021). HOBO dataloggers (HOBO U23 Pro v2 Temperature/Relative Humidity datalogger; HOBO Pendant Temperature/Light 64K datalogger; Onset Computer Corp., Pocasset, MA, USA) were deployed adjacent to the NOAA weather station for the duration of the study. We used simple linear regression to estimate the missing NOAA values from the HOBO logger data for daily mean ambient air temperature (y = 0.96*Ta + 0.50, r2 = 0.99), daily range in ambient air temperature (0.88005*TaR + 0.49969, r2 = 0.98) relative humidity (y = 1.05859*RH -9.79112, r2 = 0.98), and solar radiation (y = .0004*lux – 0.1266, r2 = 0.98).
Data preparation and analysis
All data handling and statistical analysis was performed in R version 4.1.0 (R Core Team, 2021), using RStudio version 1.4.1717. Means are reported with ± SD in parenthesis. Tv and collar activity data were filtered for the 3 days following capture to exclude values which would likely be elevated due to lingering effects of the immobilizing drugs or capture-related stress (37). Since fecal collections occurred on alternate days for the first half of the sampling period, fecal progestagen values were linearly interpolated using the “approx” function from the ”zoo” package (56) in R to generate daily values to classify luteal activity. Similarly, daily values for Tv and collar activity were generated by calculating the daily mean Tv and the total daily sum of activity values (referred to as the daily sum of activity).
To determine a baseline level of fecal progestagen considered as indicative non-luteal activity, we calculated the mean of the 2 lowest values for each individual. Further, a threshold was created by multiplying the baseline level by 2, to differentiate between non-luteal activity and luteal activity. If fecal progestagen levels were below this threshold for at least 1 day, and then increased, and stayed above it for at least 14 days, we defined the day on which the threshold was exceeded as an onset of luteal activity (OLA), and the period during which fecal progestagen levels stayed above the threshold as luteal activity (Fig. 2). Because the duration of the luteal phase is currently unknown for moose we based our OLA classification method on existing approaches for classifying normal luteal activity in cattle (57), along with the average luteal phase length in other seasonal breeding ruminants such as marsh deer (Blastocerus dichotomus), domestic sheep (Ovis aries), and domestic goats (Capra hircus; 15, 14 and 16 days, respectively (34, 58, 59)). The date of the first OLA was assigned as day 0, which was used as a scale to explore variations in Tv and collar activity for that individual in relation to luteal activity. The above calculations and classifications were done by applying a rolling function using the “zoo” R package (56).
The four moose which we were not able to determine an OLA for either had erratic fecal progestagen profiles (age = 2 years old), fecal progestagen levels which did not stay above the set threshold for the minimum required days, or an erratic fecal progestagen profile which stayed above the set threshold but was not preceded by a day in which the level was below the set threshold (age = 19 years old). Thermal and activity registrations from these individuals were therefore excluded from further statistical analyses. These individuals were either below or above prime age (2.5–10.5 years old (21)), indicating that our method is more successful when applied on prime aged moose.
Relationship between luteal activity and vaginal temperature
To determine when changes in daily mean Tv occur in relation to the onset of luteal activity, and to account for non-linear relationships between response and explanatory variables, we used Generalized Additive Mixed Models (GAMMs; using the “mcgv” R package from Wood (60)). All explanatory variables included in the model were adjusted to daily values to match up with the daily frequency of fecal progestagen values. Daily mean Tv was selected as the response variable and OLA was included as the main explanatory variable to explore the relationship to luteal activity. Additionally, explanatory variables that are known to influence core body temperature in moose were included: daily sum of activity, Julian day (to account for seasonal variation in Tv (36)), daily mean Ta (°C), daily range of Ta (°C), daily mean vapor pressure (hPa), daily mean wind speed (m/s), total daily precipitation (mm), and total daily solar radiation (W/m2) (37). Correlation and structure of the explanatory variables were checked using the function “ggpairs” from the “GGally” (61) extension to the “ggplot2” R package (62). If the Pearson’s correlation coefficient between OLA and another explanatory variable was above 0.50 or below − 0.50, it was considered a high correlation (63) and this variable was not included together with OLA in the same model (this was the case for Julian day, vapor pressure and daily ambient temperature, see Additionalfile_1, additional files). Also, if two explanatory variables other than OLA had a correlation with each other above or below this threshold, the combination of these variables was not included in the same model (this was the case for daily ambient temperature and solar radiation, and wind speed and relative humidity, see Additionalfile_2, additional files).
These steps left us with the following explanatory variables: OLA, daily range of Ta, daily mean relative humidity, total daily precipitation, daily mean wind speed, total daily solar radiation, and daily sum of activity (Table 1, additional files). Additionally, to test if the trend in Tv was attributed to seasonal variation rather than OLA, we included Julian day in one of the models. We included the individual moose (“CollarID”) as a random intercept in all models to control for inter-individual variability and repeated measurements. A gaussian distribution with an identity link function and the maximum likelihood was used. We then ran a model selection based on Akaike’s Information Criterion, adjusted for small sample sizes (AICc, function “AICctab” from the “bbmle” R package (64)), by selecting the most parsimonious model within ΔAICc ≤ 2 (65). Residuals were assessed and basis dimensions for the parameter k were checked before we inspected diagnostic plots to validate model assumptions. Finally, the selected model was refitted with restricted maximum likelihood.
Table 1
Summary of the response and explanatory variables included in the model selection process.
Variable
|
Category
|
Name
|
Type
|
Definition
|
Range
|
Response
|
Tv
|
Continuous
|
Daily mean vaginal temperature (°C)
|
37.53–38.21
|
Explanatory
|
OLA
|
Discrete
|
Number of days in relation to the onset of luteal activity
|
-21-20
|
Explanatory
|
yday
|
Discrete
|
Julian day
|
237–288
|
Explanatory
|
Ta range
|
Continuous
|
Daily range of ambient temperature (°C)
|
1.80–20.00
|
Explanatory
|
Wind
|
Continuous
|
Daily mean wind speed (m/s)
|
0.40–2.80
|
Explanatory
|
RH
|
Continuous
|
Daily mean relative humidity (%)
|
37.40–94.70
|
Explanatory
|
Prec
|
Continuous
|
Total daily precipitation (mm)
|
0.00-17.80
|
Explanatory
|
Solar_rad
|
Continuous
|
Total daily solar radiation (W/m2)
|
0.51–18.23
|
Explanatory
|
Activity
|
Continuous
|
Daily sum of activity
|
2623.00-17821.00
|
Description: The table includes vaginal temperature, activity registrations and days in relation to onset of luteal activity (OLA) from captive female moose (n = 7), with corresponding environmental data. M/s = meters per second. Mm = millimeters. W/m2 = watts per square meter.
Relationship between luteal activity and collar activity
GAMMs were also used to determine when changes in daily sum of activity occur in relation to the onset of luteal activity. Daily sum of activity was selected as the response variable and OLA was included as the main explanatory variable to explore the relationship between luteal activity and activity values. The same environmental variables mentioned previously were included as they have been shown to affect behavior and habitat selection in moose (66). In addition, we included daily mean Tv as an explanatory variable since thermoregulation influences activity levels in moose (66), and the same random intercept (individual moose, “CollarID”). We used a gamma distribution with a log link function and applied the same model selection process as described above.