Phenology, population size, and factors influencing variation in density of an endangered butterfly, the mottled duskywing Erynnis martialis

: Understanding the natural history and ecology of endangered species is critical for developing effective, evidence-based conservation and management plans. The mottled dusky-wing Erynnis martialis is a skipper butterfly inhabiting oak savanna, oak woodlands, alvars, and tallgrass prairie habitats containing their host plants, Ceanothus spp. Listed as endangered in Canada, few populations persist in Ontario and Manitoba, but there are no formal estimates of population size or what influences spatial variation in adult density, and we have limited knowledge of phenology. To address these knowledge gaps, we conducted spatially explicit mark−re-sighting over multiple years at 2 of the largest known populations in Ontario (specific locations withheld). Population sizes at Site A were estimated to be 1159 individuals (confidence limit [CL] = 845−1598) in 2021 and at Site B, 626 individuals (CL = 466−851) in 2020 and 2227 individuals (CL = 1110−4463) in 2021. Most (81%) re-sightings occurred within 5 d of marking, suggesting that adults likely live less than 1 wk. High densities of adults were generally associated with high Ceanothus spp. abundance and open canopy, but relationships varied by site and year. Peak flight across sites and years ranged from 5−13 June. We recommend that existing and prospective habitat managers prevent succession with brush-cutting and/or prescribed burns, specifically in areas of abundant Cean-othus spp. This work provides valuable information for future management and recovery efforts, including planned translocations to both historically occupied and newly restored sites.


INTRODUCTION
For many threatened and endangered species, we lack detailed species accounts required for developing effective, evidence-based conservation and management plans. For example, without information on population size, phenology, and habitat preference, decisions related to the types of habitats to protect (Bergman 2001, Westin et al. 2018, how to optimally design reserve networks (Bedward et al. 1992, Zhang et al. 2012, what times of year habitats require protection (Cayrou & Céréghino 2005, Furey et al. 2011, and where to best conserve or restore corridors that link habitats (Cabeza & Moilanen 2001, Green et al. 2015 become extremely challenging. While some have argued that, under limited financial resources, funds directed towards conservation actions should take precedence over gathering basic ecological information (Whitten et al. 2001, Grantham et al. 2009), this view has been challenged (Lyles & May 1987, Sutherland et al. 2004, IUCN SSC 2013, Walsh et al. 2015, Toomey et al. 2017, Wood et al. 2018. Such challenges are justified, given that we often do not have enough information to determine which species are declining and at what rate or whether they should be listed as threatened or endangered (van Swaay et al. 2011, Swan et al. 2016, Taylor et al. 2017, IUCN 2019.
One strategy used to prevent species extinctions is conservation translocations: intentionally moving an animal from one area to another for the purposes of conservation (Seddon et al. 2007, Swan et al. 2016). However, successful translocations require highquality information about a target species' ecology (IUCN 1987, IUCN SSC 2013, Heikkinen et al. 2015, Taylor et al. 2017. For example, to select the appropriate donor population, it is important to have population size estimates to avoid removing individuals from an already small and vulnerable population. Population size estimates and demographic information from source populations can then be incorporated into population models that predict the outcome of reintroductions or other management practices, thereby informing the decision-making process for utilizing limited conservation funds (Sutherland 2006, Armstrong & Seddon 2008, Oberhauser et al. 2017, Wood et al. 2018). Monitoring of in situ outcomes can then be compared to model predictions and, where discrepancies arise, refinements can be made to support the development of more robust population models (Soutullo et al. 2008, Lewis et al. 2012. To improve conservation science, studies describing outcomes of translocations and other conservation management strategies, including population sizes and other basic ecological data for at-risk species, are critical. The goal of this study was to acquire vital ecological information on an endangered skipper butterfly, the mottled duskywing Erynnis martialis (COSEWIC 2012), slated for reintroductions in Ontario, Canada. Mottled duskywing inhabit oak savanna, openings in oak woodlands, alvars, and tallgrass prairie habitats that support their sole larval food source: Ceanothus spp. shrubs. Due to historical land clearing for European colonial settlement, agriculture, habitat degradation and fragmentation, natural succession of woody plants, competition with prolific exotic species, ex cessive deer browsing, pesticide use, and fire suppression (Linton 2015), these early successional habitat types have become extremely rare in southern Ontario (Bakowsky & Riley 1994). There are currently 2 geographic regions where the mottled dusky wing is present in Canada: the southeast portion of boreal Manitoba and the Great Lakes Plains in Ontario (COSEWIC 2012). The populations of mottled duskywing in Ontario's Great Lakes Plains have declined from 20 known locations extending across southern Ontario to 9 extant, geographically isolated populations in the last 30 yr (COSEWIC 2012, S. Kroeze unpubl. data). Mottled duskywing populations in boreal Manitoba also appear to have de clined over the last 20 yr, and this species has been extirpated from Quebec since 1958 (COSEWIC 2012). While site-based presence−absence information is useful and appears to be comprehensive for the mottled duskywing in Ontario, population sizes are unknown for this species in any part of its range (Linton 2015). Planning reintroductions of the mottled duskywing to formerly occupied areas will require data about their existing population sizes (to inform source population selection) and phenology (to understand basic biological interactions with their environment).
We sought to gather information about 2 existing populations of mottled duskywing in Ontario to (1) estimate population size and density, (2) de termine how environmental factors might influence spatial variation in adult density, and (3) de scribe sex-specific phenology over the flight period. Adult mottled duskywing enter their adult stage, mate, and lay eggs from mid-May to early July throughout their current range in the Great Lakes Plains. Historically, where populations existed in ex treme southern Ontario, a second adult brood occur red from mid-July to late August (Layberry et al. 1998). Mottled duskywing host plants include 2 species of Ceanothus shrubs: New Jersey tea C. ameri canus and Prairie redroot C. herbaceus (Lack 1954, Thomas 1991, Curtis et al. 2015, hereafter re ferred to collectively as Ceanothus. Females oviposit on Ceanothus leaves and larvae feed exclusively on these plants, forming leaf nests on them for shelter. Adults that emerge from pupation in early spring are thought to remain closely associated with the host plants for nectaring, patrolling, basking, and oviposition (Scudder 1869, Scott 1974, Layberry et al. 1998), but no empirical evidence exists about this association during the flight period. When mature, the larva either pupates and emerges in a few weeks (in locations where a second brood occurs) or overwinters as a diapausing larva until the following spring (in locations where there are single broods; Schweitzer et al. 2011). Overwintering larvae are usually 1−2 cm beneath the leaf litter, which exposes them to early spring warmth from the sun (Olson 2002).
To accomplish our 3 primary goals, we conducted a spatially explicit mark−re-sight (SECR) study at 2 sites located in the Great Lakes Plains in southeast-ern Ontario (the names and precise locations of sites are withheld because of the endangered status of this species). An SECR approach can provide information about sex-specific phenology and survival and is more precise than non-spatial mark-re-sighting/ recapture approaches because it includes information on the animal's location as part of the detection model (Efford & Fewster 2013). Using this SECR approach, we also examined the effect of Ceanothus density and canopy openness on spatial variation in density. In addition to the reliance of mottled duskywing on their host plant, canopy openness may also impact variation in adult density through thermoregulatory processes; particularly by mediating radiant heat, which acts as a trigger for adult butterflies to become active, for larvae to break diapause, and/or to trigger eclosion (Ohwaki et al. 2017, Grundel et al. 2020, Rather et al. 2021, Rija 2022.

Study sites
This study took place at 2 sites in the Great Lakes Plains of Ontario, Canada (COSEWIC 2012). Since the mottled duskywing Erynnis martialis is an endangered species targeted by collectors, the names and precise locations of 'Site A' and 'Site B' are not disclosed. The size of the suitable habitat monitored at Site A is 8.48 ha and characterized by oak savanna and tallgrass prairie, which has been intensively managed and restored. The size of the suitable habitat monitored at Site B is 8.03 ha and confined to the edges of a linear tract situated within a complex of alvar and mixed forest types with a high proportion of oaks. While closed-canopy forest borders most of the suitable mottled duskywing habitat, Ceanothus still grows along the tracts and into the edges of the forested areas a few metres beyond. Brush cutting is conducted annually, which provides the necessary disturbance regime that encourages Ceanothus growth by preventing succession. C. americanus is the only mottled duskywing host plant found at Site A, while Site B has both C. americanus and C. herbaceus.
At both sites, mottled duskywing individuals were monitored using an SECR approach (Efford 2004). To characterize variation in suitable habitat, we divided sites into 50 × 50 m grid squares using the GIS software ArcGIS (ESRI). Grid squares were only included if they had characteristics of suitable mottled duskywing habitat, such as (1) host plants, (2) other nectar plants, (3) hilltops, (4) puddling areas, and/or (5) appropriate ecological land class (mixed savanna, deciduous savanna, dry-fresh graminoid tallgrass prairie, open and shrub alvar, sugar maple−oak deciduous communities, or cedar−pine conifer communities). This left 44 marking grid squares for Site A and 29 grid squares at Site B. Unsuitable patches, such as dense forest with closed canopies, impenetrable shrubbery, or steep drop-offs, were excised from grid squares in an ArcMap shapefile that was imported into RStudio version: 2022.02.3+ 492.pro3 'Prairie Trillium' (RStudio Team 2022), using R v.4.2.0 (R Core Team 2022) to define the habitat mask applied in statistical modelling. In practice, some grid squares, particularly those at Site B where the habitat is primarily restricted to narrow areas beside the linear tracts, were much smaller than 50 × 50 m because of unsuitable habitat. A priori, little was known about the distribution of mottled duskywing at Site B, so marking grid squares were more spread out (averaging 130 m apart) than at Site A, where grid squares were immediately adjacent to one an other. The sampling protocol involved walking a sighting path in search of mottled duskywing (see Section 2.2 and Fig. S1 in the Supplement at www.int-res.com/ articles/suppl/n050p195_supp.pdf).

Mark−re-sight sampling and phenology
At both sites, mottled duskywing individuals were monitored using an SECR approach (Efford 2004). Mottled duskywing were individually marked and re-sighted to produce estimates of adult butterfly density and population size. Data were collected from Site B in 2020 and from both Site A and B in 2021. No data were collected in 2020 from Site A because of COVID-19-related access restrictions. Marking and re-sighting were conducted on separate and alternating days over the flight period from late May to early July. Either marking or re-sighting took place between 10:30 and 16:30 h, except for rainy and windy days when butterflies were inactive. When weather disrupted either a marking or resighting day, sampling was rescheduled for the following day and the alternating cycle continued.
Marking days consisted of visits to grid squares in a randomized order without replacement. To capture variation over the flight period, each grid square was visited multiple times. Each marking visit to a grid square involved a timed 2-person wandering search lasting 15 min, excluding butterfly processing time. The search consisted of 2 peo-ple meandering within the grid square at a slow and steady pace, keeping distance from one another to maximize potential butterfly encounters. Butterflies were netted by hand, and their location (UTM coordinates), sex, and be haviour were recorded. Upon capture, mottled dusky wing individuals were chilled on ice in a cooler to reduce activity, stress, and escape behaviour be fore being marked (Otis & Linton 2017). After cooling, the butterfly was given a uniquely identifiable set of 3 marks on the forewings (Fig. 1). Wing scales at the mark positions were gently rubbed off using a blunt-tipped probe to aid in mark adhesion. We used INSTA-DRI ® nail polish (Sally Hansen) for all markings.
Re-sighting occasions followed each marking day. A sighting path was established intersecting every marking grid square to gather daily transect counts (Pollard 1977, Pollard & Yates 1993. To account for time-of-day bias, we rotated where we started monitoring the sighting path between the 4 cardinal directions. During re-sighting, mottled duskywings were recorded if they flew within 5 m of the sighting path. Butterflies tend to fly constantly during peak temperatures; thus, it was sometimes necessary to capture and cool the butterflies to determine mark status. Otherwise, only binoculars were used to identify marking codes to minimize the stress from handling. When a marked or unmarked butterfly was encountered, their colour combination, sex, wing wear, behaviour, time, and nearest marking grid square were recorded. Re-sighted individuals with a missing mark had the code inferred using a process of elimination. On some occasions, the dorsal side of the wing that was marked was not visible before the individual flew off. In such cases, these sightings were recorded as 'unknown mark status', as they may have been marked or unmarked. Unknown mark status data contributed only to sighting day counts and was not incorporated into density models. Phenology was calculated from counts of marked, un marked, and individuals of unknown mark status that were encountered along the sighting path during re-sighting days. We report how the numbers of identified sightings of either sex changed throughout the flight period.

Detection covariates
Butterfly activity, and hence a butterfly's probability of being detected by an observer, can be influenced by daily temperature, cloud cover, wind, and moisture. Mottled duskywing are well camouflaged when close to the ground; thus, anything impeding flight could potentially reduce detection probability (Satter et al. 2019). Conversely, increased activity could lead to a higher probability of an individual being detected by observers. For example, when emergence peaks in the middle of the flight period, the high density of conspecifics could result in more conspicuous activities, such as 'spinning wheel' competition and mating flights. We tested whether detection probability changed based on weather conditions, over the flight period, or between sighting and marking days.

Habitat covariates
Within each marking grid square, we measured canopy openness and host plant density. For canopy openness, data were collected using a densiometer, with 96 representing maximum canopy openness and 0 being completely closed. Every marking grid square was split into 4 sub-quadrats, with canopy measurements taken from their centre. The 4 points were then averaged to derive one value for canopy per marking grid square. Canopy openness values were analyzed in 3 categories: 'dense canopy' (values from 0−32); 'mixed canopy' (33−64); and 'open canopy' (65−96). Host plant abundance was assessed in the field using a 3-person visual, subjective, categorical system: either low (0−50% coverage) or high (51−100% coverage). In cases where only 2 of 3 people agreed, we let majority rule; otherwise, there were no cases of all 3 people assessing differently. Confidence limits (CLs) for all estimates are reported at the 95% interval.

Statistical methods
SECR models use individual encounter history data to estimate population density while accounting for imperfect detection of marked animals during sampling (Royle & Young 2008). SECR models are hierarchical models, comprising an observation model that estimates the distance-dependent probability of detecting an animal nested within an ecological process model which estimates the number of activity or home range centres of individuals within the study area. In other words, the observation model identifies the environmental and human-observer factors that limit or promote the detectability of a butterfly given that it is present, while the process model estimates animal location based on the true ecological relationships between animals and their environment after accounting for variability in the detection process (Royle & Young 2008). A sum of the estimated number of activity centres in the study area provides an estimate of population density, and the population density multiplied by the study area provides an estimate of population size. Both the observation and the process model can be regressed against covariates to improve inference about the system.
The SECR model estimates population density (D) as a Poisson spatial point process of home range centres (Royle & Young 2008). Spatial variation in density can be related to underlying covariates as D(X; ) where is a vector of parameters for a model relating density to a location, X, specified by a vector of coordinates, x,y (UTM locations). The data for the model comprise detection histories (ω i ) for the n observed individuals for each occasion as a series of 0s (if the individual was not observed) or 1s (if the individual was observed). The entries of ω i are the individual de tections of individual i on S successive occasions at a set of known grid points (K ) -in our case, the centre of the sampled grid squares. The probability of observing a particular ω i depends on the detection parameters (θ) and the unknown home range centre of individual i (X i ) (Borchers & Efford 2008): (1) where L is likelihood, ω i > 0 indicates an animal that has been de tected at least once and f is the probability of density of home range centres, given that the animal was de tected. Integration occurs only over the habitat mask denoted above. Maximizing the likelihood provides estimates of θ (detection) and (den- sity) covariates (Borchers & Efford 2008). The generic form of the above equation provides the ability to calculate the probability (p sk ) of detecting a particular animal in grid square k on occasion s, which will depend on the location of the animals' range relative to the sampling square, the efficiency of detection while the animal is present in the grid square, and the shape of the home range (Efford 2011). We modeled the instantaneous location of an individual with a probability density function h(u), where u represents a vector of x,y locations. The h is commonly Gaussian or Laplacian bivariate kernels (e.g. half normal or negative exponential, respectively) with a single scale parameter, σ, which, as noted above, relates the distance between the location of the animals' range relative to the detection location. Efficiency of detection is controlled by the parameter g0, which is commonly interpreted as the probability of detection when the home range is centered on the grid square. Unlike territory centres, which would denote the centre of a defended space, we assumed that home range centres simply represented the centre of the total space used by an individual. The shape of home ranges was assumed to be circular by convention (Efford 2011). Further details about SECR models can be found in Borchers & Efford (2008) and Efford (2011).
Here, we used SECR to estimate mottled duskywing density, population size, and detection probability in relation to habitat covariates at 2 different sites. Three parameters are estimated directly using SECR: detection probability at the activity centre (i.e. g0, aka home range centre) and the spatial scale of movement around activity centres (σ) as well as the density of individuals per hectare (D). All models contained an individual random effect (ID). SECR analyses were conducted in RStudio version: 2022.02.3+492.pro3 'Prairie Trillium' (RStudio Team 2022), using R v.4.2.0 (R Core Team 2022) and the package 'secr' (Royle & Young 2008, Efford et al. 2009, Efford 2011, 2021, which estimates observation and point process model parameters using maximum likelihood. Assumptions of this method include that the population is demographically closed, those marked represent a demographically random sample, spatial sampling is random, markings are properly identified and not lost, and grid squares are subjected to a uniform search intensity.

Detection function and covariate modelling
The modelling process occurred in 2 steps to produce a final model that was used to infer population density. The first step was to parameterize a null (in-tercept-only) observation model to determine an effective sampling area (ESA), adjust the model for overdispersion, and identify the most appropriate shape function for detection probability (Table S1). Exponential, hazard rate, and half-normal detection models were compared using Akaike's information criterion (AIC). While the exponential function performed best, it could not produce parameter estimates; therefore, the half-normal detection model was selected and used in all subsequent models. The ESA informed the selection of a buffer width around the trapping array, which explicitly defined the area to which the density estimate applied. In instances where the suggested buffer width varied between years, the larger of the 2 was applied to both years for consistency. Next, we held density constant and compared several models using several different detection probability covariates to model the parameter g0, including occasion, time trend, and type of monitoring (either mark or re-sight). Occasion refers to betweenday variation in detection (g0 ~ occasion), time trend looks at linear variation in detection over time (g0 ~ time trend), and type of monitoring considers how the probability of detection may vary between active captures and passive observation (g0 ~ type of monitoring). All models were ranked using AIC, and the model with the lowest AIC was considered the most parsimonious model to explain detection probability. The variables to explain g0 in this best-fit detection model were then carried over to a similar model selection procedure to identify the best model to explain variation in density with the habitat covariates.

Density estimation, covariate modelling, and population size estimation
The second step in inferring population density of the mottled duskywing was to formulate the best explanatory model for density given the establishment of the detection function above. The model from the previous step assumed homogeneous density across the entire site and served as a null model, but we hypothesized that adult mottled duskywing density would vary according to the abundance of its host plant (Lack 1954, Thomas 1991, Brommer & Fred 1999, Holzschuh et al. 2013, Curtis et al. 2015 and canopy cover (Ohwaki et al. 2017, Grundel et al. 2020, Rather et al. 2021). Therefore, we built models including these features (D ~ Ceanothus, and D ~ canopy) and interactions between them (D ~ Ceano thus × canopy and D ~ Ceanothus + canopy). As be fore, all models were ranked using AIC, and the model with the lowest AICc was considered the most parsimonious model to explain spatial variation in density.
Using the top model, we present the results and plotted relationships among detection probability and density against the covariates. Density surface maps were produced describing how density varied across the site according to the most parsimonious model for each site and year. To estimate population size for each site−year combination, density estimates from each of the null models were multiplied by the total site area (Site A: 8.28 ha; Site B: 8.03 ha). CLs for all estimates are reported at the 95% interval. All analyses were conducted in RStudio version 2022.02.3+ 492.pro3 'Prairie Trillium' (RStudio Team 2022) using R v.4.2.0 (R Core Team 2022).

Mark−re-sighting and phenology
At Site A in 2021, a total of 232 individuals were marked over a period of 43 d (138 males, 89 females, 5 unknown sex). There were 35 re-sightings with complete identification and 1 re-sighting with partial identification. Four marked butterflies were resighted with one mark having fallen off the wing, but mark codes were able to be inferred. Thirty-one butterflies were re-sighted only once, while 2 were resighted on 2 occasions. A total of 656 unmarked individuals were sighted, while 14 were sighted with unknown mark status.
At Site B in 2020, of 141 mottled duskywing butterflies that were marked over 39 d (88 males, 52 females, 1 unknown sex), there were 43 re-sightings with complete identification and 3 re-sightings with partial identification. Five marked butterflies were re-sighted with marks having fallen off the wing (2 butterflies with 1 mark lost, 3 butterflies with 2 marks lost), but mark codes were able to be inferred for 4 of them. Thirty-three butterflies were re-sighted once, while 5 were re-sighted on 2 separate occasions. A total of 428 individuals were sighted without markings, and 36 were sighted with unknown mark status.
At Site B in 2021, of 96 mottled duskywing individuals that were marked over 38 d (59 males, 38 females), there were 6 re-sightings with complete identification and 1 re-sighting with partial identification. One marked butterfly was sighted with 1 mark having fallen off the wing, and the code was not able to be inferred. All 6 complete identification individuals were re-sighted only once. A total of 275 individuals were sighted without markings.
Across both sites and all years, 81% of re-sightings occurred within 5 d of marking, but the 3 oldest butter flies were re-sighted 11, 12, and 18 d after marking (Fig. S2). Peak emergence at Site A in 2021 oc curred on 8 June; at Site B, peak emergence occurred on 13 June in 2020 and 5 June in 2021, corresponding to 17, 21, and 15 d after the first sighting of the season, respectively. Across both sites and years, males emerged before females, resulting in the sex ratio being skewed towards males early in the flight period and switching to roughly even or skewed towards females later in the flight period (Fig. 2).

Detection, density estimates, population size estimates, and covariates
At both sites and all years, the top model describing variation in mottled duskywing detectability was 'occasion' (Table S2) (Fig. 3).
To examine the factors influencing mottled duskywing density, host plant density (i.e. Ceanothus) and canopy cover were incorporated into models for each site and year combination (Figs. 4,5 & 6) and ranked with AICc (Table S3). Site A in 2021, 'Ceanothus' was the top model, with greater Ceanothus abundance predicting higher mottled duskywing density ( Fig. 6; low Ceanothus = 2 ind. ha −1 , high Ceanothus = 254 ind. ha −1 ). At Site B, 'Ceanothus × canopy' was the best density model in 2020, with open canopy and high Ceanothus abundance predicting the highest mottled duskywing density (354 ind. ha −1 , CL = 159−789), and dense canopy with low Ceanothus abundance predicted to have lower mottled duskywing density (42 ind. ha −1 , CL = 18−98). At Site B in 2021, Ceanothus + canopy was the best model predicting mottled duskywing density, with the highest predicted densities also found in areas of high Ceanothus abundance and open canopy (1351 ind. ha −1 , CL = 438−4170) and lower predicted densities in areas of low Ceanothus abundance and closed canopy (17 ind. ha −1 , CL = 2−187).

DISCUSSION
We report the first effort to systematically monitor mottled duskywing Erynnis martialis populations over their entire flight period anywhere in their range. Our results provide density and population size estimates for what are thought to be 2 of the largest remaining populations in Canada. Our estimates suggest that the density at Site A (140 ind. ha −1 in 2021) and total population size (1159 in 2021) was comparable to Site B (density: 78 ind. ha −1 in 2020 and 269 ind. ha −1 in 2021; population sizes: 626 in 2020 and 2227 in 2021), supported by the overlapping confidence intervals for population size (Site A in 2021(Site A in : 845−1598Site B in 2020: 466−851;Site B in 2021Site B in : 1110. Population size estimates of declining species are useful for IUCN Red List ranking criteria (Cardoso et al. 2011, Fox et al. 2019). With consideration of interannual variation in insect population sizes, it will require a few additional years of monitoring to provide an actual ranking, to compare population sizes between sites, or to be able to report population trends describing exactly how susceptible mottled duskywing may be to stochastic or deterministic events that threaten their long-term persistence.
While both sites have small populations, their vulnerability to future loss or degradation of habitat differs substantially. Over the course of sampling at Site B, significant patches of host plant were altered or destroyed due to clearing for vehicle parking, untime ly brush-cutting of host plants in rights-of-way during the larval development stage, and herbicide sprays during the flight period. There has also been consideration of spraying Bacillus thuringiensis kurstaki (commonly known as 'Btk') insecticides for invasive Lymantria dispar (spongy moth) within mottled duskywing habitat (J. Linton pers. obs.). In contrast, the habitat at Site A is actively managed as a tallgrass system to prevent succession and degradation. Management practices include targeted re moval of inva- sive species and prescribed burns that occur annually or semi-annually in different areas of the site. Site A is also protected from encroachment of development.
All these factors suggest that the mottled duskywing has a higher probability of long-term persistence at Site A compared to Site B. Protection of the partially developed habitat at Site B will require cooperation from municipal authorities and residents, through education and enforcement of laws pertaining to species at risk. Moving beyond protection, land management practices utilized at Site A, namely prescribed burns, brush cutting, and targeted invasive species removal, should be considered as tools for maintaining, enhancing, or creating mottled duskywing habitat elsewhere in their range.
Our phenology data has important implications for future conservation translocations of this species. We have shown that the flight period lasts for approximately 6 wk, with peak emergence occurring 2−4 wk from the first sighting (Fig. 2). Being able to predict peak flight could be important for the augmentation of existing populations to ensure that translocated individuals have the highest probability of mating with individuals from the local population. We also 203 Fig. 4. Habitat covariates at Site A and the predicted mottled duskywing density surface. Habitat covariates (a) Ceanothus spp. abundance and (b) canopy openness were input in 'secr' software (Efford 2021) to derive the top-performing model for describing mottled duskywing density across the site. The top model in (c) 2021 'Ceanothus' was used to produce this predicted mottled duskywing density surface. Density shown as the predicted ind. ha -1 showed that female flight phenology tends to trail males by approximately 2 wk, which is essential information for the timing of efforts to collect mated females for captive rearing. Additional years of data on emergence and peak flight timing will help us to understand the influence that climate variables, like precipitation patterns and heat degree-day accumulation, have on emergence for this species. We found that the timing of peak flight varied by up to 4 wk over the 3 yr we monitored (at Site A: 30 June in 2019 and 8 June in 2021; at Site B: 13 June in 2020 and 5 June in 2021). Understanding the role of climate variables on emergence enables us to better predict timing each year for reintroduction planning and provide insight into how climate change could be affecting this species. Our work also demonstrates that both Ceanothus and canopy were important predictors of mottled duskywing density but that the strength of these predictors varied by site and year. Host plant abundance strongly predicted mottled duskywing density at both sites and in all years. The relationship between host plant abundance and adult density has been demonstrated in numerous butterfly species, such as the similarly monophagous and endangered frosted elfin Callophrys irus, among others (Swengel 1996, Kelly & Debinski 1998, Hanski 1999, Bierzychudek & Warner 2015, Curtis et al. 2015, Hill et al. 2018. However, by providing empirical evidence that adult mottled duskywing are primarily concentrated in areas where high densities of Ceanothus occur, our results point to the importance of having large Ceanothus patches if any habitat restoration projects targeted at this species are to be successful. Differences in the drivers of density between sites (just Ceanothus at Site A and both Ceanothus and canopy at Site B) can most likely be attributed to differences in their habitat composition and configuration. Site A contains abundant open canopy with substantial variation in Ceanothus abundance (Fig. 4).  (Efford 2021) to derive the top-performing model for describing mottled duskywing density across the site. The top models in 2020 (c) 'Ceanothus × canopy' and in 2021 (d) 'Ceanothus + canopy' were used to produce the predicted mottled duskywing density surfaces. Density shown as the predicted ind. ha -1 Thus, it appears that open canopy is not limiting at Site A, resulting in Ceanothus being the primary correlate with mottled duskywing density. In contrast, the primary limiting ecological factor at Site B is likely access to areas that provide sufficient sunlight for thermoregulation and mating. While there is also variation in Ceanothus abundance at Site B, its linear layout is such that all these strip-like Ceanothus patches are bordered by forest on at least one side, resulting in a substantial portion of the site being classified as mixed or dense canopy (Fig. 5). Management actions to open the canopy at Site B might have a positive effect on overall population size. More generally, for translocations and habitat restoration, our results point to the importance of not just having abundant Ceanothus, but also significant areas of open canopy. Such a recommendation is not necessarily trivial, as Ceanothus can also grow in shaded areas, as seen at Site B, as well as other areas of the province (A. Demarse, J. Linton & D. R. Norris pers. obs.). COSEWIC (2012) reports describe at least 5 other sites in Ontario where mottled duskywing have become extirpated in the last 30 yr. The likely cause of these disappearances relates to habitat declines through succession and fragmentation, possibly driven by canopy closure by mature trees and shrubs during later stages of succes-sion preventing butterfly movements between patches of suitable habitat. Active reintroductions are currently taking place at one of these sites, Pinery Provincial Park, and more re introductions are being considered at other previously occupied sites in Ontario. Based on the results of this study, we are better able to identify which potential release sites have the habitat characteristics that could be associated with higher densities of reintroduced mottled duskywing. Our density analysis, for example, suggests there will be very low mottled duskywing densities in areas of dense canopy cover, even where the host plant is abundant. Dense canopies could be inhibiting adult thermoregulatory processes associated with oviposition, larval development, and simply slowing movement. In addition to informing reintroductions, our results suggest it may be possible to expand and grow existing populations if connections are made to neighboring patches using habitat corridors, specifically by providing open and mixed canopy and abundant host plant. Furthermore, despite our findings that open canopy increases the density of adult mottled duskywing, variation in canopy cover might benefit immature life stages by offering particular microclimates (Lane & Andow 2003). Any significant habitat alterations, such as tree removal, should, therefore, be  Table S3 for estimates). At Site A in (a) 2021, the top model was Ceanothus; in (b) Site B 2020, the top model was an interaction between Ceanothus × canopy. At Site B in 2021 the top model was (c) an additive effect between Ceanothus and canopy. Estimates of low Ceanothus are sometimes too small to be seen on the graph. Error bars: confidence limits for all estimates, reported at the 95% interval managed adaptively by reviewing downstream effects on adult density. In all existing habitats, succession should be prevented using brush-cutting and/or prescribed burns. The karner blue Lycaeides melissa samuelis, eastern Persius duskywing E. persius, and frosted elfin are all butterfly species at risk, extirpated from Canada, that are also habitat specialists of tallgrass prairie and oak savanna. Although never common in southern Ontario, these habitats were historically more widespread, having been previously maintained by burns by Indigenous people (Bakowsky & Riley 1994), a practice that ought to be considered by habitat managers and studied further for its ability to prevent succession and canopy closure, promote growth of Ceanothus and other host plants, and for its impacts on overwintering larval survival. The results of our density analysis describe the habitat characteristics to which this habitat specialist butterfly species at risk responds: canopy cover and host plant abundance. It seems plausible that habitat management actions taken to improve tallgrass prairie and oak savanna habitats for mottled duskywing will also prove beneficial to other butterfly species at risk.
It is important to acknowledge our study may violate some of the assumptions of a spatial mark− recapture model. First, mark−re-sighting models rely on the assumption that marks are not lost. Unfortunately, 1 or 2 of 3 total marks put on each butterfly were seen having sloughed off the wings in about 10% of re-sightings at Site A, 18% of re-sightings at Site B in 2020, and 14% in 2021. When the assumption of no mark loss is violated, the data suggest that there are fewer marked individuals present than there really are, which can potentially lead to overestimates of population size (Arnason & Mills 1981, Malcolm-White et al. 2020). The application of 3 marks likely attenuated any real impacts of this issue on our study, since we were still able to identify marked individuals even if 2 of the 3 marks were lost. A second assumption that may have been violated is that the population is closed, which can also result in over-estimation of the actual population size. Any possible bias in over-estimation due to emigration is likely attenuated by the habitat matrix surrounding the study sites. The areas surrounding both study sites are unsuitable habitat, such as agriculture, development, and dense forests, which all lack the host plant, likely preventing high rates of immigration and emigration.Therefore, we suspect that any positive bias in population size estimates arising from the assumption of closed population is likely minimal.