Size-at-maturity shift in a male-only ﬁshery: factors affecting molt-type outcomes in Newfoundland and Labrador snow crab ( Chionoecetes opilio )

Size-at-maturity shift in a male-only ﬁshery: factors affecting molt-type outcomes in Newfoundland and Labrador snow crab ( Chionoecetes opilio A sex-asymmetric downward shift in size-at-terminal-molt has recently occurred in males in some portions of the Newfoundland and Labrador (NL) snow crab stock range, a ﬁrst known occurrence for such processes in snow crab ( Chionoecetes opilio ) stocks. This study exam-ines plausible factors promoting the shift in size-at-terminal-molt [synonymous with size-at-maturity (SaM)] including individual size, temperature, population density, and sex ratio. Analyses highlight expanse of cold water and large male density as being signiﬁcant predictors of molt-type outcomes. A conﬂuence of cold conditions and low density of large males promoted the SaM shift. In turn, the low male density was associated with recently elevated ﬁshery exploitation rates under quota-controlled management. It remains unknown the extent to which the reduction in terminal size reﬂects a phenotypic vs. genotypic process. Factors affecting skip-molting in male snow crab are investigated, and we ﬁnd that skip-molting occurs most frequently under extreme cold and high population density conditions. Potential complications arising from altered growth dynamics are discussed. Overall, the results advance knowledge on intraspeciﬁc competition processes within snow crab populations and inform ﬁsheries management systems that male-only harvest strategies do not provide full protection from biological harm to aquatic resources through ﬁshing.


Introduction
Context Snow crab (Chionoecetes opilio) has formed the basis of lucrative commercial fisheries in northern basins of the Atlantic and Pacific Oceans for decades and is currently undergoing range expansions into portions of the Arctic Ocean. In recent decades, the single largest fishery for snow crab has occurred in the Province of Newfoundland and Labrador (NL), Canada, where a fleet size numbering thousands of vessels has landed 30 000 to 70 000 tonnes of snow crab annually since 1994 (DFO, 2019a).
A broad-scale decline in size-at-terminal-molt, referring to the final molt and the point at which the animal stops growing, has been demonstrated in snow crab males during recent scientific assessments of the NL snow crab stock, with no such consistent or corresponding pattern in females (DFO, 2019b). From a biological perspective, this is a disconcerting observation as such sexasymmetric downward shifts in size-at-terminal-molt have never knowingly been documented in any other snow crab fishery resource.
Size-at-maturity (SaM) is a common metric of biological fitness used in the assessment of fish and shellfish resources. In snow crab, size-at-morphometric maturity is synonymous with size-at-terminal-molt (Conan and Comeau, 1986). Males, but not females, become sexually mature before the terminal-molt . Both sexes undergo morphometric changes during the terminal-molt, including enlarged abdomens in females and enlarged claws in males. The focal metric of SaM examined herein refers to "size-at-morphometric maturity" as V C International Council for the Exploration of the Sea 2020. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/ licenses/by/4.0/), which permits unrestricted reuse, distribution, and reproduction in any medium, provided the original work is properly cited.
opposed to "size-at-sexual maturity" more common in analyses of finfish.
The snow crab fishery targets only large males, so an obvious short-term outcome of the SaM shift is a decreased level of fishery yield per-recruit, but other biological and socioeconomic ramifications are possible should the SaM shift persist. Accordingly, it is important to investigate possible causes of the SaM shift in NL snow crab males.
In finfish, SaM is known to vary both within and across marine stocks in general. Populations in relatively stable or equilibrium situations normally show larger SaM than depleted populations (Trippel et al., 1997). In turn, population depletion is often associated with overfishing (Policansky, 1993;Trippel, 1995;Hutchings, 2004). Given that the fishery and SaM shift in NL snow crab are both restricted to males, the fishery is invoked as a likely factor contributing to the biological outcome of altered growth and maturity dynamics.

Growth and maturity dynamics
In species such as snow crab where competition among males for mating opportunities is high, evolution has promoted sexual size dimorphism, with males growing larger than females (Parker, 1992). This feature underpins management strategies for snow crab by enabling a male-only harvest. Other convenient biological management attributes include male polygyny (Sainte-Marie and Hazel, 1992), female sperm storage (Elner and Beninger, 1995;Sainte-Marie and Carriere, 1995), and sub-legal-sized sexually mature males in the population that help maintain stock reproductive capacity in the presence of size-selective harvesting (Ennis et al., 1990;Kruse, 1993;Orensanz et al., 2005;Sato, 2012). All factors considered, snow crab can theoretically be harvested aggressively relative to many aquatic resources, particularly for fisheries that more holistically target components of the spawning stock biomass including females.
SaM is ultimately an important regulator of reproductive output and population growth in brachyuran crab species (Hines, 1982(Hines, , 1989. Large adult males are generally superior to small adult males as breeders because they are most capable of clasping females (Claxton et al., 1994), carry larger sperm reserves, and afford better protection to females during the mating process (Rondeau and Sainte-Marie, 2001;Sainte-Marie et al., 2008;Sato, 2012).
Growth in crab species is a stepwise process associated with molting, with three possible molt-type outcomes in snow crab. Only adolescent (male) or immature (female) crab can molt. A "regular-molt" entails a crab remaining adolescent(m)/immature(f), a "skip-molt" entails a crab foregoing molting at an expected period of time and remaining adolescent(m)/immature(f), and a "terminal-molt" involves a crab undergoing a final molt into adulthood (morphometric maturity).
SaM has been shown to vary geographically in female snow crab in the eastern Bering Sea (Somerton, 1981). Moreover, in the closely related tanner crab (Chionoecetes bairdi) in the eastern Bering Sea, SaM has been shown to vary with geographical location (Somerton, 1981). Furthermore, there has been a previously described downward shift in SaM over time in both males and females of eastern Bering Sea tanner crab (Zheng, 2008). A spatial and bathymetric redistribution of crab was associated with the downward shift in SaM in the eastern Bering Sea tanner crab, suggesting that it was an environmentally driven phenotypic response, particularly as it occurred in both sexes.
Spatial variation in SaM in brachyuran crab may reflect several abiotic factors including latitude or temperature (Hines, 1989). Snow crab occupying cold areas often terminally molt at smaller sizes than those in warm areas. This phenomenon has been demonstrated in both sexes in West Greenland (Burmeister and Sainte-Marie, 2010) and in males in Atlantic Canada (Dawe et al., 2012). The smaller terminal size in cold areas more broadly reflects an overall inhibiting effect of temperature on molt frequency. For example, along with promoting small terminal size, cold conditions also appear to promote a higher incidence of skip-molting within snow crab populations (Taylor et al., 1994;Dawe et al., 2012;Murphy, 2019). Accordingly, the literature on molting dynamics suggests that thermal conditions are a plausible factor contributing to the recent SaM shift, particularly with a priori knowledge that recent years have been relatively cold (Mullowney et al., 2019).
Beyond temperature, molt-type probabilities in snow crab males are known to be affected by individual size, with large adolescents most prone to either skip-or terminal-molt in any given year (Dawe et al., 2012). Some literature also suggests that population density may affect molt-type probabilities (Comeau et al., 1998), with high densities associated with high skip-molt incidence (Hébert et al., 2017). Thus, by extension, population density was deemed a mechanism plausibly affecting SaM reaction norms, operating under a hypothesis that high population density would promote large SaM due to elevated intraspecific competition. This hypothesis is consistent with a suggestion by Elner and Beninger (1995) that lack of mating success as adolescents could promote further molts to large size. Furthermore, subsequent observations by Comeau et al. (1998) from an analysis of mating patterns in a western Newfoundland fjord from 1984 to 1993 showed large adult males to out-compete smaller adult and adolescent males for mating opportunities.
Population density of large males is directly negatively affected by fishing. Long-standing biological theory suggests top-down intraspecific competition release through size-selective harvesting can promote downward shifts in body size and maturation dynamics within excessively harvested animal populations in general (Miller, 1957;Policansky, 1993;Fenberg and Roy, 2008). Indeed, this phenomenon has occurred in NL fisheries resources before, specifically during the collapse of the iconic Northern Cod (Gadus morhua) stock in the early 1990s (Hutchings, 2004;Olsen et al., 2004). Such "unnatural selection" outcomes have never knowingly occurred but are theoretically possible in male-only crab fisheries (Carver et al., 2005;Sainte-Marie et al., 2008;Sato, 2012).

Objectives
The primary objective of this study is to investigate plausible factor(s) affecting male SaM in NL snow crab, with particular emphasis on events occurring during the recent SaM shift. The study focuses on demonstrated or theoretically supported determinants of molt-type outcomes. Given the focus on molting dynamics, a secondary objective is to investigate factors affecting skip-molting incidence in male snow crab. Skip-molting constitutes an intermediary interruption process between molts and was recently identified as an important confound to consider in population modelling of the Alaskan snow crab stock (Murphy, 2019). It is therefore important to better understand factors affecting skipmolting in NL snow crab. The outcomes of this work could have important implications for fisheries management with respect to how to best manage snow crab and similar male-only crustacean fisheries resources both in respect of climate regimes and harvest strategies.

Study area
The study was carried out on snow crab distributed in the easternmost Atlantic Canadian territorial waters along the Newfoundland and Labrador Shelf. The area ranges from southern Labrador in the north to the south coast of the island of Newfoundland in the south in spatial units known as the Northwest Atlantic Fisheries Organization (NAFO) Divisions 2H, 2J, 3K, 3L, 3N, and 3O ( Figure 1) and represents the northern range limit of Atlantic Canadian snow crab. The resource within this area encompassing 160 000 square nautical miles of continental shelf constitutes one genetic stock that more broadly spans throughout Atlantic Canada (Puebla et al., 2008). Although broadly synchronous, intrinsic spatiotemporal differences in population trajectories do occur across divisions within the study area (Mullowney et al., 2019).
Spatial units conforming to Assessment Divisions (ADs) used in the annual stock assessment were used in analyses. These ADs loosely reflect unique and contrasting bathymetric and thermal habitat conditions across ADs. The two northernmost ADs (2HJ and 3K) are characterized by nearshore shallow and cold plateaus bordering deeper, warmer offshore peripheries. In both ADs, large male snow crab are most typically found in channels between the nearshore plateaus and offshore banks. This distribution of large males reflects down-slope ontogenetic movements following settlement in shallower, colder areas (Mullowney et al., 2018b). AD 2HJ is overall colder than AD 3K, normally having a higher proportion of shallow grounds, particularly in the offshore. Throughout the region, shallow bottom areas are covered by the Cold Intermediate Layer (CIL), a body of <0 C water that sits intermediate in the water column. The southern AD (3LNO) comprised NAFO Divisions 3L, 3N, and 3O and is overall larger, shallower, and colder than the two northern ADs. This AD is known as the Grand Bank. The bottom in this AD is consistently and extensively covered by the CIL and this AD supports the highest biomasses of snow crab in NL (Mullowney et al., 2019). AD 3LNO has accounted for 75-80% of both snow crab biomass and landings in NL over the past decade (Mullowney et al., 2019;DFO, 2019a,b).
Along with habitat differences, fishery exploitation rates provide contrast across the three ADs, with annual exploitation rates typically highest in AD 2HJ and lowest in AD 3LNO (Mullowney et al., 2019;DFO, 2019b). Annual quotas and associated fishery exploitation rates are set by the Department of Fisheries and Oceans Canada (DFO) through a co-management consultation process (DFO, 2019a;Mullowney et al., 2020). The southern (NAFO Division 3P) and western (NAFO Division 4R) areas of NL marine shelves were not considered in analysis due to data deficiencies.

Surveys
Crab data were obtained from annual autumn offshore multispecies bottom trawl surveys in each AD. The surveys began in 1995. The surveys follow a depth-stratified random sampling design. A Campelen 1800 shrimp trawl is towed at a speed of 3 knots for 15 min wherever possible, with species-specific catches standardized to the swept area of each tow. Three Canadian Coast Guard (CCG) research trawlers (Teleost, Alfred Needler, Wilfred Templeman) have conducted the surveys since inception. The surveys occurred from September to December in most years. Two trawlers approach the survey from different ends (Divisions 2H and 3O) in September each year and converge towards 3K and 3L, typically finishing up in mid-December. In a routine year, about 600-800 trawl tows occur over the 2HJ3KLNO survey area.
Historically, inshore strata (within bays) were occupied, but that has seldom occurred in the past decade; thus, the data were restricted to offshore strata routinely occupied over the time series. Strata included in the survey footprint have a minimum depth of 55 m. Areas shallower than 55 m would not be expected to have high densities of snow crab (Mullowney et al., 2018b). The omission of inshore strata eliminated 0, 17, and 4% of the stratified (i.e. "trawlable") area within ADs 2HJ, 3K, and 3LNO, respectively. To assess potential habitat-related bias introduced by omission of inshore strata from the survey footprint, we calculated area-weighted mean depths for unoccupied vs. occupied areas. Further survey details can be found in Mullowney et al. (2018aMullowney et al. ( , 2019.

SaM estimation
Catches from each trawl set were either fully or sub-sampled atsea depending on their magnitude. A total of 152 000 males and 60 000 females measured during the 25-year period from 1995 to Size-at-maturity shift in a male-only fishery 2019 were included in the study. For all samples, carapace widths (CW) were measured (1 mm) and shell condition was subjectively assigned to a five-scale index that approximated time since molting: soft, new, intermediate, old, very-old (Mullowney et al., 2019). Female maturity status was determined based on visual examination of the abdomen, with enlarged abdominal flaps indicative of mature animals. To determine maturity status of males, the right claw of each crab was measured for height (CH, 0.1 mm) for subsequent discrimination of adult (terminally molted) vs. adolescent (not terminally molted) individuals based on a claw-carapace allometric function (1) used in the annual stock assessment. In the following equation, claw heights above d were classed as adult and those below d were classed as adolescent: d ¼ ð0:0806 Â ðCW 1:1999 Þ: (1) Following the first few years of life when multiple molts occur annually, crab are expected to be on a near-annual molting schedule . Towards focusing on the most recent molt in SaM estimation, an assumption that carapaces would not advance beyond a new-shell stage within a year was invoked and data were restricted to soft-and new-shell crab.
It is recognized that shell condition classification is a subjective process with inherent error (Fonseca et al., 2008;Murphy, 2019). It is possible that crab identified as new-shell molted more than a year ago and that crab identified as intermediate-shell molted less than a year ago. We deemed potential confounds of including intermediate-shell crab in the analysis to be greater than potential violations of our assumption that soft-and new-shell crab had most recently undertaken a molt. The omission of intermediateand old-shell crab also removed skip-molters (adolescent males or immature females with an intermediate or old shell) from SaM estimation.
A Generalized Additive Mixed Model (GAMM) was used to estimate size-specific proportions of crab that either terminally or regularly molted in any given year. Smallest (<40-mm CW males and <15-mm CW females) and largest (>150-mm CW males and >80-mm CW females) crab were excluded from the model due to a virtual absence of terminally molted crab at smallest sizes or adolescent (male) or immature (female) crab at largest sizes.
To control for small size-specific sample sizes, data were partitioned into 8-mm CW bins with a bin-centred size increment used to inform the model. The response variable was size-specific proportions of crab that were terminally molted and a binomial family distribution was assumed. Separate AD and sex-specific models were run using the mgcv package (Wood, 2017) in R version 3.6.0 (R Core Team, 2019). The model is expressed in the following equation: where the response variable M is the maturity (adult vs. adolescent/immature) of an individual of a given CW in a given year, b o is the intercept, s denotes a thin-plate smoothing spline, ti denotes an interaction (tensor) spline, f Year denotes a random factor effect of year, and e i is error. For female model runs featuring fewer size bins, the complexity of the smoothing matrix for the size term was set to a lower k ¼ 4 value than for males (k ¼ 6). Similarly, the complexity of the interaction term (CW Â Year) was set to a lower level (k ¼ 8) for females than for males (k ¼ 12). Model estimation was conducted via restricted maximum likelihood. Sufficiency of model fits was assessed based on distributions of size-specific residuals with years pooled as well as explained deviance statistics. For presentation, annual point estimates of the predicted sizes at which 50% of the crab terminally molted (mat50) were plotted and fit with loess regression curves. Similar estimates of predicted sizes at 33 and 67% maturity were calculated but not presented due to consistency in trend with the more conventional mat50 measure.

Factors affecting SaM
Annual mat50 size estimates for males from (2) were used in exploratory modelling to identify factors affecting SaM. Females were not included due to the observation of less synchrony in mat50 trends across ADs than in males and towards maintaining focus on the primary research objective of investigating reasons for the male SaM shift highlighted at recent stock assessments. Indices of abiotic and biotic variables were tested against the response variable of mat50 to try and explain processes affecting SaM. The first variable considered was a Habitat Index (HI) representing the areal extent of cold bottom water (<2 C) as per the stock assessment (Mullowney et al., 2019). For each survey set, bottom temperature data were collected by a conductivity, temperature, depth canister mounted on the headrope of the trawl. These temperature data were used along with all other bottom temperature data available from a suite of resource and oceanographic surveys conducted by DFO to linearly interpolate temperatures over space. Data capture in each AD reflected the summer to autumn period. The HIs are strongly negatively correlated with raw temperature means (not shown) but are preferred due to the spatial aspect of thermal conditions inherent in them.
Two measures of male population density were considered as possible explanatory variables affecting SaM; the densities of prerecruit and exploitable males. A pre-recruit is defined as an adolescent male of 60-94 mm CW and is expected to achieve legal size after another one or two molts. An exploitable male is an individual of !95-mm CW regardless of maturity status. We hypothesized that high densities of both groups would be positively associated with mat50 due to increased intraspecific competition that would be enhanced either in the present (i.e. by exploitable crab density) or in the near future (i.e. by pre-recruit crab density).
Annual indices of both exploitable and pre-recruit crab abundance were taken from the most recent stock assessment and divided by the area of bottom covered by the surveys to derive densities. Annual survey catches of each group were initially spatially expanded with Ogive Mapping (Evans et al., 2000) to derive indices of relative abundances. Subsequently, AD-specific catchability (q) adjustment factors were applied to increase knowingly under-estimated survey estimates into values closer to reality. These q factors (Mullowney et al., 2019) were developed based on time-series comparisons of survey exploitable biomass estimates to a complementary index developed from fishery catch rate depletion (Delury) modelling. A time-series median of annual q estimates in each AD was used to re-scale initial survey indices.
The q adjustment factors were calculated based solely on exploitable-sized crab (!95-mm CW). Accordingly, their application to pre-recruit males would likely result in an underestimation of true densities for pre-recruits because trawl catchability decreases as crab size decreases (Dawe et al., 2010). Nonetheless, despite neither the exploitable nor pre-recruit density indices being deemed absolute, and the pre-recruit index likely being under-estimated, application of the q adjustment factors brings indices for both groups into directly comparable scales across ADs (1000 crab per square nautical mile).
A fourth explanatory variable considered in modelling mat50 size was the ratio of exploitable males to mature females in the population. All exploitable male and mature female data were used regardless of shell condition. This sex ratio index was considered based on a supposition that if the sex ratio was skewed towards females, males would not have to grow as large to compete for mating opportunities as they would if the sex ratio was skewed towards males [i.e. in extension of Sainte-Marie et al. (1999)]. No survey q adjustment factors were applied to either group of crab prior to calculating the index due to concerns of using q factors from exploitable males to re-scale the much smaller females. Accordingly, the presented indices are directly comparable in scale within AD but not across ADs due to ADspecific differences in trawl catchability (Mullowney et al., 2019). The indices are also systematically skewed towards males due to the higher survey catchability for larger crab (Dawe et al., 2010).
All four explanatory indices were smoothed to 2-year moving averages ending in the terminal year. This was done for several reasons. First, there are concerns of "year effects" within trawlsurvey derived indices (Mullowney et al., 2018a(Mullowney et al., , 2019. Second, this smoothing was deemed to partially address subjectivity in shell aging concerns associated with calculating the mat50 response variable, specifically in the event that a new-shell crab molted more than 1 year ago (Fonseca et al., 2008). Finally, the 2year smoothing was deemed to consider that the "decision" to terminally molt vs. remain adolescent could be affected by ambient conditions not only in the most recent year but also by conditions during the preceding year.
Time series of the four explanatory indices were correlated with one another both within and across ADs using simple Pearson correlations to investigate redundancy prior to undertaking exploratory modelling on their relationships with mat50. An Exploitation Rate Index (ERI) of the fishery, defined as fishery landings divided by the most recent 2-year smoothed Exploitable Biomass Index, was also included in correlation analysis to investigate its relationship with density indices in particular.
A GAMM was used to investigate factors affecting mat50. This model, defined in (3), examined the response variable of mat50 size against continuous explanatory variables of crab density and thermal habitat indices. The input variables were lagged by 1 year to match the SaM outcome resulting from presumed molting in the spring or winter of the present year with conditions occurring in the preceding 1 and 2 years. Year and AD were included as random factor effects. The model including Exploitable Density (EXPDEN) as an input variable alongside HI was selected from three candidate models that individually used Pre-recruit Density (PREDEN) and Sex Ratio (SEXRAT) in lieu of EXPDEN towards formulating a parsimonious final explanatory model to investigate factors affecting molt-type outcomes. Model runs with PREDEN and SEXRAT showed no significance of either variable in affecting Mat50; thus, the EXPDEN model was chosen as the most appropriate model to use in attempting to explain molting and maturity dynamics.
where the response variable Mat50 is the CW at 50% maturity in a given AD and year for crab subjected to conditions of exploitable crab density (EXPDEN i ), HI (<2 C) coverage (HI i ). b o is the population intercept, s denotes a thin-plate smoothing spline, f Year denotes a random effect of year, f AD denotes a random effect of AD, and e i is error.

Factors affecting molt-type outcomes
A multinomial GAMM was used to regress and analyse molt-type outcomes in males (regular-, skip-, and terminal-molt) against plausible explanatory variables. As per all analyses, regular-(adolescent) and terminal-(adult) molters were restricted to soft-and new-shell males to identify individuals that most likely undertook a given type of molt in the most recent year. Skip-molters were defined as adolescent males in an intermediate-or older-shell condition. Variables selected for inclusion as explanatory inputs were those that emerged with most explanatory power from (3), specifically the HI and EXPDEN indices. A random spatiotemporal factor effect of a year and AD interaction was also included. Again, given the main effect indices were calculated as 2-year moving averages and subsequently lagged forward 1 year, it would be expected the analysis inherently captured signals from the conditions experienced by crab both during the preceding year and the year prior. The molt-type outcome model is defined in the following equation: where the response variable is the molt-type outcome (M) for a given crab of size CW subjected to preceding inter-molt conditions HI and EXPDEN. b o is the intercept, s denotes a thin-plate smoothing spline, ti denotes a tensor-interaction spline, f YearÂAD denotes an interactive random effect of year and AD, and e i is error. AD denotes AD. Model estimation was conducted via restricted maximum likelihood and the distribution of size-specific molt-type residuals at the AD level was used to determine the adequacy of fits.
To investigate directional impacts of the three main effects included as predictors in the molt-type outcome model (4), marginal effect plots were examined. The predicted values for each main effect were estimated by fitting the model with the other explanatory variables held constant. For the thermal habitat and exploitable male density indices, these constants were times-series means for each AD, while CW was held at 75-mm CW because this was a size near, which there was a high level of variability associated with proportions of crab undergoing a given type of molt in any given year.
Log-linear regression models incorporating the variable of male mat50 in relation to 1-year lagged HI and EXPDEN indices were used to qualitatively assess the extent to which observed outcomes within and across ADs have reflected the model-derived marginal effects of the explanatory variables. Similar models of EXPDEN in relation to fishery ERI were also constructed to Size-at-maturity shift in a male-only fishery explore the relationship between crab density and the controlling mechanism of fishing, with no lag applied in the comparison due to the lag inherent in the calculation of ERI (spring-summer fishery landings in year B/autumn survey biomass in year A).
Finally, the annual ratios of the total catch of sub-legal to legal adult males in the surveys within each AD were calculated and examined relative to time-series averages to qualitatively examine demographic structural changes that could potentially affect stock reproductive capacity moving forward. The intent of this investigation was to specifically examine if the ratio of adult males has skewed towards small adults in recent years, which are known to be inferior competitors for mates and likely of less benefit to breeding (e.g. Comeau et al., 1998;Sainte-Marie et al., 1999, 2008.

Surveys
Crab were captured in the trawl surveys throughout most of the offshore portions of the study region over the time series (Figure 1). The omitted inshore survey strata had area-weighted mean depths (61 standard deviation) of 236 6 65 and 201 6 101 m in ADs 3K and 3LNO, respectively. This compared to area-weighted mean depths of 492 6 310 and 279 6 336 m for surveyed strata in the two ADs, respectively. Accordingly, omitted areas (17% of area omitted in AD 3K and 4% of area omitted in AD 3LNO) systematically discriminated against shallow areas within each AD. Besides coastal areas and inshore areas omitted from analysis, only the shallow southeast portion of the Grand Bank in AD 3LNO was devoid of samples.
The bottom temperature profile for 1 May 2019 ( Figure 2) is representative of the general distribution in most years. Water on the shallow banks and nearshore plateaus is generally near or below 0 C while areas between banks and over the slope edges are warmer, typically 3-4 C. The shallow cold areas most closely reflect core habitat for small snow crab in particular, with warmer adjacent peripheries more commonly associated with large males (Mullowney et al., 2018b).

SaM estimation
A consistent pattern occurred in data clusters used to discriminate adolescent vs. adult males in each AD (Figure 3). In all ADs, smallest adult male crab were $40-mm CW. Adult males were routinely captured over the entire size spectrum above 40-mm CW, with 135-mm CW consistently representing a large male and 150-mm CW representing an exceptionally large animal. Adolescent males were infrequently found larger than about 110mm CW. As with all methods used to differentiate maturities from morphometric features in male crab, it is recognized that there is ambiguity in this assessment of maturities, particularly for observations occurring near the discrimination line.
The AD-and sex-specific GAMMs for determining SaM (2) fit the data well in all cases, with deviances explained ranging from 93.2% to 98.5% and a uniform distribution of size-specific residuals in both sexes in all ADs throughout the time series (Figure 4). The recent SaM shift in males is most apparent in AD 2HJ and least apparent in AD 3LNO, with AD 3K showing an intermediary pattern ( Figure 5). In AD 2HJ, mat50 estimates in males were at or above 80-mm CW from 1995 to 2014, with only one exception, but have since dropped substantially to range between 61-and 78-mm CW in the past 4 years. Similarly, following long-term consistency near the 95-mm CW legal-size mark, mat50 in AD 3K males has declined since 2014 and ranged from 77-to 90-mm CW over the past 5 years. AD 3LNO also underwent a reduction in mat50, with the six most recent years  consistently below the 95-mm CW legal size but is now on an upward trajectory. This recent pattern in AD 3LNO closely reflects that seen from 2005 to 2009, after which mat50 recovered and does not suggest that a major shift in mat50 in males has occurred in AD 3LNO. Converse to males, there has been no consistent broad-scale spatiotemporal declines in mat50 in females, with annual estimates in each AD fluctuating without pattern over the time series in each AD.

Factors affecting SaM
The magnitude of the HI was consistently higher in AD 3LNO than in the other two ADs throughout the time series, with >55% of AD 3LNO consistently covered by <2 C bottom water ( Figure 6). Trends in HI were synchronous in ADs 2HJ and 3K, with a correlation of r ¼ 0.70 (Table 1). However, the areal extent of cold water was typically higher in AD 2HJ in any given year, with a time-series mean of about 40% of the area covered in cold bottom water in AD 2HJ vs. a time-series mean of about 25% areal coverage in AD 3K. The trend in AD 3K was also significantly correlated with the trend in AD 3LNO (r ¼ 0.55, Table 1). In all ADs, record warm conditions occurred in 2011 but conditions cooled abruptly to the extent that the HIs were at or near timeseries highs during 2013-2017, with subsequent warming now once again occurring in ADs 2HJ and 3K.
Fishery ERIs have been highest and significantly correlated (r ¼ 0.64) in ADs 2HJ and 3K throughout the time series (Table 1 and Figure 6). However, in recent years, the fishery has been more aggressive in AD 2HJ, with the index (based on 2-year moving average biomass) averaging 51% since 2012 and featuring a minimum of 32% in 2015 and a high of 61% in 2018. Meanwhile, AD 3K ERIs have fluctuated from 36 to 48% since 2012. A progressively increasing trend in fishery ERI over the time series in AD 3LNO has become accelerated in recent years, with annual increases from 22% in 2012 to 51% in 2017.
Ratios of exploitable males to mature females are biased by higher trawl catchability on the larger males; thus, the ratios are relative and not absolute measures. Trends in the relative ratio of exploitable males to mature females in the population (SEXRAT) Figure 4. Residuals from SaM model (2). Predicted vs. observed (estimated) proportions of annual size-specific groups of crab (8-mm CW bins) that were identified as either being soft-shell or new-shell, by AD and sex.
Size-at-maturity shift in a male-only fishery were significantly correlated between ADs 2HJ and 3K (r ¼ 0.53) and between ADs 3K and 3LNO (r ¼ 0.76) ( Table 1). A dominant temporal pattern occurred in all ADs, featuring periods strongly skewed towards males from 1997 to 1999 and within the 2011-2015 time span depending on AD ( Figure 6). In recent years, all ADs have shown ratios less strongly skewed towards mature males.
Exploitable male density (EXPDEN) has oscillated in a declining trend since the late 1990s and fluctuated near time-series lows in all ADs for the past 5 years ( Figure 6). Exploitable crab have consistently been densest in AD 3LNO and least dense in AD 2HJ throughout the time series. The trends in exploitable crab density were significantly correlated (r ! 0.73) across all ADs.
Pre-recruit male crab density trends were significantly correlated across ADs 3K and 3LNO (r ¼ 0.79) but the trend in AD 2HJ was not correlated with either of the other ADs (Table 1). Pre-recruit crab densities have been consistently lowest in AD 2HJ and highest in AD 3LNO and remained below time series means in all ADs since 2011 ( Figure 6).
The exploitable and pre-recruit density indices were significantly correlated with one another in ADs 3K and 3LNO (r ! 0.74), but not in AD 2HJ (Table 1). Given their correlation in the largest ADs, these variables were approached separately in modelling. Expectedly, the fishery ERI was significantly negatively correlated with exploitable crab density in all ADs (r À0.43); thus, it was removed from further modelling in lieu of the biotic variable it directly affected, the exploitable male density.
The GAMM on variables affecting mat50 size of males (3) revealed that the habitat and exploitable male density indices were significant predictors of mat50 size (p < 0.05, Table 2). The model explained 54.4% of the deviance associated with the mat50 process.

Factors affecting molt-type outcomes
The GAMM on molt-type outcomes (4) showed that the interactions of CW with HI and EXPDEN were highly significant (p < 0.001) in affecting molt-type outcome, both in respect of the probabilities of either terminally-or skip-molting relative to regular-molting (Table 3). Together, the three variables in association with a significant year-AD interaction (p < 0.0001, Table 3) were able to reproduce the pattern of a shift in SaM found in the mat50 analysis, with a clear shift towards terminallymolting at smaller sizes in ADs 2HJ and 3K in recent years and a less pronounced shift in AD 3LNO (Figure 7). Note that CW range on the plot panels does not extend beyond 120 mm because few adolescent males were found above that size.
Notwithstanding occasionally sporadic high years, the model showed that there has been no clear systematic shift in proportions of skip-molting crab over the time series and revealed a spatial pattern of skip-molting incidence normally being lower in the northern ADs 2HJ and 3K (typically <15% of crab of any given size) and highest in southernmost AD 3LNO (often 10-30% of crab of any given size) (Figure 7). The model explained 24.8% of the deviance associated with the molt-type outcome process.
The direction of marginal effects in the GAMM molt-type outcome model showed that CW was consistently negatively associated with regular-molt and positively associated with terminalmolt in all ADs (Figure 8). The probability of regular-molt generally decreased as the areal expanse of cold water grew, although it became elevated at coverages beyond about 70% extent of <2 C water. The effect of thermal habitat on terminal vs. skip-molt was intriguing. Up to about a 55% areal extent of <2 C bottom water, crab progressively underwent terminal-molt in lieu of regular-molt. However, under most extreme conditions of >60% of the bottom being covered by <2 C water, skip-molting was equally or more common in lieu of terminally-molting, with >25% of crab predicted to skip-molt under conditions of 80% areal coverage of <2 C water in ADs 3K and 3LNO. Finally, in invoking outcomes other than regular-molt, the overall directional effects of high crab density was towards increased skip-molting and decreased terminal-molt as density increased. The marginal effect smooths were systematically linear (terminal-molt) or power shaped (skip-molt) as density changed for both molt-type outcomes.
Overall, relationships between mat50 and explanatory variables of either HI or EXPDEN were not fully consistent across ADs ( Figure 9). However, despite generally low coefficients of determination (R 2 0.28) for pairings at both the AD and overall levels, the systematic directional influence of the effects was evident, with a log-linear negative relationship of SaM with HI and loglinear positive relationship of SaM with EXPDEN. The unusually low mat50 estimates in AD 2HJ in recent years consistently fell outside the bounds of 95% confidence intervals in both analyses, suggesting that other factors could be affecting the recent major SaM shift in that AD, or potentially that the synergistic interaction of these two influential processes is substantially larger than either process can explain alone. Interestingly, in the consistently coldest AD 3LNO, the relationship of mat50 with EXPDEN was clearer than the relationship between mat50 and the HI.
In all ADs, the relationship between the ERI and EXPDEN was moderate to very strong (R 2 range 0.27-0.91), demonstrating the direct controlling influence of fishing on regulating the exploitable biomass (Figure 9).
Finally, since 2015, the ratio of sub-legal-sized to legal-sized adult males in the population has for the first time in the survey series become simultaneously skewed towards sub-legal-sized adult males (relative to time-series averages of the ratio) in all ADs (Figure 10).

Key findings
The results of this study are consistent with existing literature on factors affecting molting dynamics in snow crab, specifically that individual size is an important underlying determinant of molttype outcomes in males (Dawe et al., 2012) and that downward pressures on terminal size are promoted by cold conditions (Burmeister and Sainte-Marie, 2010). Support for existing theory Size-at-maturity shift in a male-only fishery of cold conditions promoting increased incidence of skip-molting in males (Taylor et al., 1994;Dawe et al., 2012;Murphy, 2019) is also present, with our results more specifically suggesting that it is in circumstances of extreme broad-based cold conditions (when the majority of an area is covered by <2 C bottom water) that skip-molting is most likely to occur.
Most importantly, this study advances new knowledge on the impacts of population density and by extension fisheries exploitation over molting dynamics in snow crab. Low densities of large males are shown to promote small terminal size, and it is via this route that fisheries can affect molting dynamics of male snow crab. This study constitutes the first known demonstration of evidence for the application of biological theory that "unnatural selection" associated with size-selective harvesting may affect demographic structure in a snow crab resource beyond immediate effects associated with removal of individuals.

Temperature
Overall, the study results support what is known in the literature about effects of temperature on molting dynamics in male snow crab. Our results are consistent in showing that cold conditions promote increased probabilities of both terminally-molting and skip-molting in males. Moreover, our finding that temperature was an important factor promoting the recent SaM shift in NL snow crab males supports the suppositions of Zheng (2008) in tanner crab that ambient environmental conditions affect growth dynamics in Chionoecetes crab species.
Our results on the effects of temperature on molting and growth dynamics add to the literature by demonstrating that temperature interacts with other factors to either buffer against or exacerbate changes in demographic structure within snow crab populations. For example, we showed a generally consistent broad-based phenomenon of cooling across a large spatial range in recent years, with the areal extent of cold bottom water near time-series highs in all examined ADs. Yet, despite general consistency in the cooling pattern, the extent of the downward shift in male SaM was spatially variable, arguably being far beyond (AD 2HJ), slightly beyond (AD 3K), or within the range of (AD 3LNO) what may have been reasonably expected from climateforcing alone based on historical relationships between the variables. The spatial differences in the magnitude of the SaM shift detailed herein were shown to be associated with differences in the population densities of exploitable males, consistently ranking from lowest to highest with the extent of SaM shift across the three ADs. The results show that the factor of population density consistently and systematically operated alongside temperature to regulate the extent of the SaM shift within each AD.
Our assertion that other factors interact with temperature to affect molt and growth dynamics is supported by the fact that no such SaM shifts occurred under similar or more extreme and prolonged cold periods historically, when population densities of   large males remained high within any given AD. The assertion is further supported by a lack of coupling of systematic trends in mat50 in the unfished females.
Ultimately, our results suggest that temperature is an important determinant affecting size and growth of snow crab but also that other factors interact with temperature to regulate molting dynamics. Size-at-maturity shift in a male-only fishery  Size-at-maturity shift in a male-only fishery

Population density
We suggest that the most important novel result of our study is a linkage between the density of large males and altered growth dynamics in snow crab. To date, relative to temperature, the regulatory influences of population density on molt and maturation dynamics have been less concretely advanced in the snow crab literature. Given the direction of the effect, with low densities of large male crab associated with small SaM, we interpret the recent shift in ADs 2HJ and 3K as reflecting intraspecific competition release that has lessened the need for any given individual to grow large to acquire sufficient food or spatial resources or possibly to successfully compete for mating opportunities. Large adult males are the most competitive crab in the population and are likely essential towards maintaining sufficient intraspecific competition to promote large SaM within populations.
Unlike temperature, population density of large males is at least a partially controllable variable in management of the resource. It is clear that fishing is a strong and direct controlling mechanism of exploitable male density within all examined ADs. Given that the mechanistic outcome of excessive fishing in some ADs in recent years, low population density, is invoked as a significant factor promoting a downward shift in SaM, this study constitutes the first known demonstration for the theory that unnatural selection can affect biological traits in male-only crab fisheries (e.g. Carver et al., 2005;Sainte-Marie et al., 2008;Sato, 2012). Application of unnatural selection theory extends beyond snow crab in more broadly challenging management theory of inherent precaution associated with intrinsic biological safeguards in male-only crab fisheries (Ennis et al., 1990;Kruse, 1993;Orensanz et al., 2005).
Beyond density of exploitable males, we also found evidence to support the idea that pre-recruit crab density could be an important variable affecting male SaM. Although only briefly addressed herein, this observation requires elaboration, as it is known that the impacts of fishing extend beyond exploitable males in this fishery. For example, in the ADs most consistently fished at relatively high ERIs, 2HJ and 3K, wastage in the form of capture and discards of pre-recruit soft-shell crab is a persistent concern (Mullowney et al., 2018a(Mullowney et al., , 2019. This is thought to reflect a lack of "buffering" competition by large adult males to prohibit softshell crab from successfully competing to access baited traps. ADs 2HJ and 3K are routinely associated with a depleted residual biomass (i.e. intermediate-or older-shell exploitable crab) and coincidental high levels of discarding in the fishery, primarily in the form of soft-shell pre-recruit males (Mullowney et al., 2018a(Mullowney et al., , 2019. Discard rates in these aggressively fished ADs routinely exceed 20% of the catch and can reach as high as >50% of the catch in extreme years. As many discarded soft-shell crab are believed to die after being released (Mullowney et al., 2019), pre-recruit density and subsequent recruitment into the exploitable biomass are impacted.
Finally, in the context of skip-molting, the evidence shown herein suggests that, when the density of large males remains high, there is an increased probability that adolescent crab foregoing a regular-molt will "choose" to skip-molt rather than terminal-molt. Notwithstanding a delay in progression of individuals through size within the population, this operates as another avenue via which maintaining a high population density benefits long-term recruitment potential for the fishery as well as the safeguarding of stock reproductive potential. As per general biological theory (i.e. Policansky, 1993;Trippel, 1995;Trippel et al., 1997), based on our results, we advocate that it is necessary to maintain sufficient intraspecific competition in the population at all times to maintain resiliency to exploitation of this important resource.

Fishing
Given the direct avenue by which fisheries exploitation affects both exploitable and pre-recruit male density, and the demonstration that offsetting impacts of low population density can affect biological functioning of the resource, it is reasonable to ask, particularly from a management perspective, what level of fisheries exploitation is too high? This is particularly relevant given that the other known factor affecting male SaM, temperature, cannot be directly controlled.
Overfishing is not an easy concept to define. Conventional iterations of the concept revolve around growth overfishing, whereby individuals are fished at smaller than optimal sizes, recruitment overfishing, whereby reproductive capacity and future recruitment to the stock are impacted (i.e. Ricker, 1975;Pauly, 1983), economic overfishing, whereby harvest levels are beyond those that yield optimal efficiency in long-term yields (i.e. Gordon, 1954), and ecosystem overfishing, whereby offsetting ecosystem impacts occur from overexploitation of a targeted species (i.e. Murawski, 2000). Based on the consistency of observations across our three major ADs, we propose 50% exploitation rate of legalsize males as a yardstick to define overfishing in this and potentially other snow crab stocks. AD 2HJ has frequently been fished above 50% of the estimated exploitable biomass during the past decade and AD 3K is consistently near this level. Where the phenomenon of mat50 shifts in males is most subtle, in AD 3LNO, a level of 50% ERI (based on 2-year moving average biomass) was only met in a single year.
Fisheries-induced shifts in SaM in male snow crab potentially meet aspects of all overfishing definitions. Beyond the recruitment overfishing associated with increased incidence of soft-shell pre-recruits in the catch, the phenomenon epitomizes growth overfishing, may have implications on long-term stock reproductive capacity through depletion of prime male breeders, and the rendering of a population of smaller individuals is likely to be of less economic interest and requires more individuals to be removed from the population to capture given quotas.
For context with our 50% exploitation rate definition of overfishing, in the precautionary approach management frameworks developed for other snow crab populations in Atlantic Canada, upper removal references are below 50% annual harvest rate (DFO, 2018a, b). In line with this, Siddeek et al. (2004) undertook exploratory modelling of snow crab in Alaska and the Gulf of St. Lawrence and found harvest rates in the range of about F ¼ 0.3-0.55 (annual exploitation rates of about 26-42%) to be best associated with Fmsy fishing strategies. Finally, Mullowney et al. (2018a) showed that 63% annual exploitation rate in NL snow crab populations was a level beyond, which fisheries performance was effectively guaranteed to fall below a level predicted by a simplistic climate model.

Potential outcomes
Potential problematic outcomes associated with reduced SaM may not be exclusive to the AD where fishing occurs, particularly with respect to upstream/downstream population connectivity. For example, Le Corre et al. (2019) recently demonstrated that a high level of upstream-sourced pelagic larvae of Northern Shrimp (Pandalus borealis) settle to bottom in southern portions of the NL shelf, in association with the dominant southerly flowing Labrador Current. Downstream reproductive concerns stemming from reduced SaM in males in northern portions of the NL shelf would be particularly concerning if reproductive potential of females has been compromised, such as could occur via sperm limitation. ADs 2HJ and 3K are the northernmost Canadian populations of snow crab within a stock range spanning not only NL waters, but all Atlantic Canadian waters. Potential downstream connectivity issues of degraded population health within the stock range may not even be exclusive to NL waters if they emerge.
Beyond potential immediate larval connectivity concerns, given that both sexual and harvest selection act on the same trait in this species (terminal size in males), it is possible that if excessive harvests and low densities of large males persist that more rapid genotypic selection could ensue than in species with less pronounced sexual size dimorphism, as truncation of males reduces variability in the sexual selection trait of size (Hutchings and Rowe, 2008;Sørdalen et al., 2018). This could compound long-term issues for the stock not only in areas such as ADs 2HJ and 3K, but throughout the stock range.
From a fisheries perspective, forthcoming yield potential can only be dampened with increased losses through males terminally molting at smaller sizes. This is not to suggest that increases to the exploitable biomass cannot occur moving forward. Indeed, there have been increases in pre-recruit abundance within localized areas of most ADs of the NL Shelf in recent years (Mullowney et al., 2019;DFO, 2019b). The SaM shift will impact biomass available to the fishery in the short-term, with potential to diminish long-term yield should it persist. Indeed, an example of foregone yields may have already occurred. In AD 2HJ, a large pulse of small crab <50-mm CW was present during 2013-2014 that have not subsequently progressed into either pre-recruit of exploitable crab despite sufficient time now elapsed for them to achieve these stages (Mullowney et al., 2019). The specific reasons for the loss of these crab are unknown and could include emigration, predation, high natural or fisheries mortality or the present demonstrated phenomenon of a high loss of crab to early terminal-molt.
From a reproductive capacity perspective, endogenous responses to small male SaM, such as sperm limitation (Jivoff, 2003;Carver et al., 2005;Sato, 2012), must now be considered. Mate selection is an important determinant for subsequent population-level expression of heritable traits (Sainte-Marie et al., 1999, 2008Hutchings and Rowe, 2008;Sørdalen et al., 2018), and the current trend towards reduced sexual size dimorphism in the northern ADs suggests the potential for reproductive interference. Smaller adult males are known to have inferior copulation and protection abilities for breeding females and carry smaller sperm reserves (Claxton et al., 1994;Sainte-Marie et al., 1999, 2008Rondeau and Sainte-Marie, 2001;Sato, 2012). A study on potential sperm limitation within the stock is presently on-going.
Finally, issues surrounding the increased dominance of small adult males as primary breeders in the northern portions of the stock range, where the major components of the northern cod (G. morhua) stock reside (Brattey et al., 2018), could become further exacerbated through increased losses to natural mortality, with small crab being most susceptible to predation (Chabot et al., 2008;Mullowney et al., 2014).

Biases and uncertainties
Despite linking reduced SaM to a synergistic combination of cold conditions and low large male density, the results do not suggest that these are the only factors that matter as important determinants of molt-type outcomes. The suite of predictive variables examined herein only explained a low amount of the underlying process of factors affecting molt-type outcomes. There are undoubtedly other unaccounted for factors such as food availability Size-at-maturity shift in a male-only fishery or predation impacts that need to be considered to more fully understand snow crab molting dynamics.
From a methodological and analytical perspective, the omission of some inshore areas by the surveys in recent years introduces uncertainty. These uncertainties are associated with direct crab measurements as well as calculation of habitat and density indices. Nonetheless, given these omitted areas only constituted a small overall proportion of the survey and stock ranges, and that both deep and shallow water areas were omitted from coverage, the effects are not deemed likely to affect overall results or interpretation. This is particularly relevant to AD 2HJ, where the most major SaM shift has occurred, in that it has remained fully covered by the surveys.
Finally, perhaps the most important uncertainty is whether recent SaM reductions in the northern ADs are phenotypic or genotypic. We cannot conclude that changes in male SaM are genotypic at this point in time, particularly as the present outcome in SaM remains consistent with a plastic response to recent cold conditions and there is evidence of recovery towards larger SaM having either already (AD 3LNO) or beginning to occur (ADs 2HJ, 3K).
Despite incomplete resolution regarding the extent of genetic variability in the SaM process, it is likely that continued low densities of large males can only increase the probability of genotypic changes occurring as early-maturing individuals contribute more genes to subsequent generations. Indeed, the extent of plasticity in the SaM process could become better revealed over the next few years, as fishery ERIs in the dominant AD 3LNO were reduced in 2019, bottom temperatures warmed in most areas of the NL shelf in 2019, and a moderately strong pulse of pre-recruit adolescent male crab is approaching legal size in ADs 3K and 3LNO. Accordingly, this analysis will be important to revisit within the next few years.

Conclusions
The study concludes that a combination of cold conditions and low population density of large males contributed to a recent SaM reduction in male snow crab in the northern ADs of the NL snow crab stock range. The low population density was directly affected by persistent high levels of fisheries exploitation, particularly in the northernmost AD (2HJ). It cannot be determined if SaM shifts are of short-or long-term duration, with the extent of plasticity in the process requiring further examination in forthcoming years. The study advances upon existing knowledge of climate effects on snow crab molt and growth dynamics by demonstrating interacting pathways of effects in conjunction with low population density. The identified male SaM shift will result in dampened fisheries yield potential in the short term, with offsetting long-term outcomes to stock growth or reproductive dynamics and associated fisheries yield unknown. The important observation of an SaM shift in a male-only fishery informs global fisheries management systems that male-only harvest strategies do not always provide full protection from biological harm to aquatic resources through fishing. The key piece of management advice for snow crab and similar male-only fisheries is that it is advisable to focus management strategies on maintaining a high density of large males in populations at all times and particularly so when additional stressors are likely to be exerting downward pressures on growth or maturation dynamics. Finally, our findings support existing knowledge on skip-molting in male snow crab through a conclusion that the behaviour is largely an outcome of extreme cold environmental conditions and also finds evidence that reduced density of large males may contribute to decreased levels of skip-molting in favour of early terminal-molt.