Modelling spatio-temporal patterns of fish community size structure across the northern Mediterranean Sea: an analysis combining MEDITS survey data with environmental and anthropogenic drivers Mediterranean demersal resources and ecosystems: 25 years of

Summary: The state of marine systems subject to natural or anthropogenic impacts can be generally summarized by suites of ecological indicators carefully selected to avoid redundancy. Length-based indicators capture the status of fish community structure, fulfilling the Marine Strategy Framework Directive (MSFD) requirement for Descriptor 3 (status of commercial fish species). Although the MSFD recommends the development of regional indicators, a comparison among alternative length-based indicators is so far missing for the Mediterranean Sea. Using principal component analysis and dynamic factor analysis, we identified the most effective subset of length-based indicators, whether or not based on maximum length. Indicator trends and time series of fishing effort and environmental variables are also compared in order to highlight the individual and combined capability of indicators to track system changes across geographical sub-areas. Two indicators, typical length and mean maximum length, constitute the smallest set of non-redundant indicators, capturing together 87.45% of variability. Only in combination can these indicators disentangle changes in the fish community composition from modifications of size structure. Our study supports the inclusion of typical length among the regional MSFD Descriptor 3 indicators for the Medi- terranean Sea. Finally, we show dissimilarity between the western and eastern-central Mediterranean, suggesting that there are sub-regional differences in stressors and community responses.


INTRODUCTION
Operational ecological indicators are generally used for summarizing the status of marine communities and ecosystems in a comprehensive and accessible way (Cury andChristensen 2005, Shin andShannon 2010). Indicators of anthropogenic pressure are useful for monitoring the impact of fishing on target species, communities (Greenstreet et al. 2012b) and ecosystems (e.g. Shin and Shannon 2010), thus helping to assess the effectiveness of management measures (e.g. Fulton et al. 2005, Houle et al. 2012. Ecological and impact indicators are crucial in the design of multiannual management plans, which, according to the Common Fishery Policy (EU Reg. 1380/2013, 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 (European Commission 2010).
Because no single indicator can capture the diversity of dynamics and processes within a system, suites of indicators are indeed required (Greenstreet et al. 2012a), though the selection of indicators to retain should be carefully considered: indicators included should be sensitive to specific impacts (Houle et al. 2012) or be unambiguous to monitor single species, communities or ecosystems Trenkel 2003, Shin andShannon 2010). 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 (Rochet et al. 2010, Shin andShannon 2010). The selected suite should also avoid redundancy (Methratta andLink 2006, Greenstreet et al. 2012a) and ensure applicability across different geographical areas (Shin et al. 2018).
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 (Shin et al. 2018). Most studies have focused on North Atlantic fisheries (Methratta andLink 2006, Greenstreet et al. 2012a) or compared multiple systems (e.g. Shin and Shannon 2010), while only a few have assessed the situation for the Mediterranean fish communities (Rochet et al. 2010, Brind'Amour et al. 2016.
Within the MSFD, suites of indicators of GES are outlined for commercial fish (Descriptor 3) and food webs (Descriptor 4) (European Commission 2010). 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; Shephard et al. 2011, ICES 2014, in particular due to the required assumption about threshold level, which is area-specific. In the Mediterranean Sea, for example, Spedicato et al. (2014) 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 Modelado de patrones espacio-temporales de la estructura de tallas de la comunidad de peces a través del norte del mar Mediterráneo: un análisis combinando datos de campañas MEDITS y factores ambientales y antropogénicos Resumen: Generalmente, el estado de los sistemas marinos sujetos a impactos naturales o antropogénicos puede ser resumido mediante un conjunto de indicadores ecológicos, cuidadosamente seleccionados para evitar la redundancia. Los indicadores basados en la talla reflejan el estado de la estructura de la comunidad de peces, cumpliendo el requisito de la Directiva Marco de la Estrategia Marina (MSFD) para el Descriptor 3 (estado de las especies de peces comerciales). Si bien MSFD recomienda el desarrollo de indicadores regionales, en el Mar Mediterráneo no se ha hecho hasta ahora una comparación entre los distintos indicadores disponibles basados en la talla. Mediante el análisis de componentes principales y el análisis de factores dinámicos, identificamos el subconjunto más eficaz de indicadores basados en la talla, estén o no basados en la talla máxima. Las tendencias de los indicadores y las series temporales del esfuerzo de pesca y las variables ambientales también son comparadas para resaltar la capacidad individual y combinada de los indicadores para detectar los cambios del sistema a través de las subáreas geográficas. Dos indicadores, Longitud Típica (TyL) y Longitud Máxima Media (MML), constituyen el conjunto más pequeño de indicadores no redundantes, captando juntos el 87.45% de variabilidad. Solo si se combinan, estos indicadores pueden discernir entre los cambios en la composición de la comunidad de peces y las modificaciones de la estructura de tallas. Nuestro estudio respalda la inclusión de TyL entre los indicadores regionales del descriptor 3 de MSFD para el mar Mediterráneo. Finalmente, mostramos diferencias entre el Mediterráneo occidental y el Mediterráneo centraloriental que sugieren diferencias subregionales en cuanto a factores impactantes y las respuestas de la comunidad. which might provide similar results and avoid problematic assumptions (ICES 2014, Lynam and Rossberg 2017, ICES 2018. 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 (ICES 2014, Lynam and Rossberg 2017, ICES 2018. In this paper we estimated these six indicators for 17 geographical sub-areas (GSAs) in the northern Mediterranean Sea, using data of lengthfrequency 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.

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]) following a random depth-stratified sampling and using a gear with a common design (Spedicato et al. 2019). 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; Anonymous 2017). 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 ca-pacity (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 (CFFR 2019) 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 lengthfrequency distributions (N km -2 ) according to Souplet (1996). 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.

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 1974). 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 2013) 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. 2003). 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: where y t is a N×1 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, a t 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 2002), following Zuur Table 1. -Indicators of size structure at community level. 1 , The threshold used for all the GSAs was 30 cm total length (TL) i.e. the value that returned the highest R 2 in the fitting of the 6 th degree polynomial in 8 out of 17 GSAs (for 6 GSAs it was 40 cm TL and for 3 was 35 cm TL); 2 , In this study g was set = 45 cm TL, which is the 25 th 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. ∑ = > LFI t y t y t ( ) ( ) ( ) l l lbig lbig threshold 1 ; y l (t) catch per length class l; y(t) total catch (measured species). Shephard et al. 2011OSPAR 2017 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. 2011

Evenness
Measure of the equitability in relative abundance among the length classes.
S is the number of lengths in the fish assemblage; p i is the proportion of individuals in the i-th length class. Pielou 1966 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. 2012b).

Mean maximum length
Mean maximum length in the community.
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 2012
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 2017, Lynam and Rossberg 2017) caused by increasing fishing pressure (Shin et al. 2005).

Large species indicator
Biomass proportion of large species in a community.

Shephard et al. 2011
A decline in LSI should indicate a decrease in the biomass of the predator species in the fish community (Shephard et al. 2011).
Mean weight Mean weight in the community.

Rochet and Trenkel 2003
A decreasing trend could be a signal of removal of larger individuals Trenkel 2003). et al. (2003) and Holmes et al. (2018).
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 2000) was applied (Supplementary Material Section 5).
The investigation of the canonical correlations between detected common trends and explanatory variables followed Zuur et al. (2003), 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  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 2019). The FactoMineR R library (Lê et al. 2008) was used for PCA, the Hmisc R library (Harrell et al. 2006) for correlation analyses and the MARSS R library (Holmes et al. 2018) for DFA.

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).
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).
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.

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.
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.

DISCUSSION
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) and by the ICES Working Group on Ecosystem Effects of Fishing Activities (WGECO) (ICES 2014(ICES , 2018. 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. 2012b) and ordination techniques such as non-metric multidimensional scaling (Blanchard et al. 2010) 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). The selected models were characterized by only one common trend, while Zuur et al. (2003) 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), 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) 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; . 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. 2017). 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 1992). 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 2017), 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 2017).
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 2017).
A complementary indicator that can account for changes in size structure has traditionally been identified in the LFI (European Commission 2010, OSPAR 2017), 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. 2011, Spedicato et al. 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), 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 2014, Lynam and Rossberg 2017, ICES 2018. 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), 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 andRossberg 2017, OSPAR 2017), 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) 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 2018). In the Mediterranean, these values do not yet seem to have been determined (Dupont et al. 2018), and future studies should therefore focus on this issue at different spatial scales.

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  , related to the species under investigation. For deepwater species (indicated by *), SST is not considered a reliable indicator of temperature preference. Table S5.

Adequacy of data to the analysis
According to Kaiser-Meyer-Olkin test, the available data were considered relatively adequate to the analysis (Table S2). Globally, GSAs 2, 17 and 22+23 showed the lowest KMO-Criterion, indicating widespread variance that could be scarcely interpreted as common variability. For all the GSAs the Bartlett test returned a p<0.05, thus we always rejected the null hypothesis that the correlation matrix was equal to the identity matrix. The results of the tests carried out ensure the validity of the hypothesis of the PCA allowing to perform the analysis on the available data.

Temperature preference by species
This section reports a table containing the available information of the temperature preference for the species under consideration. The mean temperature of the catch is obtained from , where available. For selected species, information from the online source Aquamaps (Kaschner et al. 2019;http://www.aquamaps.org) or from other sources was gathered.
The species with preference temperature higher than or equal to 17°C were considered thermophile.

Normality of residuals in DFA
The Shapiro-Wilk test detected deviations from normal distributions in the residuals of the DFA best fits for TyL and MML; this was investigated by scale finite mixture model. This analysis indicated that for TyL the 90% of residuals are represented by a normal distribution with mean -0.17 and standard deviation 0.69, while the second normal component (10% of residuals) is characterized by mean of 1.72 and standard deviation 0.77. For MML the residuals are divided in 96% due to a normal distribution with mean -0.13 and standard deviation 0.81 and another normal distribution (4% of residuals) with mean -1.82 and standard deviation 0.41 (Fig. S1).

Correlations among the indicators by GSAs
Before carrying out the PCA, the hypotheses needed to apply any factor analysis were verified. The correlations among the six indicators for each GSAs were statistically tested using Pearson's coefficient and the results are reported in Table S4. Time series Figure S2 shows the trend of the TyL for each GSA; Figure S3 the MML by GSA. Figure S4 reports the biomass indices (summing up the different GSAs) by species related to the GSAs where an increasing trend was found in the Typical length (Fig. S2), namely in GSAs 1, 10, 11, 15, 16, 25, 6, 7, 8 and 9. This Figure highlights the species responsible of the trend in the TyL, focusing on the areas where a stronger relationship between the indicator and the biomass (which the definition of Typical length is based on) was found. Figure S5 shows the abundance (density) indices by species, summing up the GSAs where a decreasing trend in the MML was observed (Fig. S3). Figure 5S has the aim to highlight the species responsible of the trend in the MML, focusing on the areas where a stronger relationship between the indicator and the density (which the definition of MML is based on) was found.

Z. faber
Length-based indicators of Mediterranean fish community • S11 SCI. MAR. 83S1, December 2019, S1-S14. ISSN-L 0214-8358 The number of common trends detected (M) is equal to 1 in all models selected.