sm83s1-5015

INTRODUCTIONTop Operational ecological indicators are generally used for summarizing the status of marine communities and ecosystems in a comprehensive and accessible way (, ). Indicators of anthropogenic pressure are useful for monitoring the impact of fishing on target species, communities () and ecosystems (e.g. ), thus helping to assess the effectiveness of management measures (e.g. , ). Ecological and impact indicators are crucial in the design of multiannual management plans, which, according to the Common Fishery Policy (), should consider ecosystem descriptors in addition to fishing mortality targets at stock level. Indicators of fish community size structure can reveal the fishing effects caused by removing certain sizes and altering the abundance of different-sized species. Given their ability to summarize such complex dynamics, indicators of fish community size structure are increasingly used by national and transnational legislation, such as the European Marine Strategy Framework Directive (MSFD), which seeks to achieve Good Environmental Status (GES) by 2020 for all European seas (). Because no single indicator can capture the diversity of dynamics and processes within a system, suites of indicators are indeed required (), though the selection of indicators to retain should be carefully considered: indicators included should be sensitive to specific impacts () or be unambiguous to monitor single species, communities or ecosystems (, ). The suite of indicators should be based on criteria such as availability and accessibility of data, conceptual basis, ability to track different dynamics and processes, responsiveness, communicability and relevance to management (, ). The selected suite should also avoid redundancy (, ) and ensure applicability across different geographical areas (). Several studies have assessed existing indicators to identify suites that best capture the diversity of impacts and dynamics for fish communities or for whole ecosystems (). Most studies have focused on North Atlantic fisheries (, ) or compared multiple systems (e.g. ), while only a few have assessed the situation for the Mediterranean fish communities (, ). Within the MSFD, suites of indicators of GES are outlined for commercial fish (Descriptor 3) and food webs (Descriptor 4) (). The suite provided by the MSFD is considered a guideline that Member States must integrate with indicators relevant at regional and national level. Previous studies highlighted issues with one MSFD indicator, the large fish indicator (LFI; , ), in particular due to the required assumption about threshold level, which is area-specific. In the Mediterranean Sea, for example, found that the LFI is sensitive to the choice of species inclusion and to the threshold level. Recent research identified typical length (TyL), an indicator of the size structure of fish community, as a potential alternative which might provide similar results and avoid problematic assumptions (, , ). This paper aims to explore six indicators representing the size structure of the fish community: TyL, mean maximum length (MML), the LFI, the large species indicator (LSI), mean weight and evenness. This set includes the size structure indicators currently recommended by the MSFD (LFI and MML) and a selection of potential alternative or complementary length-based indicators for monitoring GES at national, regional or European level (, , ). In this paper we estimated these six indicators for 17 geographical sub-areas (GSAs) in the northern Mediterranean Sea, using data of length-frequency distributions obtained from the MEDITS experimental trawl survey. Given the emerging role of length-based indicators and their demand for complementing the MSFD at regional level, this study sets three goals. First, to evaluate which length-based indicators are essential and sufficient to explain the status and impacts of the Mediterranean Sea demersal fish community; second, to determine whether hidden common trends detected through dynamic factor analysis (DFA) across Mediterranean GSAs can be explained by basin-scale pressure indicators; and third, to determine whether TyL, proposed by ICES as an alternative to LFI, could be useful in the Mediterranean to complement or replace existing MSFD indicators. MATERIALS AND METHODSTop Data The MEDITS bottom trawl survey data used in this paper were collected in 17 GSAs (according to the GFCM classification, GSAs 1, 2, 5, 6, 7, 8, 9, 10, 11, 15, 16, 17, 18, 19, 20, 22-23 and 25) following a random depth-stratified sampling and using a gear with a common design (). Data relevant for this study were stratum, depth, latitude, longitude, number and weight of individuals caught by species, individual length, sex and maturity stage for a list of target species. The biological data were collected following a common protocol (MEDITS Handbook; ). We used the time series from 1999 to 2015, to include the largest number of GSAs and target species. The pressure indicators explored to explain common trends were environmental and exploitation proxies. In particular, we used the North Atlantic Oscillation (NAO) index, sea surface temperature (SST), the annual average SST anomaly of the Mediterranean Sea (Med anomaly), and as a proxy for fishing intensity, the fleet capacity. The NAO index was obtained through the Climatic Research Unit NAO website (https:/crudata.uea.ac.uk/~timo/datapages/naoi.htm). Med anomaly data were obtained from the European Environment Agency website (http://www.eea.europa.eu). The SST time series was obtained as average of the SST reprocessed in the Mediterranean Sea from the Copernicus website (http://www.copernicus.eu/). Data to derive fleet capacity (expressed as N of vessels×GT×Kw, where GT is the gross tonnage in tons and Kw is the engine power in kilowatt) were obtained from the fishing capacity time series of the Community Fishing Fleet Register (http://ec.europa.eu/fisheries/fleet/) after merging ports and landing sites at country and GSA level. The use of a combined metric for fleet capacity was applied to better differentiate pressure at sub-regional level and because it was more sensitive to changes over time. Fish community indicators based on size structure The analysis included data on bony and cartilaginous fish species, mostly but not limited to the commercially important ones, which have been MEDITS target species since 1999 (Supplementary material Table S1). The same species were used throughout the survey time series to avoid potential bias due to changes of target species in the MEDITS reference list over time. The indicators were calculated from the length-frequency distributions (N km^–2) according to . The size structure indicators considered in the analyses are summarized in Table 1. The resulting dataset include a time series for each of the six indicators in each of the 17 GSAs.

Table 1. – Indicators of size structure at community level. ¹, the threshold used for all the GSAs was 30 cm total length (TL) i.e. the value that returned the highest R² in the fitting of the 6th degree polynomial in 8 out of 17 GSAs (for 6 GSAs it was 40 cm TL and for 3 was 35 cm TL); ², in this study g was set = 45 cm TL, which is the 25th percentile of the L_max observed in the 17 GSAs under study for the species in Table S1 over the years.

Indicator	Semantic definition	Equation	Reference	Interpretation
Large fish indicator	Proportion of fish biomass larger than a set threshold.	$L F I (t) = \sum_{l > l b i g} y_{l} (t) / y (t)$ lbig threshold ¹; yl(t) catch per length class l; y(t) total catch (measured species).	ICES 2014ICES. 2014. Report of the Working Group on the Ecosystem Effects of Fishing Activties (WGECO), 8-15 April 2014, Copenhagen, Denmark. ICES CM 2014/ACOM:26. 174 pp. Shephard et al. 2011Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972. OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.	A decrease in the LFI could be due to increasing exploitation. Exploitation reduces the biomass contribution to the community of the larger individuals/species (Shephard et al. 2011Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972.).
Typical length	Geometric mean length of fish community, weighted by body mass.	$T y L = \exp [\frac{\sum_{i = 1}^{N} M_{i} \ln (L_{i})}{\sum_{i = 1}^{N} M_{i}}]$ N number of individuals in the distribution; M_i body mass of the i-th fish; L_i denotes the length of the i-th fish.	ICES 2014ICES. 2014. Report of the Working Group on the Ecosystem Effects of Fishing Activties (WGECO), 8-15 April 2014, Copenhagen, Denmark. ICES CM 2014/ACOM:26. 174 pp. OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.	TyL measures the size structure of fish and elasmobranch communities and decreases under high fishing pressure (OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019., Lynam and Rossberg 2017Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7.).
Evenness	Measure of the equitability in relative abundance among the length classes.	$E = \frac{H^{'}}{{H^{'}}_{\max}}$ ; $H^{'} = - \sum_{i = 1}^{S} p_{i} \ln p_{i}$ H′_max = ln(S) S is the number of lengths in the fish assemblage; p_i is the proportion of individuals in the i-th length class.	Pielou 1966Pielou E. 1966. The measurement of diversity in different types of biological collections. J. Theor. Biol. 13: 131-144.	Reduction in evenness should reflect increasing dominance of the community by small-bodied, fast-growing, highly productive species caused by increased fishing mortality (Greenstreet et al. 2012bGreenstreet S.P.R., Fraser H.M., Rogers S.I., et al. 2012b. Redundancy in metrics describing the composition, structure, and functioning of the North Sea demersal fish community. ICES J. Mar. Sci. 69: 8-22.).
Mean maximum length	Mean maximum length in the community.	$\bar{L_{\max}} = \sum_{j} (L_{\max j} N_{j}) / N$ L_{max j} is the maximum length obtained by species j; N_j is the number of individuals of species j; N is the total number of individuals.	ICES 2012ICES. 2012. Report of the Working Group on the Ecosystem Effects of Fishing Activities (WGECO), 11-18 April, Copenhagen, Denmark. ICES CM 2012/ACOM:26. 192 pp.	A decline in the MML indicates that the abundance of the most vulnerable fish and elasmobranch species is decreasing, leading to a change in the species composition (OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019., Lynam and Rossberg 2017Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7.) caused by increasing fishing pressure (Shin et al. 2005Shin Y-J., Rochet M.J., Jennings S., et al. 2005. Using size-based indicators to evaluate the ecosystem effects of fishing. ICES J. Mar. Sci. 62: 384-396.).
Large species indicator	Biomass proportion of large species in a community.	$L S I = \frac{\sum_{i} B_{i} (L_{\max} > g)}{\sum_{i} B_{i}}$ B_i is the biomass of species i having L_max > g ².	Shephard et al. 2011Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972.	A decline in LSI should indicate a decrease in the biomass of the predator species in the fish community (Shephard et al. 2011Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972.).
Mean weight	Mean weight in the community.	$M W = \frac{\sum_{i} B_{i}}{\sum_{i} N_{i}}$ B_i is the biomass index of species i; N_i is the abundance.	Rochet and Trenkel 2003Rochet M.J., Trenkel V.M. 2003. Which community indicators can measure the impact of fishing? A review and proposals. Can. J. Fish. Aquat. Sci. 60: 86-99.	A decreasing trend could be a signal of removal of larger individuals (Rochet and Trenkel 2003Rochet M.J., Trenkel V.M. 2003. Which community indicators can measure the impact of fishing? A review and proposals. Can. J. Fish. Aquat. Sci. 60: 86-99.).

Principal component analysis (PCA)

PCA was applied to select the minimum number of indicators to detect common trends among the GSAs and to reduce any redundancy in the six metrics considered. Redundancy may occur because several indicators are related to the relative abundance of large species. Data were first explored to check that the assumptions required to apply factor analysis were satisfied. The correlations among the six indicators for each GSA were tested using the Pearson coefficient. Next, data adequacy for the analysis was tested by the Kaiser-Meyer-Olkin (KMO) test (Kaiser 1974Kaiser H.F. 1974. An index of factor simplicity. Psychometrika 39: 31-36.). This test measures sampling adequacy for the model and for all variables. The statistic is a measure of the proportion of variance among variables that might be common variance. The higher the KMO, the more suited the data are considered (Table S2).

The assessment of redundancy was performed through the Bartlett test, which compares the observed correlation matrix with the identity matrix. If the variables are perfectly correlated, only one factor is sufficient. If all the variables are orthogonal, we need as many factors as variables. In order to select the number of factors for the PCA, the Kaiser criterion (Yong and Pearce 2013Yong A.G., Pearse S. 2013. A Beginner’s Guide to Factor Analysis: Focusing on Exploratory Factor Analysis. TQMP 9: 79-94.) was applied, retaining only the eigenvalues at least equal to one. Finally, the PCA was carried out including GSA and year as qualitative variables in order to detect separation among the GSAs according to the estimates of the six indicators and possibly a year effect.

Dynamic Factor Analysis (DFA)

DFA is a multivariate time-series technique used to detect M hidden common trends in a set of N time series (Zuur et al. 2003Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552.). This analysis was applied to the standardized (mean 0 and a standard deviation 1) indicators linked to the principal components of the PCA. For each indicator, a DFA was carried out among the estimates of the 17 GSAs. Each group of 17 times series was modelled according to a generalized formulation in terms of linear combinations of M common trends (from 1 to 3 trends), 4 explanatory variables (NAO index, Med anomaly, SST and fleet capacity), factor loadings (on trends and covariates), an offset and an observation error. Up to three explanatory variables were combined simultaneously. SST and Med anomaly were never combined due to their strong correlation. The mathematical formulation for this model is:

y_t = Zx_t + Dd_t + a_t + v_t

where yt is a Nx1 vector containing the values of the indicator for the 17 GSAs (N=17), x_t is an M×1 vector containing the values of the M common trends at time t, Z is an N×M matrix containing the factor loadings on the M common trends, d_t is an M×1 vector with the values of the four covariates at time t, D is a N×M matrix containing the factor loadings on covariates, at is the offset and v_t is the observation error distributed as a multivariate normal distribution with mean 0 and covariance matrix R.

Four different hypotheses on the covariance matrix of the observation error were tested:
– same variance and no covariance (diagonal-equal);
– different variances and no covariance (diagonal-unequal);
– same variance and covariance (equalvarcov); and
– different variances and covariance (unconstrained).

All combinations of covariates (included at their original scale) were considered for a total of 12 combinations (with 0, 1, 2 or 3 covariates). A total of 12×4 (covariance models) ×3 (number of trends) models was tested. The performance of the models was evaluated according to the corrected Aikake information criterion (AICc; Burnham and Anderson 2002Burnham K.P., Anderson D.R. 2002. Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. 2nd ed.. Springer, New York. ), following Zuur et al. (2003)Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552. and Holmes et al. (2018)Holmes E., Ward E., Wills K. 2018. MARSS: multivariate autoregressive state-space modeling. R package version 3.10.10..

The normality of residuals of the best model was tested through the Shapiro-Wilk test. In cases of deviation of residuals from the normal distribution, a scale finite mixture model (McLachlan and Peel 2000McLachlan G., Peel D. 2000. Finite Mixture Models. John Wiley & Sons. Hoboken, New Jersey.) was applied (Supplementary Material Section 5).

The investigation of the canonical correlations between detected common trends and explanatory variables followed Zuur et al. (2003)Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552., as another way of detecting their effects on dependent variables. For this reason, the relationship between the trends in the selected model (without covariates) and the covariates was tested with Pearson canonical cross correlations, also considering a time lag, set alternatively at 1, 5 or 10 years to account for lags between drivers and response due to biological processes. Ten years is the average lifespan (Table S3) of the demersal fish community weighed by their abundance (N km^–2). Since individuals subject to an impact are expected to be no longer in the population after their average lifespan, this metric can identify patterns attributable to a lasting change in the community, i.e. a trans-generational effect of the given impact. Biomass and abundance trends by species were compared with indicator trends to detect the species driving such trends. The mean temperature in the catch from Cheung et al. (2013)Cheung W.W.L., Watson R., Pauly D. 2013. Signature of ocean warming in global fisheries catch. Nature 497: 365-368. was used as a proxy for thermal preference of the species driving the trends to help the interpretation of the results (Table S4).

All the analyses were carried out in the R environment (R Core Team 2019R Core Team. 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.). The FactoMineR R library (Lê et al. 2008Lê S., Josse J., Husson F. 2008. FactoMineR: An R Package for Multivariate Analysis. J. Stat. Softw. 25: 1-18.) was used for PCA, the Hmisc R library (Harrell et al. 2006Harrell F.E. Jr. 2006. The Hmisc package.) for correlation analyses and the MARSS R library (Holmes et al. 2018Holmes E., Ward E., Wills K. 2018. MARSS: multivariate autoregressive state-space modeling. R package version 3.10.10.) for DFA.

RESULTSTop

Principal component analysis

The Pearson correlation coefficient showed several strong correlations among variables (Table S5). In almost all the GSAs the evenness was significantly and positively correlated with mean weight (p<0.05), while the LFI was significantly (p<0.05) and positively correlated with TyL, with the only exception of GSA 2 (Table S5).

For all GSAs the Kaiser criterion suggested the use of two factors, except for GSAs 1, 5, 15, 19 and 25, for which the use of only one factor was suggested. Thus, two factors were used for the PCA. The variance explained by the first component was 66.4%, while the variance explained by the second one was 20.9%, for a cumulative variance of 87.4% (Fig. 1 and Table 2).

Full size image

Fig. 1. – PCA results. A, distribution of the six indicators, TyL, LSI, MW, LFI, evenness and MML, on the principal component axes. The variance explained by each axis in % is reported along the two dimensions. B, scattered data points grouped by the 17 GSAs and 17 years used as qualitative covariates.

Table 2. – PCA results: correlation of indicators, GSA and years with two principal components. Explained deviance: first component 66.4%, second component 20.9%, cumulative 87.4%.

Variables	Correlation
Variables	Dim1	Dim2
TyL	0.94	–0.27
LSI	0.91	0.54
MW	0.84	–0.45
LFI	0.81	–0.25
Evenness	0.76	Not significant
MML	0.57	0.79
GSA	0.74	0.76
YEARS	Not significant	–0.17
% explained variance	66.4%	20.9%

Year did not seem to influence the grouping of the data (p<0.05), with all years clustered in the centre of the PCA space (Fig. 1B). Conversely, GSA was significantly correlated with both dimensions (Fig. 1B). The eastern and central Mediterranean (in particular GSAs 17, 18, 19, 20, 22+23 and 25) were characterized by the first coordinate negative and the second one between –2 and 1, while the western Mediterranean (specifically, GSAs 1, 2, 5, 7, 9 and 11) was characterized by positive first coordinate and second coordinate widely ranging (Fig. 1B). TyL had the highest correlation with the first principal component, while MML had the highest correlation with the second component.

The PCA results separated the six indicators into two groups: one represented by indicators sharing a definition explicitly based on the L_max parameter that focuses on the fraction of the community with large species (MML and LSI). The other group included indicators not based on L_max. The subsequent analyses were carried out considering TyL and MML, the indicators most correlated with the principal components and the least correlated with each other. Each indicator represented one of the groups detected in the PCA.

Dynamic factor analysis

The DFA results showed that the TyL trend decreased between 1999 and 2002 and increased steadily afterward, while the MML trend increased until 2004 before declining sharply between 2008 and 2012 (Fig. 2).

Full size image

Fig. 2. – Underlying common trends over the time period 1999-2015 detected by DFA for the indicators TyL and MML.

The comparison among the models tested for TyL and MML showed that for both indicators the covariates did not improve the goodness of fit, giving an AICc higher than the baseline model. In Table S8, the best models are reported for each type (defined by the included covariates). The best model is characterized by one underlying common trend (M=1) for each indicator with the covariance matrix “diagonal and equal”.

The common trends detected can be decomposed to effects at each GSA level through the factor loadings (Fig. 3), representing the contribution of each GSA to the common trend. In GSAs 2, 17, 19 and 22+23, TyL was poorly influenced by the common trend (with factor loading below 0.2 in absolute value), while GSAs 1, 5, 10, 16, 18 and 25 had a factor loading higher than 0.2 (Fig. 3). The factor loadings estimated for TyL on the common trend showed that in the western Mediterranean this indicator is associated with the trend (with the exception of the Balearic Islands, GSA 5), while in the central Mediterranean TyL was inversely associated with the common trend in GSAs 18 and 20, and positively in GSAs 16 and 25.

Full size image

Fig. 3. – Factor loadings on the hidden common trend for TyL and MML. Top panel: values of factor loadings for the two indicators per GSA. Bottom panel: barplots with individual bars representing GSAs with values above ±0.2 (marked by dotted horizontal lines) are considered indicative of a positive/negative effect (reported in bold in the table).

The factor loadings of the MML indicated that most GSAs were associated with the hidden common trend, except GSA 1, 15, 16 and 25 (Fig. 3). The MML in GSAs 2, 8, 11 and 22+23 was poorly correlated with the common trend, while GSAs 1, 5, 6, 16, 17, 18, 19, 20 and 25 had a factor loading higher than 0.2 in absolute value. The factor loadings for the two indicators showed a significant negative relationship (r=–0.68; p<0.01; Fig. 4).

Full size image

Fig. 4. – Factor loadings for the two indicators by GSA. Linear smoothing overlapped to show the linear relationship. Values above ±0.2, marked with black dashed lines, are considered indicative of a positive/negative effect.

Correlations between trends and covariates

The covariates SST and Med anomaly were highly correlated with each other, and both showed a strong negative correlation with the fleet capacity (Fig. S8 and Table S9), with opposite patterns in time (Fig. S9). NAO showed no correlation with other covariates and no clear temporal patterns. The common trends for TyL and MML detected by the DFA were compared with the covariates under consideration (SST, NAO, Med anomaly and fleet capacity) through canonical correlations between the detected trends and the covariates. A significantly strong positive correlation was found between the trend of TyL and the covariates Med anomaly and SST, while a negative correlation was found between TyL and fleet capacity (Table 3). The hidden trend detected for MML was significantly correlated with the same covariates, but inversely.

Table 3. – Canonical Pearson correlation between the common trend detected by DFA for the two indicators, TyL and MML, with no time lag and a 1-year, 5-year and 10-year lag. The significant (p<0.05) correlations are indicated in bold.

	Covariates	no lag	1-year lag	5-year lag	10-year lag
TyL	NAO	0.12	–0.11	–0.33	–0.37
	MedAnomaly	0.72	0.7	0.59	0.76
	Fleet_cap	–0.94	–0.9	–0.6	–0.72
	SST	0.62	0.58	0.6	0.75
MML	NAO	–0.23	0.02	0.4	0.2
	MedAnomaly	–0.69	–0.65	–0.32	–0.61
	Fleet_cap	0.92	0.89	0.53	0.45
	SST	–0.5	–0.44	–0.36	–0.59

Correlations between changes of species abundance and selected indicators

In GSAs 1, 7, 8, 9, 10, 11, 15, 16 and 25, where TyL increased consistently with the common trend, TyL was strongly influenced by three species characterized by higher biomass: Galeus melastomus, Pagellus erythrinus and Mullus barbatus (Figs S2 and S4).

In GSA 5 the increase in S. smaris (Tables S6 and S7) produced a decrease in MML (Fig. S3), this species being very important in terms of abundance (more than in biomass) in this GSA, while in GSA 1 the increase in G. melastomus and the decrease in P. acarne led an increase in MML. The abundance trend of M. merluccius in GSA 16 and 17 (respectively increasing and decreasing, Tables S6 and S7) influenced the trend of MML, while in GSA 18 the increasing trend of density of M. barbatus and S. flexuosa returned a decreasing trend of MML.

In GSA 19 the decreasing trend in MML was influenced by the increase in S. smaris, S. flexuosa, M. barbatus and P. acarne.

Contrasting patterns between areas were observed in the comparison between TyL and MML. In some GSAs, increasing trends of species showing both high abundance and biomass led to an increase in both indicators. This was the case for G. melastomus (in GSAs 1 and 10) and M. merluccius (in GSAs 10 and 16), both species with high L_max. Conversely, in GSA 18 and 19 the increase in small-sized and thermophile species such as M. barbatus and P. acarne led to a decrease in both MML and TyL, these species being the most abundant in both abundance and biomass and characterized by low L_max (Table S5 and S6). A more complex situation was found in GSA 6, where TyL increased but MML decreased. This could be attributed to the fact that both abundance and biomass of M. merluccius declined with a sharper decline in abundance (Figs S2 and S3), entailing an increase in mean weight, and hence length.

DISCUSSIONTop

We estimated a set of six size structure indicators of the fish community which are considered useful to detect changes at population and community level and are among those assessed by OSPAR (2017)OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019. and by the ICES Working Group on Ecosystem Effects of Fishing Activities (WGECO) (ICES 2014ICES. 2014. Report of the Working Group on the Ecosystem Effects of Fishing Activties (WGECO), 8-15 April 2014, Copenhagen, Denmark. ICES CM 2014/ACOM:26. 174 pp. , 2018ICES. 2018. Report of the Working Group on Ecosystem Effects of Fishing Activities (WGECO), 12-19 April 2018, San Pedro del Pinatar, Spain. ICES CM 2018/ACOM:27. 69 pp. ).

The first goal of this study was to address the redundancy within this set of proposed indicators in order to select those that effectively summarize the fish community dynamics over time. Previous studies used factor (Greenstreet et al. 2012bGreenstreet S.P.R., Fraser H.M., Rogers S.I., et al. 2012b. Redundancy in metrics describing the composition, structure, and functioning of the North Sea demersal fish community. ICES J. Mar. Sci. 69: 8-22.) and ordination techniques such as non-metric multidimensional scaling (Blanchard et al. 2010Blanchard J.L., Coll M., Trenkel V.M., et al. 2010. Trend analysis of indicators: a comparison of recent changes in the status of marine ecosystems around the world. ICES JMS. 67: 732-744.) to address indicator redundancy and compare the community size structure of different ecosystems. We, instead, combined different statistical techniques. Firstly, we used a factor technique (PCA) to select a non-redundant subset of two size-structure indicators of fish community. TyL and MML were found to be the indicators most linked to the principal components. Secondly, we explored the occurrence of common trends among the GSAs through DFA on the two selected indicators, following Zuur et al. (2003)Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552.. The selected models were characterized by only one common trend, while Zuur et al. (2003)Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552. found three common trends in their analysis. The capability of DFA to identify more common trends in our case might be limited by the spatial and temporal scale of the analysis: the fine-scale geographical patterns shown by PCA are not captured by DFA. Moreover, the relatively short time series lacking contrasted patterns and the fact that the indicators smooth over the trends of individuals species might limit the capability of capturing potential additional trends that may be masked by the main common trend. While the two selected indicators are expected to show community-wide improvement upon reduction of fishing pressure, here we found asynchronous and contrasting patterns. This result could be explained by the second goal of our study: to explore through DFA the relationships between the common trends detected among the GSAs for the two selected indicators and the times series of environmental and anthropogenic drivers at basin scale. While DFA showed no significant link between the indicators and the covariates, correlation analysis revealed significant relationships. As observed by Zuur et al. (2003)Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552., the fact that the tested models incorporating covariates were not selected as the best models does not necessarily mean that these covariates had no influence on the time series. The higher AICc for the models including covariates could indicate that the explanatory variables did not produce an overall improvement in the model, the covariates being related to only some of the time series. Indeed, for both selected indicators, the inclusion of each covariate implies the estimation of 17 extra parameters (one for each GSA), which could make the model worse for covariates not linked to all the time series. This implies that the relationship between trends and covariates is area-specific, as also confirmed by the fact that the factor loadings seem to show geographical patterns: TyL is related to the common trend positively in most western Mediterranean GSAs and negatively in the eastern Mediterranean ones (Fig. 3). Different patterns between the eastern-central Mediterranean and the western Mediterranean also emerged from the PCA (Fig. 1B). Similarly, sub-regional patterns were observed in the North Sea, where OSPAR (2017)OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019. highlighted that the increasing TyL trend observed in northern areas was not observed in southern areas.

The fleet capacity showed a strong negative correlation with the trend detected for TyL. This result indicates that the increase in TyL (Fig. S2) observed in 10 out of 17 GSAs (mostly in the western Mediterranean) could be associated with a decrease in fleet capacity. This effect could be confounded with that determined by the increasing SST: TyL is indeed linked to the biomass of the fish community and thus strongly influenced by the species with higher biomass, which in the western Mediterranean include thermophile species such as P. erythrinus and M. barbatus (Fig. S4; Cheung et al. 2013Cheung W.W.L., Watson R., Pauly D. 2013. Signature of ocean warming in global fisheries catch. Nature 497: 365-368.). Notably, in the last few decades, landings of thermophile species in the Mediterranean have increased, contrasting with a general decrease in landings and a widely reported impact of increasing temperature on species and ecosystems (e.g. Vasilakopoulos et al. 2017Vasilakopoulos P., Raitsos D.E., Tzanatos E. 2017. Resilience and regime shifts in a marine biodiversity hotspot. Sci. Rep. 7: 13647.). However, the increasing trend of biomass in our case is observed also for non-thermophile species such as G. melastomus, S. smaris and, to a smaller extent, L. piscatorius and L. budegassa. Therefore, the increasing common trend detected in TyL is most likely due to a decrease in fishing pressure, although a synergic effect with SST cannot be excluded for some species.

Conversely, MML showed a positive correlation with fleet capacity. Based on its definition, MML is expected to decline with an increase in fishing pressure, so the reduction in fishing effort observed in all GSAs is expected to result in an increase in MML. Instead, a decrease was observed for the general trend (Fig. 3) and for the majority of sub-areas (GSA 5, 6, 7, 9, 10, 17, 18, 19, 20; Fig. S7), indicating an asymmetric response of MML to increasing and decreasing pressures. MML is related to the abundance of species characterized by large maximum size (L_max), which generally also show slow growth and late maturation (Winemiller and Rose 1992Winemiller K.O., Rose K.A. 1992. Patterns of life-history diversification in North American fishes: implications for population regulation. Can. J. Fish. Aquat. Sci. 49: 2196-2218.). We can therefore expect a lag between a decrease in fishing pressure and an increase in MML. The decrease in the observed correlation with a time lag of 5 and 10 years (Table 3) is consistent with our hypothesis. The observed pattern could also indicate an increase in the abundance of species characterized by small L_max. These generally show a faster positive response to reduction in fishing pressure than species with large L_max, and thus an increase in abundance of the former leads in the short term to a decrease in MML. This pattern, known for other areas (OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.), emerges in individual examination of GSAs (Table S5 and S6), where the dynamics of the most abundant species, such as M. barbatus and S. smaris, explain indicator patterns. The trend of MML therefore seems to be driven by these two contrasting factors, which act synergistically to produce the observed asymmetric response to decreasing pressure. Contrary to the expectation, a reduction in fishing pressure leads to an initial decline of MML driven by an increase in small, short-living species caused by release of fishing pressure. Not before several years, an increase in MML can be expected as the large fish species recover. In our case, with a lag of 10 years, the trend shows signs of inversion, with loss of the positive correlation between decreased fishing pressure and decreased MML. While this pattern might be driven by an increase in small species with preference for high SST, in our case the small species that showed an increase in abundance include both thermophile species (e.g. M. barbatus) and non-thermophile species (e.g. S. smaris; Table S3), indicating that SST does not explain the increasing patterns for all the small species.

The declining pattern in MML observed in GSA 6 and attributed to the faster decline in abundance than in biomass of M. merluccius could be indicative of fishing pressure targeting smaller individuals (for example, if larger individuals are distributed in areas not reached by trawlers), but could also be driven by a succession of poor recruitments (Lynam and Rossberg 2017Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7.).

The two indicators selected by our analysis correspond to those currently proposed by OSPAR and ICES WGECO. MML, in particular, is recommended by the MSFD to assess objective 3.3, “population age and size distribution”, to account for the relative abundance of old, large fish in the community. Maximum length is considered a proxy for species vulnerability, as species with larger maximum length given their life history traits tend to be more vulnerable to fishing. These species are also expected to be the first to decline under high fishing pressure, with a decline in MML implying a decrease in abundance of the most vulnerable fish species (OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.).

A complementary indicator that can account for changes in size structure has traditionally been identified in the LFI (European Commission 2010European Commission. 2010. EU Commission Decision of 1st September 2010 on criteria and methodological standards on good environmental status of marine waters (2010/477/EU). Off. J. Eur. Union, L232: 12-24. , OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.), currently included for MSFD Descriptor 4 (namely, 4.2.1). Recent research has however identified some difficulties related to the implementation of the LFI for the MSFD. Critically, the determination of the threshold for what constitutes a “large fish” is prone to arbitrary choices. Furthermore, methods developed for estimation of this parameter have revealed sensitivity to the fish communities studied, to the species included and to the length of the time series (Shephard et al. 2011Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972., Spedicato et al. 2014Spedicato M.T., Bitetto I., Carbonara P., et al. 2014. Methodological support for exploring the LFI: “proportion of large fish - indicator” properties for interpreting the changes in the fish community of the South Adriatic area (GSA18). In: First MedSuit Regional Workshop on indicators and targets to ensure GES of commercially exploited marine populations in the GFCM area. FAO HQs, Rome, Italy, 6-7 November 2014. ), rendering a standardized definition between areas problematic. Despite the fact that we calculated LFI by applying a 6^th degree polynomial for optimal threshold determination, as suggested by Shephard et al. (2011)Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972., the LFI was not selected in our analysis. TyL has been suggested to replace or complement the LFI. While possessing similar statistical properties to the LFI, TyL is considered more robust to the LFI’s limitations (ICES 2014ICES. 2014. Report of the Working Group on the Ecosystem Effects of Fishing Activties (WGECO), 8-15 April 2014, Copenhagen, Denmark. ICES CM 2014/ACOM:26. 174 pp., Lynam and Rossberg 2017Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7., ICES 2018ICES. 2018. Report of the Working Group on Ecosystem Effects of Fishing Activities (WGECO), 12-19 April 2018, San Pedro del Pinatar, Spain. ICES CM 2018/ACOM:27. 69 pp. ).

Our analysis improves present comprehension about the performance of TyL and LFI for Mediterranean demersal fish community, showing that the TyL has lower overlap with other indicators, emerging as a potential alternative. Additionally, our study confirms the complementarity of TyL to MML, already highlighted by Lynam and Rossberg (2017)Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7., because the joint application of both indicators can help disentangle changes in community composition from changes in size structure.

Finally, we highlight three general aspects pertaining to the indicators’ capability to track the responses of the system to different exogenous and endogenous drivers. First, indicators are generally developed and defined to track community response to increasing pressure, whereas the effect of pressure release is rarely assessed and generally assumed to be inverse to pressure increase. We propose that indicators can behave asymmetrically to increasing and decreasing pressures: while most indicators can generally track a decline in community status due to increasing pressure, they might differ in their capability to capture community response to pressure release. In our study, for example, MML patterns might be indicative of a time-lagged response to decreasing pressure. There is limited understanding of indicators’ responsiveness to decreasing pressures (Lynam and Rossberg 2017Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7., OSPAR 2017OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.), and the choice of indicators rarely takes this characteristic explicitly into account, as shown by their definitions (Table 1). Our study suggests that MML might perform worse (i.e. slower) than TyL in capturing a community’s recovery after fishing pressure reduction. Nominal effort reduction is assumed to correspond to a reduction in fishing pressure in our study. Effort and fishing mortality, however, cannot be linearly related and the observed delay might be due to reduction in fishing mortality declining less rapidly than decline in effort. Nevertheless, our metric of fleet capacity emphasizes the capacity of larger, more powerful boats, which exert the bulk of the pressure on the stocks. Combined metrics of fleet capacity is thus a better proxy of fishing impact than other metrics of fishing effort.

Second, the complementarity between indicators in responsiveness should be accounted for. We show that the indicators’ complementarity should extend to the capability to account for impacts across different time lags: MML shows a likely time-lagged response to pressure reduction, while the trend of TyL is associated with fast response of large-sized species to changes in fishing pressure and/or environmental conditions. This aspect is critical, as the responsiveness of indicators is therefore vulnerable to the length of time series, masking relationships or leading to spurious relationships and potentially confounding the interpretation of indicators (see for example the positive correlation between MML and fleet capacity).

Third, we highlight how definitions do not consider multiple stressors, so combinations of SST and effort can have antagonistic/synergistic effects, leading to an unclear interpretation of the indicators’ trends. OSPAR (2017)OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019. noted that the influence of SST and fishing on TyL can be hard to disentangle. In our study, we were able to tell the two drivers apart by taking into account the thermal preference of increasing species. We noted that the increase in TyL could mainly be attributed to decline in fishing pressure rather than increase in SST or, at best, to synergistic effects for a few species. Similarly, also for MML we determined that the observed changes could not be univocally driven by increase in SST: in fact, the increase in abundance potentially driving the MML pattern was observed in small-sized thermophile and non-thermophile species alike. This would reinforce the hypothesis of the increase in small species following fishing pressure release. Only through knowledge of the underlying biology of the species can we disentangle drivers’ effects on species patterns.

Overall, these considerations reinforce the importance of assessing the indicators’ capability to capture trends across diverse ranges of time lags and to interpret their patterns according to multiple drivers, especially in light of changing communities driven by climatic change.

Our results strongly support the inclusion of TyL as one of the regional MSFD indicators for detecting changes in the demersal fish community throughout the Mediterranean Sea. The results of this study will therefore be of support to regional application of MSFD targets by reinforcing scientific knowledge of length-based indicators in Mediterranean Sea. Additionally, these results can be applied in the context of multiannual management plans, allowing community and ecosystem dynamics to be considered, as required by the Common Fishery Policy (EU Reg. 1380/2014). The identification of reference levels for these indicators remains a critical point to be addressed. For North Atlantic areas, reference levels for TyL, MML and the LFI have been identified (ICES 2018ICES. 2018. Report of the Working Group on Ecosystem Effects of Fishing Activities (WGECO), 12-19 April 2018, San Pedro del Pinatar, Spain. ICES CM 2018/ACOM:27. 69 pp. ). In the Mediterranean, these values do not yet seem to have been determined (Dupont et al. 2018Dupont C., Barsoumian S., Moreira G., et al. 2018. Article 16 Technical Assessment of the MSFD 2015 reporting on Programme of Measures. Italy Report.), and future studies should therefore focus on this issue at different spatial scales.

ACKNOWLEDGEMENTSTop

This study was conducted in the framework of the MEDITS Project. We are grateful to Cristina Garcia, Enric Massutì and Antonio Esteban from the Instituto Español de Oceanografía (IEO) for providing the data for the analyses and to all participants of the MEDITS project for collecting these data. We gratefully acknowlege two anonymous reviewers for their comments, which helped improve the manuscript.

REFERENCESTop

Anonymous. 2017. MEDITS Handbook. Version n. 9, MEDITS Working Group: 106 pp.
http://www.sibm.it/MEDITS%202011/principaledownload.htm

Blanchard J.L., Coll M., Trenkel V.M., et al. 2010. Trend analysis of indicators: a comparison of recent changes in the status of marine ecosystems around the world. ICES JMS. 67: 732-744
https://doi.org/10.1093/icesjms/fsp282

Brind’Amour A., Rochet M.J., Ordines F., et al. 2016. Environmental drivers explain regional variation of changes in fish and invertebrate functional groups across the Mediterranean Sea from 1994 to 2012. MEPS 562: 19-35.
https://doi.org/10.3354/meps11912

Burnham K.P., Anderson D.R. 2002. Model Selection and Multimodel Inference. A Practical Information-Theoretic Approach. 2nd ed.. Springer, New York.

Cheung W.W.L., Watson R., Pauly D. 2013. Signature of ocean warming in global fisheries catch. Nature 497: 365-368.
https://doi.org/10.1038/nature12156

Cury P.M., Christensen V. 2005. Quantitative ecosystem indicators for fisheries management. ICES J. Mar. Sci. 62: 307-310.
https://doi.org/10.1016/j.icesjms.2005.02.003

Dupont C., Barsoumian S., Moreira G., et al. 2018. Article 16 Technical Assessment of the MSFD 2015 reporting on Programme of Measures. Italy Report.

European Commission. 2010. EU Commission Decision of 1st September 2010 on criteria and methodological standards on good environmental status of marine waters (2010/477/EU). Off. J. Eur. Union, L232: 12-24.

European Union. 2013. Regulation (EU) No 1380/2013 of the European Parliament and of the Council of 11 December 2013 on the Common Fisheries Policy.

Fulton E.A., Smith A.D.M., Punt A.E. 2005. Which ecological indicators can robustly detect effects of fishing? ICES J. Mar. Sci. 62: 540-551.
https://doi.org/10.1016/j.icesjms.2004.12.012

Greenstreet S.P.R., Rossberg A.G., Fox C.J., et al. 2012a. Demersal fish biodiversity: species-level indicators and trends-based targets for the Marine Strategy Framework Directive. ICES J. Mar. Sci. 69: 1789-1801.
https://doi.org/10.1093/icesjms/fss148

Greenstreet S.P.R., Fraser H.M., Rogers S.I., et al. 2012b. Redundancy in metrics describing the composition, structure, and functioning of the North Sea demersal fish community. ICES J. Mar. Sci. 69: 8-22.
https://doi.org/10.1093/icesjms/fsr188

Harrell F.E. Jr. 2006. The Hmisc package.
https://cran.r-project.org/web/packages/Hmisc/Hmisc.pdf

Holmes E., Ward E., Wills K. 2018. MARSS: multivariate autoregressive state-space modeling. R package version 3.10.10.
http://CRAN.R-project.org/package=MARSS

Houle J.E., Farnsworth K.D., Rossberg A.G., et al. 2012. Assessing the sensitivity and specificity of fish community indicators to management action. Can J Fish Aquat Sci. 69: 1065-1079.
https://doi.org/10.1139/f2012-044

ICES. 2012. Report of the Working Group on the Ecosystem Effects of Fishing Activities (WGECO), 11-18 April, Copenhagen, Denmark. ICES CM 2012/ACOM:26. 192 pp.

ICES. 2014. Report of the Working Group on the Ecosystem Effects of Fishing Activties (WGECO), 8-15 April 2014, Copenhagen, Denmark. ICES CM 2014/ACOM:26. 174 pp.

ICES. 2018. Report of the Working Group on Ecosystem Effects of Fishing Activities (WGECO), 12-19 April 2018, San Pedro del Pinatar, Spain. ICES CM 2018/ACOM:27. 69 pp.

Kaiser H.F. 1974. An index of factor simplicity. Psychometrika 39: 31-36.
https://doi.org/10.1007/BF02291575

Lê S., Josse J., Husson F. 2008. FactoMineR: An R Package for Multivariate Analysis. J. Stat. Softw. 25: 1-18.
https://doi.org/10.18637/jss.v025.i01

Lynam C.P., Rossberg A.G. 2017. New univariate characterization of fish community size structure improves precision beyond the Large Fish Indicator. arXiv:1707.06569: 1-7.
http://arxiv.org/abs/1707.06569

McLachlan G., Peel D. 2000. Finite Mixture Models. John Wiley & Sons. Hoboken, New Jersey.
https://doi.org/10.1002/0471721182

Methratta E.T., Link J.S. 2006. Evaluation of quantitative indicators for marine fish communities. Ecol. Ind. 6: 575-588.
https://doi.org/10.1016/j.ecolind.2005.08.022

OSPAR. 2017. Intermediate Assessment. Fish and food web. Last accessed: 12/07/2019.
https://oap.ospar.org/en/ospar-assessments/intermediate-assessment-2017/biodiversity-status/fish-and-food-webs/

Pielou E. 1966. The measurement of diversity in different types of biological collections. J. Theor. Biol. 13: 131-144.
https://doi.org/10.1016/0022-5193(66)90013-0

R Core Team. 2019. R: A language and environment for statistical computing. R Foundation for Statistical Computing, Vienna, Austria.
https://www.R-project.org/

Rochet M.J., Trenkel V.M. 2003. Which community indicators can measure the impact of fishing? A review and proposals. Can. J. Fish. Aquat. Sci. 60: 86-99.
https://doi.org/10.1139/f02-164

Rochet M.J., Trenkel V.M., Carpentier A., et al. 2010. Do changes in environmental and fishing pressures impact marine communities? An empirical assessment. J. Appl. Ecol. 47: 741-750.
https://doi.org/10.1111/j.1365-2664.2010.01841.x

Shephard S., Reid D.G., Greenstreet S.P.R. 2011. Interpreting the large fish indicator for the Celtic Sea. ICES J. Mar. Sci. 68: 1963-1972.
https://doi.org/10.1093/icesjms/fsr114

Shin Y.-J., Shannon L.J. 2010. Using indicators for evaluating, comparing, and communicating the ecological status of exploited marine ecosystems. 1. The IndiSeas project. ICES J. Mar. Sci. 67: 686-691.
https://doi.org/10.1093/icesjms/fsp273

Shin Y-J., Rochet M.J., Jennings S., et al. 2005. Using size-based indicators to evaluate the ecosystem effects of fishing. ICES J. Mar. Sci. 62: 384-396.
https://doi.org/10.1016/j.icesjms.2005.01.004

Shin Y.J., Houle J.E., Akoglu E., et al. 2018. The specificity of marine ecological indicators to fishing in the face of environmental change: A multi-model evaluation. Ecol. Ind. 89: 317-326.
https://doi.org/10.1016/j.ecolind.2018.01.010

Souplet A. 1996. Calculation of abundance indices and length frequencies in the MEDITS survey. In: Bertrand J.A. et al. (eds), Campagne internationale du chalutage démersal en Méditerraneé. Campagne 1995. EU Final Report, Vol. III.

Spedicato M.T., Bitetto I., Carbonara P., et al. 2014. Methodological support for exploring the LFI: “proportion of large fish - indicator” properties for interpreting the changes in the fish community of the South Adriatic area (GSA18). In: First MedSuit Regional Workshop on indicators and targets to ensure GES of commercially exploited marine populations in the GFCM area. FAO HQs, Rome, Italy, 6-7 November 2014.

Spedicato M.T., Massutí E., Merigot B., et al. 2019. The MEDITS trawl survey specifications in an ecosystem approach to fishery management. Sci. Mar. 83S1: 9-20.
https://doi.org/10.3989/scimar.04915.11X

Vasilakopoulos P., Raitsos D.E., Tzanatos E. 2017. Resilience and regime shifts in a marine biodiversity hotspot. Sci. Rep. 7: 13647.
https://doi.org/10.1038/s41598-017-13852-9

Winemiller K.O., Rose K.A. 1992. Patterns of life-history diversification in North American fishes: implications for population regulation. Can. J. Fish. Aquat. Sci. 49: 2196-2218.
https://doi.org/10.1139/f92-242

Yong A.G., Pearse S. 2013. A Beginner’s Guide to Factor Analysis: Focusing on Exploratory Factor Analysis. TQMP 9: 79-94.
https://doi.org/10.20982/tqmp.09.2.p079

Zuur A.F., Tuck I.D., Bailey N. 2003. Dynamic factor analysis to estimate common trends in fisheries time series. Can. J. Fish. Aquat. Sci. 60: 542-552.
https://doi.org/10.1139/f03-030

SUPPLEMENTARY MATERIAL

The following supplementary material is available through the online version of this article and at the following link:
http://scimar.icm.csic.es/scimar/supplm/sm05015esm.pdf

Table S1. – List of the 23 species included in the analysis for the calculation of community indicators

Table S2. – Kaiser-Meyer-Olkin test: measure of sampling adequacy by GSA. In bold the KMO criterion values <0.5.

Table S3. – Lifespan by species used to estimate the average lifespan of the fish community considered.

Table S4. – Mean temperature of the catch, as found in Cheung et al. (2013), related to the species under investigation. For deep-water species (indicated by *), SST is not considered a reliable indicator of temperature preference.

Table S5. – Correlations (Pearson) among the indicators by GSAs. In bold, the significant correlations (P<0.05) are reported.

Table S6. – Trends of biomass (kg km^-2) by GSA/species over the time series (1999-2015) represented. Species not considered in areas are indicated as NA.

Table S7. – Trends of abundance (N km^-2) by GSA/species over the time series (1999-2015) represented. Species not considered in areas are indicated as NA.

Table S8. – Best models fitted by DFA for TyL and MML. The characteristics of the best models are described: the observation error covariance matrix R; the list of covariates included (none, one, two or up to three) and the Akaike information criteria corrected for small sample size (AICc), used to identify the best model for the indicators TyL and MML. The best model for both indicators, shown in bold font, has diagonal and equal matrix and no covariates.

Table S9. - Correlations between covariates. The significant correlations (p<0.05) are indicated in bold.

Fig. S1. – Residuals histograms of the best DFA fit for TyL and MML broken down by normal mixed distributions through a scale mixture model.

Fig. S2 – TyL time series by GSA.

Fig. S3. – MML time series by GSA.

Fig. S4. – Trend in biomass (kg km^–2) by species. The index is related to GSAs 1, 7, 8, 9, 10, 11, 15, 16 and 25.

Fig. S5. – Trend in abundance (N km^–2) by species. The index is related to GSAs 5, 6, 7, 8, 9, 11, 17, 18 and 19.

Fig. S6. – Fitting of the best model selected by DFA for TyL. The blue lines represent TyL time series for each GSA. The smoothed black line represents the common trend for TyL.

Fig. S7. – Fitting of the best model selected by DFA for MML. The blue lines represent MML time series for each GSA. The smoothed black line represents the common trend for MML.

Fig. S8. – Correlation between the four covariates used in the analysis: fleet capacity (GTxKwxN of vessels), NAO, Med anomaly and SST (°C). Correlation values are reported for each plot.

Fig. S9. – Temporal patterns of the four covariates used in the analysis: fleet capacity (GT x Kw x number of vessels), NAO, Med anomaly and SST (°C).