Fishing and temperature effects on the size structure of exploited fish stocks

Size structure of fish stock plays an important role in maintaining sustainability of the population. Size distribution of an exploited stock is predicted to shift toward small individuals caused by size-selective fishing and/or warming; however, their relative contribution remains relatively unexplored. In addition, existing analyses on size structure have focused on univariate size-based indicators (SBIs), such as mean length, evenness of size classes, or the upper 95-percentile of the length frequency distribution; these approaches may not capture full information of size structure. To bridge the gap, we used the variation partitioning approach to examine how the size structure (composition of size classes) responded to fishing, warming and the interaction. We analyzed 28 exploited stocks in the West US, Alaska and North Sea. Our result shows fishing has the most prominent effect on the size structure of the exploited stocks. In addition, the fish stocks experienced higher variability in fishing is more responsive to the temperature effect in their size structure, suggesting that fishing may elevate the sensitivity of exploited stocks in responding to environmental effects. The variation partitioning approach provides complementary information to univariate SBIs in analyzing size structure.

frequency (L 95 ) to test fishing effect 35,36 , while the other used length class diversity to investigate the stability of population through time 37 . Recently, European Commission Marine Strategy Framework Directive required regional (or local) fishery reports to provide the information of basic SBIs (e.g. the mean length) in order to improve the management and maintain the sustainable development 38 . However, it remains unclear whether these univariate SBIs could represent the entire size structure and the status of a population. It has been suggested that no single SBI can represent an effective overall indicator for external forces 39 . Also, SBIs need to be selected carefully based on their implications. For example, L 95 can only reflect the variety of large fish in fish population. The analysis with the North Sea cod, herring and plaice found L 95 failed to reveal the effects of external forces on fish populations, as it was rather insensitive in responding to fishing mortality 36 .
To overcome the limitation of existing SBIs, we need an alternative approach to (1) analyze complete information of size structure and (2) examine how external forces affect the size structure. Here, we employed the variation partitioning approach to conduct a size structure-based analysis that examines the variation of size class composition in response to external forcings. Variation partitioning can be best understood as a method for extending multivariate regression. In multivariate regression (y ~ x), y represents a univariate response vector and x represents multiple predictors, x 1 , x 2 , etc. (each is a vector) and possibly their interactions; the contribution of each predictor variable (x i ) can be evaluated by partial R-square. Whereas in variation partitioning (Y ~ X), Y represents a multivariate matrix and X represents multiple predictors, X 1 , X 2 , etc. (each is a matrix); the contribution of each predictor matrix (X i ) and their interaction is also evaluated by partial R-square 40 . Variation partitioning is commonly used in community ecology to examine the relationship between species composition (Y matrix) and various sets of explaining variables (e.g. 2 or 3 predictor matrices) 40 . This method has also been extend to analyze temporal and space-time variation of community composition data 41 . Here, we borrow this concept to analyze temporal variation of size composition data in responding to fishing, temperature, and the interaction, with the simplification that fishing and temperature is just a vector. Specifically, for a given fish population, we apply the variation partitioning to quantify how the temporal variation of their size composition responds to fishing, temperature and their interaction (see Methods). The explained fraction of variation (partial R-square) by each factor then allows us comparing their relative contribution in affecting length composition through time.
Next, we perform a cross-stock meta-analysis linking the relative contribution of fishing or temperature (the output of variation partitioning as explained fraction of variation) to the life history traits of fishes, as well as long-term mean and variability of fishing or temperature (see Methods). This meta-analysis aims to examine which factor can explain the relative contribution of fishing, temperature and their interaction across stocks. This meta-analysis is motivated by the fact that life history traits are associated with the size structure of population 42 . We hypothesize that the large, slow growth, and late-matured species is more likely to be impacted by fishing in their size structure because size-selective removal (i.e. size truncation) may be more severe in these species 43 and their recovery will take longer time 42,44 . We also hypothesize that small species is more vulnerable to temperature effects, because smaller species are more sensitive to temperature changes due to the constrains from metabolic allometries 45 .
Furthermore, we expect that fishing and temperature might exhibit interactive effects via multiple ways 46,47 . For example, the long-term fishing effects, such as long-term mean and variability of mortality ratio (fishing mortality divided by natural mortality, F/M) may affect the relative contribution of temperature in explaining temporal variation of size structure. Here, we standardize fishing mortality by natural mortality in order to have a fair cross-stock comparison. Motivated by previous studies showing that fishing elevated sensitivity of exploited stocks to environmental changes 18,19,43 , we hypothesize that the fish stocks experienced higher fishing pressure is more responsive to temperature effect in their size structure. We also hypothesize that habitat conditions, including mean and variability of temperature, affect the relative contribution of fishing effect. For instance, Wang et al. 48 found that temperature affects the cod's life history trait, making the cod population more vulnerable to fishing.
Our objectives are, first, to apply variation partitioning to quantify how the variation of size structure responded to fishing, temperature and the interactive effects for 28 exploited stocks (Table 1) living in a wide range of habitats, including the west coast of US, Alaska, and North Sea ( Figure S1). Secondly, we linked the fraction of explained variation by fishing (or temperature) to life history traits (including von Bertalanffy growth rate (K), length infinity (L inf ), age at maturation (A 50 ), length at maturation (L 50 )), as well as long-term mean and variability of fishing and habitat temperature conditions. Finally, to demonstrate the efficacy of our approach, we compared the performance of variation partitioning approach with the traditional univariate SBIs analyses. This comparison is straightforward, as both the univariate SBIs analysis and variation partitioning are computed using similar linear modeling of variance/covariance, with the difference only in the response variable-the response variable is a vector (y) in the univariate SBIs analysis whereas the response variable is a matrix (Y) in the variation partitioning; two methods have the same explaining variables (i.e. fishing and temperature).

Results and Discussion
Our results of variation partitioning showed that the variance of size structure could be appreciably explained by fishing (on average of 10.9%), which is significantly higher than that of temperature and interaction (P = 0.019 and P = 0.00023 in ANOVA with paired t-test) (Fig. 1). Specifically, 12 out of 28 stocks were significantly affected by fishing while 7 stocks were significantly associated with temperature (Table 2); whereas, the interactive effect is small in most of stocks (Table 2).
We also found difference in fishing and temperature effect among regions (Fig. 2) and habitat types (Fig. 3). The fraction of variation explained by fishing is significantly different from temperature (P = 0.005) and interaction (P = 0.00046) in the west US. Both fishing and temperature effects are not significantly different from interaction in Alaska (P = 0.74 and P = 0.73) and North Sea (P = 0.27 and P = 1) (Fig. 2)  Results of ANOVA indicate that the fraction of variation explained by fishing is significantly higher than the temperature (P = 0.019) and interactive effect (P = 0.00023), but the fraction of variation explained by temperature is not significantly different from interaction (P = 0.344).  effect on size structure was significantly different from both temperature and interaction for demersal species (P = 0.0014 and P = 7.5 × 10 −5 , ANOVA with paired t-test) (Fig. 3). We employed the variation partition approach to analyze the time-series data of size structure of exploited stocks and found that fishing effect on size structure prevails ( Table 2). Such commensurate result goes with the study that analyzed the stock biomass of 28 stocks in the Northeast Atlantics 32 ; that is, in a large-scale, fishing is the most critical factor. Our analyses echo the increasing concern over effects of fishing on size structure of exploited stocks 13,26,47,49 .
We note however, many of these studies applied univariate size-based indicators (SBIs) as proxies for the change of fish size structure. Our approach incorporates the full size structure information in the analysis without assuming distribution of size data. Comparison between variation partitioning and univariate SBIs suggests that variation partitioning is more efficient in rejecting the null hypothesis (p-value is smaller) than the univariate SBIs (Fig. 4), particularly in detecting temperature effects (Probability of success = 0.63, P = 0.003) although only marginally in fishing (Probability of success = 0.56, P = 0.10). In other words, variation partitioning can be a useful complementary method to investigate the external forcings on the size structure of fish populations, in addition to univariate SBIs.
Our multi-stock meta-analytical framework also allows us to investigate how life history traits, and long-term mean and variability of fishing and temperature influence on the explained variation of size structure responding to fishing or temperature. The stocks experienced higher fishing variability (CV of mortality ratio) is more   (Table 3, Fig. 5d, P = 0.038). This supports our hypothesis that fishing elevates the sensitivity of exploited stocks in responding to environmental changes. Surprisingly, we found none of the life history traits is able to explain the fishing and temperature effect ( Figures S4, S6). Also, the total adjusted R 2 (Table 2) suggests that fishing plus temperature explained at most 44% of variation among 28 stocks. There may be other factors associated with body size, such as oxygen limitation of thermal tolerance 50 , affect the response to the fishing/temperature effect on the size structure.
Through the application of variation partitioning, we had expected the efficacy to identify interactive effect of fishing and temperature on the size structure. However unexpectedly, our results indicate that the interaction effect is very weak (Fig. 1). This finding may superficially be interpreted as evidence for lack-of interaction of fishing and temperature on size structure because fishing effects have dominated. However, we caution the interpretation of this finding, as variation partitioning is a linear variance decomposition method, which cannot account for nonlinear interactions.
While we demonstrate the efficacy of our size-structure based, meta-analytical framework to examine fishing and temperature effects on size composition of exploited stocks, we shall point out some caveats in our study. First, fishery-dependent data may lose some information due to the discard of small-size individuals. Second, there were fewer pelagic species than demersal and benthic species in our dataset, which might cause biased interpretation. Third, we assume an instantaneous response or a fixed lag in order to relate changes in size structure to changes in temperature, and we cannot provide detailed information concerning size class-specific responses.  Table 3. Results of univariate linear regression analysis on % variation explained by fishing or temperature versus each life history trait, mean mortality ratio (meanF_M), CV of mortality ratio (cvF_M)), mean temperature (meanTemp), and CV of temperature (cvTemp). This may be an important concern, because, for example, a warm year may lead to a strong year class and therefore affect the size structure of this year 25,51,52 , whereas the effect of temperature on the asymptotic body size may occur for several years.

Final Remark
We introduced the size-structure based approach relying on variation partitioning to quantify fishing and temperature effects on size composition of exploited fish stocks, instead of focusing on univariate size-based indicators. Through our multi-stock meta-analytical framework, we found that fishing explained most of the variation (Fig. 1), but difference existed between different regions (Fig. 2) and habitats (Fig. 3). We acknowledge that our analytical method still assumes linear responses of size structure to external forcings (as all univariate SBIs analyses do), because the nonlinear response forms remain unknown. Nevertheless, our analytical framework is a step toward better quantification of fishing and temperature effects on the size structure of exploited stocks in the context of life history theory. The information gained here may be useful for ecosystem-based approaches to fisheries.

Material and Methods
Size structure data of commercial stocks. We collected length frequency data from 28 exploited stocks, which contain temporal coverage over 20 years and annual fishing mortalities (or exploitation rates) estimated by stock assessment are available (Tables 1, S1). These stocks came from 3 regions in the northern hemisphere ( Figure S1): (1) the west coast of US (West US), which is part of the Northeast Pacific; (2) the Alaska region, which separates into 3 fishing areas-the Aleutian Islands (AI), Gulf of Alaska (GOA), and Bering Sea; and (3) the North Sea. Collectively, these 28 stocks in 3 regions were well studied, spanning distinct distribution of size structure, with a wide range of life-history traits and habitats, and therefore are representative of a compilation of global-scale fish stocks (Table S2).
We primarily use length frequency from the fisheries-independent surveys for each stock. These are the length frequency per size range of the given year as arranged in the stock assessment reports. For the Alaska region, the survey data include the bottom trawl surveys in the Aleutian Island (AI), East Bering Sea shelf (EBS) and Gulf of Alaska (GOA). For the North Sea, we compiled the length frequency data of the ICES International Bottom Trawl Survey in the North Sea (NS-IBTS). We used the 1 st quarter (Q1) in the NS-IBTS for consistency because there was only annual Q1 survey prior to 1991. For the west coast of US, we used fisheries-dependent length compositions instead of bottom trawl survey because the data from fishery-independent surveys during 1980-1990 were limited.
Fishing and natural mortality. To quantify the fishing effect, we used time series of annual fishing mortality from the stock assessment reports (Table S1). We focus on the single fishery whenever possible to minimize the uncertainty in fishing selectivity due to changing gears. We noted that analyses in the west coast of US use exploitation rate instead of fishing mortality in stock assessment. To make a fair comparison for meta-analysis, we transformed the exploitation rate into fishing mortality through the relationship between mortality and survival rate in fisheries 53 . Here, we first assume that these fisheries are type II fishery, in which the fishing and nature mortality operate concurrently. The exploitation rate (μ), fishing mortality (F), natural mortality (M), instantaneous total mortality rate (Z), and actual total mortality (A) have the following relationships: The equation (1) can also be written as: With the known exploitation rate and natural mortality, this equation can be solved numerically and yields the fishing mortality. The natural mortality here is mostly compiled from the value of the preferred model in the stock assessment reports (Tables S1, S2).
Temperature. In analysis, we primarily used empirical measurements of temperature along with the trawl surveys (see Figure S2). The North Sea (53-59°N, 3°W-10°E) near-bottom temperature is the station observations of hydrochemical measurements from the ICES Oceanographic database (http://ocean.ices.dk/HydChem/ HydChem.aspx) at Q1 (January-February), with coverage of almost the entire North Sea. For the Alaska region (Aleutian Islands, East Bering Sea and Gulf of Alaska) the bottom trawl surveys provide bottom temperature measurements (https://www.afsc.noaa.gov/RACE/groundfish/survey_data/data.htm). We took average for all the stations that the species occurred. Although many demersal species examined here have a pelagic phase at larval-juvenile stages and some species (e.g. Atlantic cod) forage on the whole water column, the preliminary analysis on surface (see Figure S3) and bottom temperature suggests that they are highly correlated (AI: 0.92, EBS: 0.72, GOA: 0.67, North Sea: 0.76). This suggests even the water is stratified, the interannual variability is very similar in the temperate-subarctic ocean. Therefore, we use only the bottom temperature collected in trawl surveys for these 4 areas in our analysis. For the west coast of US, we used the sea surface temperature from the NCEP/ NCAR Re-analysis monthly mean (http://www.esrl.noaa.gov/psd/data/timeseries/) and then took aerial average ( Variation partitioning to quantify fishing and temperature effect on size structure. To determine how much of the temporal variation in size structure is explained by the fishing and temperature effect, variation partitioning (aka. redundancy analysis) was used to decompose the total variation of length composition of each stock through years. The response variable here is the matrix of size composition data through time. The explanatory variables were time series of annual fishing mortality and temperature. The unbiased estimation of adjusted R 2 (accounting for sample size effect) provides a test with correct Type I error rate and good power for redundancy analysis 40 . We reported the adjusted R 2 for the pure effect of fishing, temperature and their interaction. Next, we perform permutation test (1000 times) to evaluate the significance of each fraction. Although the significance test cannot be done for the interactive component 40 , we still report the interactive component for the sake of comparison.
To incorporate the lagged effect of temperature, we additionally used 1 year-and 3 year-lag temperature as explanatory variables in variation partitioning. For the stocks in the Alaska region, sampling interval is 2-or 3-year, and thus we only considered 2 or 3-year lag. We did not investigate lagged effect of fishing; we assumed that fishing instantaneously affects the adult population while temperature may affect the size structure through influencing the future recruitment at early life stages (e.g. egg and early larvae). Because variation partitioning does not estimate log-likelihood in the procedure, the best model is selected according to the largest effective size (highest total adjusted R 2 (see Table 2)). The fraction of variation explained by fishing and temperature from the best model was used in further analyses.
Variables affecting the relative contribution of fishing and temperature effect. To investigate what determines the relative contribution of fishing and temperature across stocks, we considered the following variables (see Table S2): (1) life history traits: von Bertalanffy growth rate (K), length-at-infinity (L inf ), age at maturation (A 50 ), and length at maturation (L 50 ); (2) indices of fishing effect: long-term mean and variability of mortality ratio (i.e. fishing mortality/natural mortality, F/M); and (3) indices of temperature effect: long-term mean and variability of temperature. We used the mortality ratio (F/M) to reflect the fishing strength among different stocks. Because fisheries management usually sets optimal fishing mortality (F opt ) proportional to natural mortality (M) 54 , here we defined mortality ratio as fishing mortality normalized by natural mortality (F/M). When comparing multiple stocks with different natural mortalities, the mortality ratio can better reflect the impact of fishing pressure for cross-stock analyses.
To examine whether life history traits, long-term mean and variability of fishing and temperature have influenced the relative contribution of fishing and temperature effect, we first used simple univariate regression analysis to test each variable. We then build a linear mixed effect model (LMM) with the same variable as fixed effect and habitat as random effect to check if the observed pattern remains.
Comparing the performance of univariate SBIs to variation partitioning. To demonstrate the efficiency of variation partitioning, we compared the P value of explained fraction in the variation partitioning to the P value in regression with univariate SBIs. The univariate SBIs considered in the analysis were: 95% percentile of length class (L 95 ), mean length, Shannon diversity, Pielou's evenness index, and skewness. After calculating the univariate SBIs, each SBI was fitted to a regression model where SBI ~ temperature + fishing + temperature*fishing. If the P value of variation partitioning was lower than that of multilinear regression, it would suggest that variation partitioning is more efficient than univariate SBIs in terms of rejecting the null hypothesis.
Computation and Data availability. All computation was done with R version 3.3.3. Variation partitioning was carried out using R package vegan 55 . The linear mixed effect model was done with R package lme4 56 . Further model testing and model diagnostics were done with R package lmerTest 57 . The original data (except few stocks via personal communications, see Table S1 for details), the R scripts to carry out all analyses, and the outputs are available at: https://zenodo.org/record/1211120.