\documentclass[10pt]{article}

\usepackage{fullpage}
\usepackage{setspace}
\usepackage{parskip}
\usepackage{titlesec}
\usepackage[section]{placeins}
\usepackage{xcolor}
\usepackage{breakcites}
\usepackage{lineno}
\usepackage{hyphenat}





\PassOptionsToPackage{hyphens}{url}
\usepackage[colorlinks = true,
            linkcolor = blue,
            urlcolor  = blue,
            citecolor = blue,
            anchorcolor = blue]{hyperref}
\usepackage{etoolbox}
\makeatletter
\patchcmd\@combinedblfloats{\box\@outputbox}{\unvbox\@outputbox}{}{%
  \errmessage{\noexpand\@combinedblfloats could not be patched}%
}%
\makeatother


\usepackage{natbib}




\renewenvironment{abstract}
  {{\bfseries\noindent{\abstractname}\par\nobreak}\footnotesize}
  {\bigskip}

\titlespacing{\section}{0pt}{*3}{*1}
\titlespacing{\subsection}{0pt}{*2}{*0.5}
\titlespacing{\subsubsection}{0pt}{*1.5}{0pt}


\usepackage{authblk}


\usepackage{graphicx}
\usepackage[space]{grffile}
\usepackage{latexsym}
\usepackage{textcomp}
\usepackage{longtable}
\usepackage{tabulary}
\usepackage{booktabs,array,multirow}
\usepackage{amsfonts,amsmath,amssymb}
\providecommand\citet{\cite}
\providecommand\citep{\cite}
\providecommand\citealt{\cite}
% You can conditionalize code for latexml or normal latex using this.
\newif\iflatexml\latexmlfalse
\providecommand{\tightlist}{\setlength{\itemsep}{0pt}\setlength{\parskip}{0pt}}%

\AtBeginDocument{\DeclareGraphicsExtensions{.pdf,.PDF,.eps,.EPS,.png,.PNG,.tif,.TIF,.jpg,.JPG,.jpeg,.JPEG}}

\usepackage[utf8]{inputenc}
\usepackage[ngerman,greek,english]{babel}



\usepackage{float}








\begin{document}

\title{Altitude, latitude, and climate zone as determinants of mountain hare
(\emph{Lepus timidus}) coat colour change}



\author[1]{Allan Stokes}%
\author[2]{Tim Hofmeester}%
\author[3]{Neri Thorsen}%
\author[4]{John Odden}%
\author[1]{John Linnell}%
\author[1]{Simen Pedersen}%
\affil[1]{Inland Norway University of Applied Sciences - Evenstad Campus}%
\affil[2]{Sveriges Lantbruksuniversitet Fakulteten for Skogsvetenskap}%
\affil[3]{Norwegian Institute for Nature Research Oslo}%
\affil[4]{Norwegian Institute for Natural Research}%


\vspace{-1em}



  
  \date{May 16, 2024}


\begingroup
\let\center\flushleft
\let\endcenter\endflushleft
\maketitle
\endgroup





\selectlanguage{english}
\begin{abstract}
Local adaptation to annually changing environments has evolved in
numerous species. Seasonal coat colour change is an adaptation that has
evolved in multiple mammal and bird species occupying areas that
experience seasonal snow cover. It has a critical impact on fitness as
predation risk may increase when an individual is mismatched against its
habitat's background colour. In this paper we investigate the impact of
landscape covariates on moult timing in a native winter-adapted
herbivore, the mountain hare (\emph{Lepus timidus}), throughout Norway.
Data was collected between 2011 and 2019 at 678 camera trap locations
deployed across an environmental gradient. Based on this data, we
created a Bayesian multinomial logistic regression model that quantified
the correlations between landscape covariates and coat colour phenology
and analysed among season and year moult timing variation. Our results
demonstrate that mountain hare moult timing is strongly correlated with
altitude and latitude with hares that live at higher latitudes and
altitudes keeping their winter white coats for longer than their
conspecifics that inhabit lower latitudes and altitudes. Moult timing
was also weakly correlated with climate zone with hares that live in
coastal climates keeping their winter white coats for longer than hares
that live in continental climates. We found evidence of some among year
moult timing variation in spring, but not in autumn. We conclude that
mountain hare moult timing has adapted to local environmental conditions
throughout Norway.%
\end{abstract}%



\sloppy


\emph{Introduction}

For species living in seasonal environments (e.g. from summer to winter
or from dry to wet season) local adaptations to annually changing
environmental conditions may evolve. Numerous species have evolved to
time their life history events to match these changes in local seasonal
conditions (Bradshaw and Holzapfel 2007, Williams et al. 2015). To time
these phenological events, one reliable ``zeitgeber'' is daylength or
photoperiod (Gwinner 2003, Hofman 2004), either on its own or in
combination with other variables, such as temperature (Watson 1963,
Jackes and Watson 1975, Larkin et al. 2001) and snow cover (Watson 1963,
Flux 1970, Zimova et al. 2014).

Many animal species use photoperiod to time breeding (Goldman 1991,
Gwinner 1996b, Dawson et al. 2001, Coppack and Pulido 2004), moulting
(Lesher and Kendeigh 1941, Lyman 1943, Bissonnette and Bailey 1944), and
migration (Gwinner 1996a), and other life history events. As photoperiod
remains constant between years at specific locations, between year
variation in local conditions could result in photoperiod timed
phenological events being mistimed against the local environment.
Fluctuations in environmental variables, such as precipitation
(Villellas et al. 2014) and temperature (Ashmore and Janzen 2003,
Kreyling et al. 2019), can result in increased within-species phenotypic
variation in a variety of plant and animals species when compared to
individuals of the same species that live in more stable habitats.
Consequently, synchrony in phenological timing of individuals within a
population is expected to increase with climate stability.

Animals occupying areas that are seasonally covered by snow live in
environments that change from dark in summer to white in winter. As a
predator avoidance strategy, at least 21 species (Mills et al. 2018,
Zimova et al. 2018) have adapted seasonal changes in colouration of fur
and feather, which provides camouflage in both a winter white and summer
dark landscapes (Wallace 1879, Cott 1940, Merilaita and Lind 2005). To
provide optimal camouflage, the timing of coat colour change should be
synchronised with the period of continuous snow cover. Individuals would
thus be expected to adapt the timing of their coat colour change to
local conditions. Indeed, mismatched timing of coat colour change is
linked to range contractions and population declines in several species
including snowshoe hares, (Diefenbach et al. 2016, Sultaire et al.
2016), mountain hares (Acevedo et al. 2012, Pedersen et al. 2017), rock
ptarmigan (Imperio et al. 2013), and white-tailed ptarmigan (Wang et al.
2002), showing the importance of correct timing. However, snow
conditions might not be stable from year to year, and there might be
seasonal differences in the predictability of the appearance and
disappearance of snow.

There are some suggestions that snow cover might be more stable in
spring compared to autumn. Snowshoe hares (\emph{Lepus americanus} )
(Zimova et al. 2014) and least weasels (\emph{Mustela nivalis nivalis} )
(Atmeh et al. 2018) exhibit limited phenotypic variation in moult timing
in parts of their distribution during the spring moult, when
transitioning from white to brown, but not during the autumn moult, when
transitioning from brown to white. Therefore, the timing of seasonal
coat colour change is expected to be more synchronised in spring
compared to autumn. This could result in clearer differences in the
timing of among years moult timing in the spring compared to the autumn.
However, this is the first long-term study over a geographical area
large enough to test these predictions.

Mountain hares (\emph{Lepus timidus} ) express seasonal coat colour
change in most of their range, except the subspecies of Irish hare
(\emph{L. t. hibernicus} ) found in Ireland (Mills et al. 2018). They
are a generalist herbivore inhabiting boreal and alpine areas that
occupy a wide range of climatic, latitudinal, and altitudinal gradients,
experiencing large variations in winter snow cover duration. Coastal
areas in the south and south-west of Norway experience relatively short
snow cover duration compared to inland areas and areas in the north
(Schuler et al. 2006) with coastal areas in the south and south-west
receiving as little as one month of snow cover per year (Tallaksen et
al. 2018). Additionally, coastal areas experience greater between-year
variation in the depth and extent of snow cover than inland areas
(www.senorge.no).

Here, we provide the first quantitative study of variation in mountain
hare moult timing with nine years of data using 678 camera locations
along an extensive geographic gradient in Norway. We studied 1) how the
timing of moult varied with local geographical conditions and 2) how the
timing of moult varied among years and seasons, utilising camera trap
data in a Bayesian multinomial logistic regression model framework.
Snowshoe hares that live at lower altitudes and latitudes displayed
winter coats for a longer time compared to their low latitude (Grange
1932, but see Zimova et al. 2019) and altitude conspecifics (Holmgren et
al. 2001, Nowak et al. 2020, Zimova et al. 2020b). Also, increased snow
cover in continental areas is likely to result in hares living in these
areas keeping their winter coats for longer than hares in coastal areas.
Therefore, we used altitude, latitude, and climatic zone, distinguishing
between coastal and continental climates, as indicators of local
geographical conditions for our first aim. For our second aim, we
predicted that moult timing would be more synchronised among individuals
in spring compared to autumn, based on previous studies on snowshoe
hares (Mills et al. 2013, Zimova et al. 2014). Furthermore, we predict
larger among year variation in the timing of moult in spring compared to
autumn.

\emph{Methods}

\emph{Data collection}

Mountain hares predominately inhabit boreal and alpine habitats. They
have a circumpolar distribution spread across Europe and Asia from
Britain and Ireland in the west to Japan in the east (Angerbjorn and
Flux 1995). They are native to Norway and are found throughout the
country, with their distribution ranging from sea level to around 1,600
m, which is above the tree line. Population size within Scandinavia is
thought to be decreasing (Pedersen and Pedersen. 2012, Elmhagen et al.
2015, Pedersen et al. 2017) resulting in the species being classified as
``near threatened'' in the 2015 Norwegian Red List (Henriksen and Hilmo.
2015).

We utilized images from camera traps that were deployed by the SCANDCAM
project (viltkamera.nina.no), originally designed to monitor Eurasian
lynx (\emph{Lynx lynx} ) (Hofmeester et al. 2021). Camera traps were
deployed in multiple study areas in an extensive grid with approximately
one camera per 50 km\textsuperscript{2} grid cell (Figure 1). For this
study, we selected all mountain hare records from the period between
25th November 2010 and 25th June 2019. Images containing mountain hares
were recorded at 678 locations from 2010-2019 spanning a latitudinal
gradient from 58\selectlanguage{ngerman}°N to 69° N and altitudes from 0 to 841 m a.s.l. in
Norway (Figure 1). To reduce pseudo-replication, we discarded
observations recorded within 60 minutes of the previous observation.
Mountain hares and invasive European hares (\emph{Lepus europaeus} ),
which are also present in south eastern Norway (Viken County), were
differentiated using the species descriptions contained in Smith et al
(2018). When mountain hares were identified, we classified moulting
stage in three categories (Zimova et al. 2020a) estimating the
proportion of the hares' coat (excluding the long white belly) which was
white. 1) Hares with \selectlanguage{english}[?]90\% white fur were classified as ``white'', 2)
hares with \selectlanguage{english}[?]10\% white fur were classified as ``brown'', and 3) all
other hares were classified as ``moulting''. All images were accessed on
viltkamera.nina.no and were visually classified by one of two observers,
and quality controlled whenever the observer was uncertain of how to
classify the images.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image1/image1}

\end{center}
\end{figure}

Figure 1. The location of camera trap sites that recorded a mountain
hare observation during 2010-2019 (black circles). The camera traps were
deployed across an environmental gradient within Norway by the SCANDCAM
project, with latitudes varying between 58\selectlanguage{ngerman}°N and 69° N.

\emph{Covariates}

We divided the year into two seasons, ``spring'' and ``autumn''. Spring
was defined as ordinal day 1 (1\textsuperscript{st} January) to ordinal
day 212 (1\textsuperscript{st} August in non-leap years), and autumn for
the rest of the year. Ordinal day 212 was chosen as all mountain hares
had moulted to their summer brown coats by this date and had not started
moulting back to winter white. Altitude and latitude were extracted
based on camera trap positions. We obtained altitude from a digital
elevation model (DEM) with 50 m\textsuperscript{2} resolution
(Kartkatalogen 2007) (Suppl. 2) using the \emph{raster} package's
(Hijmans 2022) extract function. We obtained climate zone as a
continuous variable with a resolution of 1 km\textsuperscript{2} from
Bakkestuen et al (2008) (Suppl. 1). We converted climate zone vector
data to a raster using the \emph{fasterize} (version 1.0.4) package
(Ross 2020). Bakkestuen et al (2008) mapped climate zones by conducting
principal component analysis (PCA) using terrain data, climatic data,
hydrological data, and geological data. A positive PCA value indicates a
continental climate whereas a negative value indicates a coastal
climate.

\emph{Data analysis}

We used multinomial logistic regression to estimate the probability of a
hare being in moult category \emph{i} (white, moulting, or brown) on
each ordinal day \emph{d} . We included year \emph{k} specific
intercepts to test if moult timing varied between years. The model used
the following equation:

\begin{equation}
p\left(y=i\right)=\frac{e^{\alpha_{\text{ik}}+\beta_{1i}\times d+\beta_{2i}\times a_{j}{+\beta}_{3i}\times l_{j}{+\beta}_{4i}\times c_{j}+s_{\text{ij}}}}{1+\sum_{m=1}^{i-1}e^{\alpha_{i}+\beta_{1i}\times d+\beta_{2i}\times a_{j}{+\beta}_{3i}\times l_{j}{+\beta}_{4i}\times c_{j}+s_{\text{ij}}}}\nonumber \\
 \end{equation}

Where \emph{d} is ordinal day, and \emph{a\textsubscript{j}}
,\emph{l\textsubscript{j}} , and \emph{c\textsubscript{j}} ,
respectively, are the altitude, latitude, and climate zone at camera
site \emph{j} (as adapted from Zimova et al 2019 and 2020b).
Parameters\selectlanguage{greek}\emph{β\textsubscript{1-4j}}\selectlanguage{english} represent the slopes for the
different covariates. We included a category \emph{i} and year \emph{k}
specific intercept \selectlanguage{greek}\emph{α\textsubscript{ik}}\selectlanguage{english} as well as a site-specific
random intercept \emph{s\textsubscript{ij}} to correct for multiple
observations per camera site. The brown category was set to 0 in both
spring and autumn models to provide a baseline for comparison with
moulting and white categories.

We implemented the multinomial logistic regression models in a Bayesian
framework using JAGS (Denwood 2016) called from R V4.1.3 (R Core Team
2022) with the \emph{jagsUI} package (version 1.5.2) (Kellner 2021). We
standardized altitude, latitude, and climate zone (mean = 0, SD = 1)
before running the models. Ordinal day was included as an explanatory
covariate to enable estimation of the probability of a hare being white,
brown, or moulting between coats on specific days of the year. We
checked for collinearity between the covariates using both the variance
inflation factor (VIF) and Pearson correlation coefficient. For every
covariate combination the Pearson values were below 0.6 (Suppl 6) and
VIF values were below 2.0 (Suppl 7).

We used uninformative, normally distributed priors with a mean of 0 and
precision of 0.01 for all slopes and the year \emph{k} specific
intercepts. For the site \emph{s} random intercept, we used a mean of 0
and a standard deviation defined as a vague prior with a uniform
distribution between 0 and 100. We ran the model with three chains all
thinned by 100 for 120,000 iterations, with a burn-in of 60,000
iterations. We confirmed model convergence using traceplots and the
Gelman-Rubin convergence statistic (R-hat) (Brooks and Gelman 1998) with
all variables used in the final models having R-hat values of 1.10 or
less. Additional models in which the dataset was subset to only include
cameras located south of 61\selectlanguage{ngerman}° North were used to test if model
performance was affected by camera trap placement north of this
latitude. The results were consistent with those obtained using the full
dataset indicating that camera placement north of 61° did not affect
model performance. We produced all figures using the
\emph{ggplot2}(Wickham 2016), \emph{raster} (Hijmans 2022), and
\emph{cowplot} (Wilke 2020) packages.

\emph{Results}

Between 2010 and 2019 a total of 9,981 mountain hare observations were
obtained at 678 camera trap locations (Figure 1 and Suppl 4). Only two
observations were recorded in 2010, which we discarded due to the small
sample size. Of the remaining observations, 7,454 were recorded in
``spring'' and 2,525 were recorded in ``autumn''. The number of camera
traps deployed across Norway increased throughout the study period,
leading to an increase in the number of observations obtained in each
year (see Suppl 3 for location of traps that recorded observations in
each year). Camera traps that recorded mountain hare observations were
located between altitudes of 0 and 841 m.a.s.l., latitudes were between
58° N and 69° N, and climate zone PCA values were between -2.68 (coastal
climate) and 2.86 (continental climate).

The results from our multinomial logistic regression model indicate that
all three explanatory covariates, altitude, latitude, and climate zone,
correlated with moult timing in spring and autumn. We found strong
evidence for hares keeping their winter white coats for a longer
duration at increased altitude and latitude (Figures 2.a and b.). The
probability of being white increased with altitude in spring
(\selectlanguage{greek}\emph{β\textsubscript{2,white}}\selectlanguage{english} = 1.492, 95\% credible interval (CI) =
1.309 to 1.673; Figure 2.a) and in autumn
(\selectlanguage{greek}\emph{β\textsubscript{2,white}}\selectlanguage{english} = 0.870, 95\% CI = 0.639 to 1.117;
Figure 2.a). Similarly, the probability of being white increased with
latitude in spring (\selectlanguage{greek}\emph{β\textsubscript{3,white}}\selectlanguage{english} = 2.472, 95\% CI =
2.265 to 2.672; Figure 2.b) and in autumn
(\selectlanguage{greek}\emph{β\textsubscript{3,white}}\selectlanguage{english} = 0.642, 95\% CI = 0.443 to 0.844;
Figure 2.b). Climate zone had the smallest effect size in both spring
and autumn. The probability of being white decreased with an
increasingly inland climate (increasing PCS values) in both spring
(\selectlanguage{greek}\emph{β\textsubscript{4,white}}\selectlanguage{english} = - 0.335, 95\% CI -0.510 to - 0.163;
Figure 2.c) and autumn (\selectlanguage{greek}\emph{β\textsubscript{4,white}}\selectlanguage{english} = - 0.287, 95\%
CI - 0.525 to - 0.053; Figure 2.c).

The ordinal day where the probability of being white is 50 \% occurred
27 days earlier in spring and 23 days later in autumn for hares at sea
level compared to 700 m a.s.l. (Figure 2 a). At 59\selectlanguage{ngerman}° N the ordinal day
where the probability of being white is 50\% occurred 27 days earlier in
spring and from brown to white 14 days later in autumn compared to hares
at 65° N (Figure 2.b). Hares living in inland climates moulted 9 days
earlier in spring and 11 days later in autumn (Figure 2. c).

Mountain hare moult timing varied across mainland Norway (Figure 3).
Starting in mid-winter moving into spring, there was a travelling wave
of moulting moving from lower to higher altitudes and latitudes (Suppl
10). The opposite effect was observed in autumn. The altitudinal and
latitudinal moult timing gradients show that hares that inhabit
mountainous areas and the north of Norway kept their winter white coats
for longer than their conspecifics that inhabit low altitude and
latitude areas. The model output, which represents all of Norway,
indicates that mountain hares in coastal areas moult later in spring and
earlier in autumn (Figure 2.c.). The prediction maps show that hares in
southern Norway's coastal areas Norway moult earlier in spring and later
in autumn (Figure 3.a and 3.f.), indicating that the effect of latitude
is stronger than the effect of climate zone. This is consistent with the
climate zone effect sizes being smaller than those of altitude and
latitude. Predictions for areas outside of camera trap locations (Figure
1) should be interpreted with care as these results are extrapolated.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image2/image2}

\end{center}
\end{figure}

Figure 2. Timing of spring and autumn moult at different altitudes (a),
latitudes (b), and climate zones (c). Solid lines represent the
predicted probability of being white and shaded areas indicate the 95\%
credible intervals. Figures use the mean of all intercepts included in
the final models.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image3/image3}

\end{center}
\end{figure}

Figure 3. Prediction maps with a resolution of 1 x 1
km\textsuperscript{2} giving the probability of mountain hares being
white across Norway on ordinal days 105, 135, 165, 255, 285, and 315.
The probability of being white was predicted using the model output and
the environmental covariates contained in every cell.

Timing of spring and autumn moults varied slightly between years (Figure
4). From 2013 onwards the 95\% CIs are consistently narrower in spring
than in autumn, which probably results from increased sample sizes in
spring (Suppl 4). For the spring intercepts, multiple years have
non-overlapping CIs for 90\% white, which signifies the start of spring
moulting, and 90\% brown, which signifies the end of spring moulting.
For example, 2013 does not overlap with 2014 and 2017 does not overlap
with 2018. For the autumn intercepts, the brown CIs, signifying the
start of autumn moulting, overlapped in most years. The CIs for 90\%
white, signifying the end of autumn moulting, overlapped in all years.
The time taken for moulting in spring to finish, represented by the
number of days between 90\% of hares being white and 90\% of hares being
brown, ranged between 46 days in 2019 and 64 days in 2013 and 2018
(Suppl 9). The time taken for moulting to occur in autumn ranged between
38 days in 2011 and 72 days 2013. In autumn, both the shortest and
longest time taken were in years with limited sample sizes (Suppl 4).
Therefore, these results should be viewed with caution.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image4/image4}

\end{center}
\end{figure}

Figure 4. Between year intercepts with 95\% CI for spring and autumn at
mean latitude, altitude, and climatic zone. Grey dots with 95\% CI
represent the mean ordinal day on which 90\% of hares are white. Brown
dots with 95\% CI represent the mean ordinal day on which 90\% of hares
are brown. The supplemental information contains the sample sizes (Suppl
4), intercept beta values (Suppl 5), and the number of days between the
midpoints (Suppl 9).

\emph{Discussion}

We used camera trap by-catch observations to provide the first
quantitative assessment of how large-scale environmental variables
correlate with mountain hare moult timing. As predicted, hares at higher
altitudes and latitudes moulted later in spring and earlier in autumn,
keeping their winter white coats for longer, compared to their lowland
and low altitude conspecifics. Contrary to our prediction, we found
support for a later spring moult and earlier autumn moult in areas
characterised by coastal climates, rather than inland climates,
resulting in hares that live in areas with coastal climates keeping
their winter white coats for longer than their inland conspecifics.

The unexpected correlation between climate zone and moult timing (Figure
2.c.) could result from variation explained by climate zone also being
partially explained by altitude (Pearson correlation coefficient
\textgreater{} 0.5 in both seasons) and latitude (Pearson correlation
coefficient \textless{} -0.3 in both seasons) (Suppl 6). Additionally,
all camera traps from 63\selectlanguage{ngerman}° North are close to the coast which may
confound results. However, in southern Norway there is a coastal to
continental moult timing gradient shown in the prediction maps (Figure
3.a. and 3.f.). Furthermore, it is likely that there will be increased
variability in moult timing in coastal climates as there is increased
among year variation in snow cover and duration compared to continental
climates (senorge.no).

Both altitude (Figure 2.a) and latitude (Figure 2.b) had the predicted
correlation with moult timing. Mountain hares that inhabit areas
characterised by higher altitude and latitude moulted later in the
spring and earlier in the autumn, keeping their winter white coats
longer, likely due to an expectation of longer snow cover duration.
These results are consistent with previous studies of other mountain
hare populations and other hare species, which found that increased
elevation correlated with mountain hares keeping their winter coats for
longer (Watson 1963) and increased latitude correlated with snowshoe
hares keeping their winter white coats for longer (Grange 1932).
However, a more recent study on snowshoe hares found no evidence of
moult timing variation across a latitudinal gradient in autumn and, in
spring, hares that live at high latitudes became brown earlier than
those at lower latitudes (Zimova et al. 2019). The correlations between
moult timing and altitude and latitude indicates that mountain hare
populations have adapted to local conditions, which suggests that gene
flow between populations is insufficient to dilute local adaptations.
This is particularly true for altitude as this variable can change
significantly over a short geographic distance.

The non-overlapping CIs between some years in the spring moult (Figure
4) indicate that there is some between year phenotypic variation.
However, the between year difference in moult timing is small, which is
consistent with photoperiod rather than climate being the main driver of
mountain hare moult timing (reviewed in Zimova et al. 2018). This is
consistent with similar studies conducted on snowshoe hares (Mills et
al. 2013, Zimova et al. 2014) and least weasels (Atmeh et al. 2018),
which found evidence of between year phenotypic variation in spring, but
not in autumn. The limited moult timing variation may reduce fitness as
an inability to change moulting patterns in response to among year
variation in snow extent and duration will increase camouflage mismatch
and, consequently, decrease survival probability (Zimova et al. 2016).
Behavioural responses to mismatch are possible, micro-habitat patch
selection, changed diurnal activity patterns, or manually removing
winter fur during the spring moult. However, there is currently no
evidence of snowshoe hares modifying their behaviour in response to
being mismatched (Zimova et al. 2014).

Moult timing variation could result from population level phenotypic
variation, individual level phenotypic plasticity, or a combination of
both factors. As our methodology did not facilitate monitoring specific
individuals, we could not disentangle the relative importance of
population and individual level variation. Observations were obtained at
specific locations in multiple years (Suppl 8). This makes it likely
that some individuals were recorded in multiple years, increasing the
probability of individual phenotypic plasticity influencing results.

As the camera traps were deployed to maximise the probability of
detecting lynx their locations within the surveyed area were not
randomised and some areas of Norway were not monitored. Additionally, we
obtained three times as many observations in spring compared to autumn
(Suppl 4), which is probably caused by increased hare activity patterns
during the mating season (Pettigrew et al. 2021) and the spring dataset
containing observations for 58 more days. This resulted in the intercept
CIs (Figure 4) being larger in autumn than in spring. Despite these
issues, the prediction maps (Figure 3 and Suppl 10) indicate that the
model produces accurate results.

Understanding the role that environmental characteristics have on moult
timing is vital when assessing the impact that climate change may have
on species that express a seasonal coat colour change. Our study sites
span 1,300 km and 11 latitudinal degrees, from sea level to 841 m.a.s.l.
making this is the first study that investigates moult phenology over a
large, continuous climatic gradient spanning three biomes (temperate
forest, boreal forest, and alpine tundra). Analysing the correlation
between the explanatory variables and mountain hare moult timing will
enable us to predict how the species will react to climate change
induced reductions in snow cover extent and duration. We plan to
investigate this issue in a future paper.

\emph{Acknowledgements}

We thank Håkon Bjørnstad and Lina Lange for classifying camera trap
images. The project was supported by internal funding from Inland Norway
University of Applied Sciences, Swedish Agricultural university, and
Norwegian Institute for Nature Research. We also thank Morten Odden for
support during the project's initial stage and the numerous volunteers
that help to maintain camera traps.

\emph{Data Accessibility}

Dataset used for analysis including phenotypic and environmental data
and the location where the observations were collected: Figshare, DOI:
10.6084/m9.figshare.22560340

\emph{}

\emph{References}

Acevedo, P., A. Jiménez-Valverde, J. Melo-Ferreira, R. Real, and P. C.
Alves. 2012. Parapatric species and the implications for climate change
studies: a case study on hares in Europe. Global Change
Biology\textbf{18} :1509-1519.

Ashmore, G. M., and F. J. Janzen. 2003. Phenotypic variation in smooth
softshell turtles (\emph{Apalone mutica} ) from eggs incubated in
constant versus fluctuating temperatures. Oecologia\textbf{134}
:182-188.

Atmeh, K., A. Andruszkiewicz, and K. Zub. 2018. Climate change is
affecting mortality of weasels due to camouflage mismatch. Scientific
Reports \textbf{8} :7648.

Bakkestuen, V., L. Erikstad, and R. Halvorsen. 2008. Step-less models
for regional environmental variation in Norway. Journal of
Biogeography\textbf{35} :1906-1922.

Bissonnette, T. H., and E. E. Bailey. 1944. Experimental modification
and control of molts and changes of coat-color in weasels by controlled
lighting. Annals of the New York Academy of Sciences\textbf{45}
:221-260.

Bradshaw, W. E., and C. M. Holzapfel. 2007. Evolution of animal
photoperiodism. Annual Review of Ecology, Evolution, and
Systematics\textbf{38} :1-25.

Brooks, S. P., and A. Gelman. 1998. General methods for monitoring
convergence of iterative simulations. Journal of Computational and
Graphical Statistics \textbf{7} :434-455.

Coppack, T., and F. Pulido. 2004. Photoperiodic response and the
adaptability of avian life cycles to environmental change. Pages
131-150\emph{in} A. Moller, W. Fiedler, and P. Berthold, editors. Birds
and Climate Change.

Cott, H. B. 1940. Adaptive coloration in animals. Methuen, London.

Dawson, A., V. M. King, G. E. Bentley, and G. F. Ball. 2001.
Photoperiodic control of seasonality in birds. Journal of Biological
Rhythms \textbf{16} :365-380.

Denwood, M. 2016. runjags: an R package providing interface utilities,
model templates, parallel computing methods and additional distributions
for MCMC models in JAGS. Journal of Statistical Software\textbf{71}
:1-25.

Diefenbach, D. R., S. L. Rathbun, J. K. Vreeland, D. Grove, and W. J.
Kanapaux. 2016. Evidence for range contraction of snowshoe hare in
Pennsylvania. Northeastern Naturalist \textbf{23} :229-248.

Elmhagen, B., J. Kindberg, P. Hellstrom, and A. Angerbjorn. 2015. A
boreal invasion in response to climate change? Range shifts and
community effects in the borderland between forest and tundra.
Ambio\textbf{44 Suppl 1} :S39-50.

Flux, J. E. C. 1970. Colour change of mountain hares (\emph{Lepus
timidus scoticus} ) in north-east Scotland. Journal of
Zoology\textbf{162} :345-358.

Goldman, B. D. 1991. Mammalian photoperiodic system: formal properties
and neuroendocrine mechanisms of photoperiodic time measurement. Journal
of Biological Rhythms \textbf{16} :283-301.

Grange, W. B. 1932. The pelages and color changes of the snowshoe
hare,\emph{Lepus americanus phaeonotus} , Allen. Journal of
Mammalogy\textbf{13} :99-116.

Gwinner, E. 1996a. Circadian and circannual programmes in avian
migration. The Journal of Experimental Biology \textbf{199} :39-48.

Gwinner, E. 1996b. Circannual clocks in avian reproduction and
migration. IBIS \textbf{138} :47-63.

Gwinner, E. 2003. Circannual rhythms in birds. Current Opinion in
Neurobiology \textbf{13} :770-778.

Henriksen, S., and O. Hilmo. 2015. The Norwegian red list for species
2015. Norwegian Biodiversity Information Centre, Trondheim, Norway.
http://www.biodiversity.no/Pages/135380

Hijmans, R. J. 2022. raster: Geographic Data Analysis and Modeling.

Hofman, M. A. 2004. The brain's calendar: neural mechanisms of seasonal
timing. Biological reviews of the Cambridge Philosophical
Society\textbf{79} :61-77.

Hofmeester, T. R., N. H. Thorsen, J. P. G. M. Cromsigt, J. Kindberg, H.
Andren, J. D. C. Linnell, and J. Odden. 2021. Effects of camera\selectlanguage{english}-trap
placement and number on detection of members of a mammalian assemblage.
Ecosphere \textbf{12} :e03662.

Holmgren, N., P. Jonsson, and L. Wennerberg. 2001. Geographical
variation in the timing of breeding and moult in dunlin \emph{Calidris
alpina} on the Palearctic tundra. Polar Biology \textbf{24} :369-377.

Imperio, S., R. Bionda, R. Viterbi, and A. Provenzale. 2013. Climate
change and human disturbance can lead to local extinction of Alpine rock
ptarmigan: new insight from the western Italian Alps. PLoS One\textbf{8}
:e81598.

Jackes, A. D., and A. Watson. 1975. Winter whitening of Scottish
Mountain hares (\emph{Lepus timidus scoficus} ) in relation to
daylength, temperature and snow-lie. Journal of Zoology \textbf{176}
:403--409.

Kartkatalogen. 2007.

Kellner, K. 2021. jagsUI: a wrapper around 'rjags' to streamline 'JAGS'
analyses.

Kreyling, J., S. J. Puechmaille, A. V. Malyshev, and F. Valladares.
2019. Phenotypic plasticity closely linked to climate at origin and
resulting in increased mortality under warming and frost stress in a
common grass. Ecology and Evolution \textbf{9} :1344-1352.

Larkin, J. E., D. A. Freeman, and I. Zucker. 2001. Low ambient
temperature accelerates short-day responses in Siberian hamsters by
altering responsiveness to melatonin. Journal of Biological
Rhythms\textbf{16} :76-86.

Lesher, S. W., and S. C. Kendeigh. 1941. Effect of photoperiod on
molting of feathers. Wilson Ornithological Society \textbf{53} :169-180.

Lyman, C. P. 1943. Control of coat color in the varying hare \emph{Lepus
americanus erxleben} . Bulletin of the Museum of Comparative Zoology at
Harvard College \textbf{93} :393--466.

Merilaita, S., and J. Lind. 2005. Background-matching and disruptive
coloration, and the evolution of cryptic coloration. Proceedings of the
Royal Society B: Biological Sciences \textbf{272} :665-670.

Mills, L. S., E. V. Bragina, A. V. Kumar, M. Zimova, D. J. R. Lafferty,
J. Feltner, B. M. Davis, K. Hacklander, P. C. Alves, J. M. Good, J.
Melo-Ferreira, A. Dietz, A. V. Abramov, N. Lopatina, and K. Fay. 2018.
Winter color polymorphisms identify global hot spots for evolutionary
rescue from climate change. Science \textbf{359} :1033-1036.

Mills, L. S., M. Zimova, J. Oyler, S. Running, J. T. Abatzoglou, and P.
M. Lukacs. 2013. Camouflage mismatch in seasonal coat color due to
decreased snow duration. Proceedings of the National Academy of Sciences
of the United States of America \textbf{110} :7360-7365.

Nowak, K., J. Berger, A. Panikowski, D. G. Reid, A. L. Jacob, G. Newman,
N. E. Young, J. P. Beckmann, and S. A. Richards. 2020. Using community
photography to investigate phenology: a case study of coat molt in the
mountain goat (\emph{Oreamnos americanus} ) with missing data. Ecology
and Evolution \textbf{10} :13488-13499.

Pedersen, S., M. Odden, and H. C. Pedersen. 2017. Climate change induced
molting mismatch? Mountain hare abundance reduced by duration of snow
cover and predator abundance. Ecosphere \textbf{8} :e01722.

Pedersen, S., and H. C. Pedersen. 2012. The population status of
mountain hare in Norway-state of knowledge. Norwegian Institute for
Nature Research, Trondheim, Norway.

Pettigrew, G. W., V. Di Vita, M. Pettigrew, and J. S. Gilchrist. 2021.
The diel activity pattern of mountain hare (\emph{Lepus timidus} ) on
managed heather moorland in Scotland. Ecology and Evolution\textbf{11}
:7106-7113.

R Core Team. 2022. R: A language and environment for statistical
computing.

Ross, N. 2020. fasterize: fast polygon to raster conversion.

Schuler, D. V., S. Beldring, E. J. Forland, L. A. Roald, and T. E.
Skaugen. 2006. Snow cover and snow water equivalent in Norway - current
conditions (1961--1990) and scenarios for the future (2071-2100).

Smith, A. T., C. H. Johnston, P. C. Alves, and K. Hacklander. 2018.
Lagomorphs: pikas, rabbits, and hares of the world. Johns Hopkins
University Press.

Sultaire, S. M., J. N. Pauli, K. J. Martin, M. W. Meyer, M. Notaro, and
B. Zuckerberg. 2016. Climate change surpasses land-use change in the
contracting range boundary of a winter-adapted mammal. Proceedings of
the Royal Society B: Biological Sciences \textbf{283} :20153104.

Tallaksen, L. M., K. Gisnas, J. H. Stagge, I. B. Nilsen, and J. Rizzi.
2018. Five decades of warming: impacts on snow cover in Norway.
Hydrology Research \textbf{49} :670-688.

Villellas, J., R. Berjano, A. Terrab, and M. B. Garcia. 2014. Divergence
between phenotypic and genetic variation within populations of a common
herb across Europe. Ecosphere \textbf{5} :Article 56.

Wallace, A. R. 1879. The protective colours of animals.

Wang, G., N. T. Hobbs, H. Galbraith, and K. M. Giesen. 2002. Signatures
of large-scale and local climates on the demography of white-tailed
ptarmigan in Rocky Mountain National Park, Colorado, USA. Int J
Biometeorol \textbf{46} :197-201.

Watson, A. 1963. The effect of climate on the colour changes of mountain
hares in Scotland. Proceedings of the Zoological Society of
London\textbf{141} :823-835.

Wickham, H. 2016. ggplot2: elegant graphics for data analysis.
Springer-Verlag New York.

Wilke, C. O. 2020. cowplot: streamlined plot theme and plot annotations
for 'ggplot2'.

Williams, C. M., H. A. Henry, and B. J. Sinclair. 2015. Cold truths: how
winter drives responses of terrestrial organisms to climate change.
Biological Reviews \textbf{90} :214-235.

Zimova, M., L. S. Barnard, B. M. Davis, A. V. Kumar, D. J. R. Lafferty,
and L. S. Mills. 2020a. Using remote cameras to measure seasonal molts.
Ecosphere \textbf{11} :e03084.

Zimova, M., S. T. Giery, S. Newey, J. J. Nowak, M. Spencer, and L. S.
Mills. 2020b. Lack of phenological shift leads to increased camouflage
mismatch in mountain hares. Proceedings of the Royal Society B:
Biological Sciences \textbf{287} :20201786.

Zimova, M., K. Hacklander, J. M. Good, J. Melo-Ferreira, P. C. Alves,
and L. S. Mills. 2018. Function and underlying mechanisms of seasonal
colour moulting in mammals and birds: what keeps them changing in a
warming world? Biological Reviews \textbf{93} :1478-1498.

Zimova, M., L. S. Mills, P. M. Lukacs, and M. S. Mitchell. 2014.
Snowshoe hares display limited phenotypic plasticity to mismatch in
seasonal camouflage. Proceedings of the Royal Society B: Biological
Sciences \textbf{281} :20140029.

Zimova, M., L. S. Mills, and J. J. Nowak. 2016. High fitness costs of
climate change-induced camouflage mismatch. Ecology Letters\textbf{19}
:299-307.

Zimova, M., A. P. K. Siren, J. J. Nowak, A. M. Bryan, J. S. Ivan, T. L.
Morelli, S. L. Suhrer, J. Whittington, L. S. Mills, and P. Morellato.
2019. Local climate determines vulnerability to camouflage mismatch in
snowshoe hares. Global Ecology and Biogeography \textbf{29} :503-515.

Supplemental information\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image5/image5}

\end{center}
\end{figure}

Supp 1. Climatic zones across Norway obtained from Bakkestuen et al
(2008). Negative values indicate an oceanic climate whereas positive
values indicate a continental climate. Values were obtained via
principle component analysis (PCA).\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image6/image6}

\end{center}
\end{figure}

Supp 2. Digital elevation model (Kartkatalogen 2007) with a resolution
of 50 m\textsuperscript{2} giving the altitude across Norway. The
altitude of the camera traps used in this study varies between 0 m and
841 m above sea level.\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image7/image7}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image8/image8}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image9/image9}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image10/image10}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image11/image11}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image12/image12}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image13/image13}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image14/image14}

\end{center}
\end{figure}\selectlanguage{english}
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.70\columnwidth]{figures/image15/image15}

\end{center}
\end{figure}

Supp 3. Location of the camera traps that recorded at least one
observation in each year. The locations are laid over the DEM. Suppl 4.
The number of samples recorded by camera traps in spring and autumn from
2010 to 2019. The number of samples recorded increased as the study
progressed due to an increase in the number of traps deployed.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Year & Spring & Autumn\tabularnewline
\midrule
\endhead
2010 & NA & 2\tabularnewline
2011 & 40 & 36\tabularnewline
2012 & 113 & 80\tabularnewline
2013 & 258 & 23\tabularnewline
2014 & 308 & 165\tabularnewline
2015 & 859 & 204\tabularnewline
2016 & 602 & 358\tabularnewline
2017 & 1,517 & 644\tabularnewline
2018 & 2,035 & 1,027\tabularnewline
2019 & 1,738 & NA\tabularnewline
\bottomrule
\end{longtable}

Suppl 5. The beta intercepts taken from our Bayesian multinomial
logistic regression model. They illustrate the between season and year
variation in moult timing.\selectlanguage{english}
\begin{longtable}[]{@{}lll@{}}
\toprule
Year & Spring (95\% CI) & Autumn (95\% CI)\tabularnewline
\midrule
\endhead
2011 & 20.932 (15.809 to 29.746) & -44.525 (-48.253 to
-40.601)\tabularnewline
2012 & 21.438 (19.727 to 23.146) & -43.931 (-47.367 to
-40.340)\tabularnewline
2013 & 27.432 (25.804 to 29.129) & -43.677 (-47.693 to
-39.776)\tabularnewline
2014 & 23.819 (22.364 to 25.148) & -43.201 (-46.692 to
-39.672)\tabularnewline
2015 & 23.923 (22.530 to 25.338) & -43.494 (-46.964 to
-39.944)\tabularnewline
2016 & 24.338 (22.925 to 25.727) & -44.456 (-47.803 to
-40.869)\tabularnewline
2017 & 22.761 (21.478 to 24.114) & -44.174 (-47.477 to
-40.675)\tabularnewline
2018 & 26.066 (24.622 to 27.529) & -44.424 (-47.751 to
-40.946)\tabularnewline
2019 & 24.942 (23.549 to 26.324) & NA\tabularnewline
\bottomrule
\end{longtable}

Suppl 6. Pearson's correlation coefficient (3 SF) between altitude,
latitude and climate zone in spring and autumn. All years and camera
sites were used for the calculations. The datasets were subset so that
only one observation per site was used for the calculation.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Spring & Spring & Spring & Autumn & Autumn & Autumn\tabularnewline
\midrule
\endhead
Altitude * latitude & Altitude * climate zone & Latitude * climate zone
& Altitude * latitude & Altitude * climate zone & Latitude * climate
zone\tabularnewline
-0.0295 & 0.530 & -0.362 & -0.0649 & 0.548 & -0.306\tabularnewline
\bottomrule
\end{longtable}

Suppl 7. Variance inflation factor (VIF) (3 SF) between elevation,
latitude and climate zone in spring and autumn. All years and camera
sites were used for the calculations. The datasets were subset so that
only one observation per site was used for the calculation.\selectlanguage{english}
\begin{longtable}[]{@{}llllll@{}}
\toprule
Spring & Spring & Spring & Autumn & Autumn & Autumn\tabularnewline
\midrule
\endhead
Altitude & Latitude & Climate zone & Altitude & Latitude & Climate
zone\tabularnewline
1.47 & 1.20 & 1.69 & 1.45 & 1.12 & 1.60\tabularnewline
\bottomrule
\end{longtable}

Suppl 8. The number of years that observations were obtained in for
every camera trap location that obtained a mountain hare observation.\selectlanguage{english}
\begin{longtable}[]{@{}ll@{}}
\toprule
Number of camera trap locations & Number of years samples were obtained
in\tabularnewline
\midrule
\endhead
304 & 1\tabularnewline
180 & 2\tabularnewline
116 & 3\tabularnewline
39 & 4\tabularnewline
30 & 5\tabularnewline
4 & 6\tabularnewline
1 & 7\tabularnewline
2 & 8\tabularnewline
2 & 9\tabularnewline
\bottomrule
\end{longtable}

Suppl 9. The ordinal days in which moult was initiated and finished in
spring and autumn. In spring, moult initiation was defined as the day on
which less than 90\% of the hare was white and moult end was defined as
the day on which less than 90\% of the hare was white. In autumn, moult
initiation was defined as the day on which less than 90\% of the hare
was brown and moult end was defined as the day on which less than 90\%
of the hare was brown.\selectlanguage{english}
\begin{longtable}[]{@{}lllllll@{}}
\toprule
Year & Ordinal day when spring moult started & Ordinal day when spring
moult finished & Days taken to moult in spring & Ordinal day when autumn
moult started & Ordinal day when autumn moult finished & Days taken to
moult in autumn\tabularnewline
\midrule
\endhead
2011 & 74 & 124 & 50 & 292 & 330 & 38\tabularnewline
2012 & 75 & 135 & 60 & 278 & 329 & 51\tabularnewline
2013 & 101 & 165 & 64 & 262 & 334 & 72\tabularnewline
2014 & 88 & 142 & 54 & 266 & 326 & 60\tabularnewline
2015 & 90 & 141 & 51 & 273 & 326 & 53\tabularnewline
2016 & 92 & 143 & 51 & 284 & 332 & 48\tabularnewline
2017 & 84 & 136 & 52 & 278 & 331 & 53\tabularnewline
2018 & 100 & 150 & 50 & 281 & 332 & 51\tabularnewline
2019 & 96 & 142 & 46 & NA & NA & NA\tabularnewline
\bottomrule
\end{longtable}

\textbf{Hosted file}

\verb`image16.gif`  available at \url{https://authorea.com/users/603769/articles/633959-altitude-latitude-and-climate-zone-as-determinants-of-mountain-hare-lepus-timidus-coat-colour-change}

Suppl 10. Prediction map animation with a resolution of 1 x 1
km\textsuperscript{2} giving the probability of mountain hares being
white across Norway on every day of the year. The probability of being
white was predicted using the environmental covariates contained in each
cell and the model output.

\selectlanguage{english}
\FloatBarrier
\end{document}

