A Novel Index of the Performance of Mytilus galloprovincialis to Improve Commercial Exploitation in Aquaculture

Estuarine environments are highly heterogeneous habitats where numerous organisms interact with each other. Aquaculture systems encompass such interactions, and the eventual yields depend on how the cultivated species respond to the environmental heterogeneity. Marine mussels are calcifying organisms that rely on calcium carbonate shell and byssus filaments secreted during their lifetime to protect soft vital tissues against aggressive (abiotic and biotic) environments. Nevertheless, these protective structures can be energetically costly, depending on the environment, thus affecting the energy allocation patterns in the organism. Consequently, other important fitness parameters for the aquaculture industry, such as soft tissue condition and market value, may also be affected. Here, we present a spatial and temporal analysis of the protective and soft tissues responses in the mussel Mytilus galloprovincialis with the aim of obtaining a better understanding of the inter-location variability in individual’s performance and the corresponding main environmental drivers. Local environment with regard to food availability and hydrodynamics impact very distinctly mussel tenacity and shell thickening, e.g., higher-energy environments at the outer exposed zones caused individuals secreted stronger byssus and thicker shells. By contrast, the state of soft tissues (condition index) was similar between very distinct and distant outer exposed and inner sheltered zones due to differences in both environmental drivers. A novel, intuitive ecological index that measures the impact of both protective and soft tissues was developed and is discussed in relation to cultivation timing. The data presented provides the basis for understanding the natural variability in energetic requirements for different vital tissues in bivalve mollusk that support survival and growth. We report the scientific basis for management actions aimed to shortening the cultivation cycle in the aquaculture sector. These actions are based on the combined use of the PROFIT index and other aquaculture practices (e.g., modifying density). PROFIT helps to identify when the quality of the product, understood as PROtection and FITness strategies, would be best suited for the market.


INTRODUCTION
The morphology and performance of bivalve mollusks in general, and marine mussels in particular, are known to vary in relation to key environmental factors for different species and habitats (Bergström and Lindegarth, 2016;Kroeker et al., 2016;Telesca et al., 2018). Like many other marine invertebrates, mussels secrete exoskeletons that provide a physical barrier that protects against abiotic and biotic aggressors (Burnett and Belk, 2018). Calcified shells therefore play a primary ecological role, reducing predation impacts on mussels, protecting them from intense wave action and providing mechanical strength for natural clustering in beds, cultivation ropes or matrices (Briones et al., 2014 and references therein).
As well as producing calcified shells, many marine bivalves (especially those in the family Mytilidae) secrete byssus, a non-calcified extracellular tissue. This attachment-related tissue confers the individuals the ability to maintain a relatively sessile life mode and hence, contributes to their successful cultivation around the world. Byssus filaments must be replaced continuously as the threads decay over time (approximately 4-6 weeks; Carrington, 2002). Secretion is thus dynamic rather than static and is also affected by various environmental parameters, seasonality and the endogenous rhythms of the organisms (Carrington, 2002;Zardi et al., 2006). As animals can modify energy allocation to protective tissues, energetic trade-offs between soft and hard (byssus filaments and shell) tissues production can occasionally occur (Carrington, 2002;Zardi et al., 2007;Babarro and Carrington, 2011). Therefore, for a comprehensive analysis of the fitness of an organism in a particular environment, the energetic requirements for the secretion of these protective structures must be considered. Byssus and shell secretion may represent up to 8-15% and 25-50% respectively of the total energy expenditure in mussels (Hawkins and Bayne, 1985;Gardner and Thomas, 1987;Steffani and Branch, 2003).
The shell represents a barrier that protects soft tissues being influenced by multiple environmental and endogenous factors and is generally linked to specific life-habitat and ecological niches of the organisms (Zuschin and Stanton, 2001;Zuschin et al., 2003). Shell thickness has recently been described as a proxy for the health status of shellfish as well as an indicator of environmental change, including climate change (Martinez et al., 2018). The aforementioned study reported shell thickness as "a new reliable indicator when informing the Marine Strategy Framework Directive (MSFD) (2008/56/CE), " in order to achieve "good ecological status (GES)" in European Marine Waters by 2020 (see Martinez et al., 2018). A clear link between shell thickness and strength has been established previously (Nagarajan et al., 2006) although other factors like shell shape, mineralogy local environment or genotypes may also impact shell performance (Penney et al., 2007;Fitzer et al., 2014). Shell fragility is not only important for specimens living in highenergy environments (see Kidwell and Bosence, 1991 for a review; Beadman et al., 2003;Steffani and Branch, 2003;Babarro and Carrington, 2011), but also for those inhabiting low-energy habitats where biogenic interactions can lead to shell breakage (Powell et al., 1989). Overall, shell tissues have been used to explore the intraspecific variations in different environments and its implications in aquaculture activities have been considered (Lewis and Magnuson, 1999;Beadman et al., 2003;Nagarajan et al., 2006).
Similarly, byssus failure, understood as lack of primary physiological response to attach optimally on substrates, has been used to explain dislodgement and biomass loss in aquaculture (Babarro et al., 2019). According to this study, the local environment in terms of hydrodynamics and food availability helped to minimize dislodgement from ropes through optimal byssus attachment strength. In the specific case of farmed mussels, operational mechanical procedures such as separation of clusters, removal of the byssus and cleaning attached fouling from the shell may lead to loss of marketable product due to shell breakage. Some examples of these types of actions affect approximately 5-15% of the whole annual production (G. Sarà, personal communication in Martinez et al., 2018). Cracked mussels are unwanted products representing production loss, and factors affecting the resistance of the shell to breakage other than ecological aspects are therefore also of great concern in the aquaculture sector (Penney et al., 2007).
Coastal embayments (rias) in areas characterized by upwelling systems, such as the Atlantic Galician coast (NW Iberian Peninsula), are heterogeneous ecosystems due to the mixing of salt and fresh water. They are also very productive in relation to nutrient and phytoplankton availability, which affect both farmed and wild bivalve species. These aspects are very important for rope-grown mussel cultivation. The heterogeneous environments may significantly affect the energetics and allocation patterns for investment in protective and vital tissues, with significant consequences for survival, biomass maintenance on cultivation structures and yield at harvest. Furthermore, the outcomes of this investment may also be affected by growth. As mussel growth involves continued byssus secretion, increased shell area and thickening, juvenile mussels will have thinner and weaker shells that can be more easily attacked and broken by predatory or mechanical actions; by contrast, adults will attach less strongly to substrates (e.g., littoral coastline in intertidal environments) as a consequence of lower byssus tenacity (Babarro and Carrington, 2013). The energetic expenditure may significantly affect the performance of individual mussels as the trade-off of energy allocation to byssus secretion in relation to growth and reproduction may be critical for individual fitness (Carrington et al., 2015). The ratio between soft tissue weight and shell weight (i.e., the condition index) may represent a proxy for product quality (Watanabe and Katayama, 2010), which is of great interest for harvesting purposes, e.g., to indicate the market value (Martinez et al., 2018). The physiological implications of mussel size on shell thickening or byssus attachment strength, together with the ecological implications driven by site-specific heterogeneity, will therefore play a key role in mussel performance in a widely exploited area such as the Ría de Arousa (NW Iberian Peninsula). To our knowledge, this study represents the first attempt to investigate how both byssus and shell secretion are affected by different natural environments in aquaculture systems, and how the potential energetic trade-offs affect mussel fitness, survival and growth performance. Field responses were combined in a novel and intuitive eco-physiological index that could be used to identify how the local environment may determine aspects that matter for aquaculture activities, i.e., biomass maintenance through attachment, structural integrity of the shell and yield or market value, and also when these characteristics are optimal for harvest within the usual cultivation period. Therefore, hypotheses planned in this survey can be summarized as follows: (i) local environment within a single embayment like Ría de Arousa (NW Iberian Peninsula) would modify significantly the energy allocated to distinct hard and soft tissues that matter for the life cycle of the mussels and, (ii) apart from individual size effect associated to organism's growth, food availability resources and hydrodynamic profiles, e.g., wave height together with salinity variability during the cultivation cycle would explain differences in mussels performance in the aquaculture system. In accordance with the importance of these protective and soft tissue structures in relation to aquaculture practices and markets, an integral index of mussel eco-physiology is proposed for evaluating mussel status with the aim of minimizing losses in biomass and better practices.

Description of the Study Area
The coastline of Galicia (NW Iberian Peninsula) represents the northern limit of the North Atlantic Upwelling System. The area is one of the world's major upwelling zones and is therefore highly productive. The northerly component of shelf wind stress on the Galician coast drives the surface layer offshore, inducing the cooler and nutrient-rich sub-surface waters close to the coast to rise the surface (Arístegui et al., 2009;Álvarez-Salgado et al., 2010). These oceanographic characteristics are modulated by the seasonal cycle of the wind direction, which is dominated by northeast winds between May and October (Wooster et al., 1976;Pardo et al., 2011). The highest mussel growth rates worldwide have been reported in the four large coastal embayments on the Galician coast (Pérez-Camacho et al., 1995;Fernández Reiriz et al., 1996), collectively known as the Rías Baixas. The Ría de Arousa is the largest of these embayments, which occupy an area of 245 km 2 (Figure 1). Mussels are cultivated in the area on floating rafts (20 m × 25 m), each of which supports 500 ropes each of length 12 m. The annual mussel production in Galicia is approximately 250,000 t, representing around 12% of the worldwide production (FAO, 2018).

Experimental Design and Sampling
Four different mussel farm sites were chosen for study in the Ría de Arousa (NW Spain; Figure 1), covering the outer exposed vs. inner sheltered spatial gradient (OuN and OuS at the northern and southern edges of the outer mouth of the ría, MidW in the middle section and InC in the innermost zone). Standard handling and cultivation techniques were used. Juvenile mussels of shell length approximately 20 mm were collected for use as seed from an intertidal environment on the Galician coastline. In all locations, the cultures were initiated with 900 mussels per m of rope at the first stage, from seeding to thinning out, and the density at the beginning of the second phase, from thinning out to harvest was 800 mussels per m of rope. Special caution was taken to ensure a similar seeding density in all ropes, as any differences could significantly affect mussel growth and physiological performance (Fréchette et al., 1996;Cubillo et al., 2012). The culture ropes were suspended on the rafts between February 2015 and August 2016 to encompass the two cultivation phases. Six ropes each of length 12 m (replicates), separated by a distance of 1 m, were hung from the aft of each raft (i.e., opposite to where the mussel raft is anchored with an iron chain, at the bow). Sampling was carried out approximately once a month. At each sampling time, 60 individual mussels (n = 240 for all cultivation sites) were dislodged from the mussel clusters along the length of the rope for estimation of byssal attachment force and tenacity (see below). The mussels used in the dislodgement tests occupied the outer positions of the cluster along the rope, as this provided easier access to each individual. Mussels were sampled at different depths along the rope (1, 3, 5, 7, 9, and 11 m) (n = 10 mussels per depth) with the aim to monitor several mussel responses and obtain sufficient biological material along the rope. The mussels dislodged from different sections of the ropes were also used to determine the shell thickness index (see below). Despite the main goal of sampling at different depths was to collect enough biological material, the data was pooled into two main groups. This allowed us to collect sufficient biological material, but also test the general effects of depth. The selection 1-5 and 7-11 m was made based on previous results in the Galician mussel cultivation for the importance of depth (Fuentes et al., 2000) and the empirical knowledge of the farmers with regard to both water sections depth and mussel performance.
A second group of mussels was sampled at the same depths as before to measure tenacity. Entire sections of the rope (inner and outer sections of the cluster) of known length were extracted. The maximum length, height and width of the shell were measured as the anterior-posterior, dorsal-ventral and lateral axes, respectively, in a minimum of 200 individuals. These mussels were also used to determine the condition index.

Environmental Conditions
A network of 38 oceanographic stations has been maintained by a local government agency in Galicia (Instituto Tecnolóxico para o Control do Medio Mariño 1 ) since 1992 (11 of them at Ría de Arousa), providing weekly hydrographic data. Data on the water column from the four stations closest to the mussel experimental rafts in the Ría de Arousa was obtained for a total of 147 vertical profiles, between February 2015 and August 2016 (Figure 1). Vertical hydrographic profiles were obtained throughout the water column, from the sea surface to the bottom at the four stations, with an SBE 25 Sealogger CTD (Sea-Bird Scientific). Chlorophyll a (Chla) was sampled at three depth intervals (0-5 m; 5-10 m; and 10-15 m) according to the Lindhal technique (Sutherland et al., 1992). Chlorophyll a concentration was measured by spectrofluorometric analysis (Zapata et al., 1994).
Current data were obtained from the output of a numerical configuration of the Regional Ocean Modelling System (ROMS) model (Schepetkin and McWilliams, 2005), for the three southernmost Galician Rías: Vigo, Pontevedra and Arousa. This configuration was established by a system of nested configurations used to describe the two-way exchange between the rías and the coastal ocean. Circulation in the rías is measured at a spatial resolution of 1/450 • (about 190 m) and includes wind-induced currents, tides and estuarine circulation due to inflowing fresh water. The results were validated with Intecmar data 1 and ADCP sampling of currents in the Ribeira buoy, also managed by Intecmar. The configuration is forced by wind stress, extracted from the Meteogalicia weather service through the threaded web server mandeo.meteogalicia.es, with a resolution of 4 km. River inflow from the tributary rivers in the Rías was downloaded from the measurements of gauging stations closer to the mouths of the Umia and Ulla rivers distributed by Meteogalicia 2 .
Wave parameters were calculated from the forecast power spectra of the SWAN model, developed at Delft University of Technology. The model computes random, short-crested wind-generated waves on a curvilinear grid on the Galician coast and was forced by Meteogalicia operational WW3 model. The output of the model estimated with an hourly frequency was averaged daily.

Mussel Tenacity
The attachment strength of mussels aggregated in clusters was measured by connecting a single mussel to a spring scale (Digital Force Gauge DN431, 0.01 N resolution) with the aid of custom-made forceps (see Babarro and Comeau, 2014;Babarro et al., 2019 for details of the procedure). Care was taken to avoid disturbing neighboring mussels when dislodging individual specimens. Individuals immediately adjacent to those selected for dislodgement were not considered in the trials if they had interconnected byssus threads. Individual mussel attachment strength in the field can be overestimated as a consequence of the dense network of fibers that several individuals secrete within the cluster. However, this functional value is valid for comparisons as it represents the force required to pull an individual from the cultivation ropes as done by currents or predators (Babarro et al., 2019). Dislodgement measurements were made with wet mussels to prevent modification of the mechanical properties of the byssus. Attachment force (F) was normalized by mussel size in order to calculate tenacity (TEN, N m −2 ), as follows: TEN = F/AP, where AP is the projected area of the individuals pulled for dislodgement, approximately an ellipse obtained by the product of width and height values of shells (Bell and Gosline, 1997).

Shell Thickness Index
After dislodgement of mussels, the anterior-posterior (length), dorso-ventral (height) and lateral axis (width) of the shell were measured (n = 60 per rope) with the aid of Vernier calipers (±0.1 mm). Live weight of the whole individual and dry weight of the shells were also recorded. Shell thickness index (STI) values were obtained, according to the following formula: (1) where L, H and W are respectively the length, height and width of the shell (Freeman and Byers, 2006;Freeman et al., 2009). To measure shell weight (wt), mussels were sacrificed, the tissue was dissected, and the shells were patted dry with paper towels and weighed on a Sartorius digital balance (±0.01 mg). Organic residues were removed from the shell surfaces and the shells were dried in a muffle furnace at 40 • C for 24 h.

Shell Length Growth and Condition Index
The second group of mussels sampled was used to determine growth and condition index. Mussels of different sizes can coexist within clusters on the ropes. Therefore, once the individual size (shell length) of all mussels from the cluster was determined, and with the aim of establishing the total weight of tissues (soft and shell dry weights), mussels assigned to specific interquartile range, e.g., 25, 25-50%, 50, 50-75%, and 75% of the latter sizefrequency range of the whole cluster were considered.
Shell length growth curves were fitted to Gompertz models: L t = L ∞ (e (−e(−k(t−t ))) ) where L t is the shell length (mm) at time t (months), L ∞ the maximum size, k is the growth parameter indicator of the speed at which maximum size is attained, and t is the inflexion point of the curve (Ratkoskwy, 1990). The Gompertz model parameters were estimated by non-linear regression, using the Levenberge-Marquardt algorithm and least squares as loss function.
The condition index (CI) of the mussels was calculated as follows: CI = (DW tissue /DW shell ) × 100, where DW tissue is the dry weight of the soft tissue and DW shell is the dry shell weight (Freeman, 1974). Both shell and soft tissues were dried at 110 • C for 48 h before being weighed.

Statistical Analysis
Mussel tenacity (TEN) and shell thickness index (STI) values obtained for the whole mussel cultivation period were analyzed using allometric functions in relation to individual shell length, SL (TEN or STI = aSL b ). ANCOVA was used to explore the effect of both cultivation site and depth as independent factors, and mussel size was considered a co-variable for TEN and STI variability.
TEN and STI were standardized to 60 mm shell length in order to determine the temporal variability throughout the cultivation cycle, allowing for site comparisons and seasonal patterns. For this purpose, specific allometric functions (TEN or STI = aSL b ) obtained for each cultivation site and depth were used to standardize the size effect.
Gompertz parameters of mussel growth were compared using a set of pairwise contrast and F statistical tests, according to Peteiro et al. (2006). Condition index variability was examined by ANCOVA, to test the effect of both cultivation site and depth as independent factors and mussel size as a co-variable.
Considering the previously described individual responses, two distinct eco-physiological strategies were explored for inclusion in the new PROFIT (PROtective -FITness) index: (i) shell thickness index (STI) and byssus functionality (TEN) as primary protective responses to aggressive abiotic (hydrodynamic) and biotic (predatory actions, inter-specific competition) factors, and (ii) condition index (CI), which as an indicator of soft tissue state provides an idea of the potential market value. The PROFIT index thus integrates both the quality of the shell (fragility) and attachment strength (maintenance on substrate), together with the marketable product (meat yield). Responses of a standard 60 mm-SL mussel were used to determine the PROFIT index and considering specific percentage values that each metabolic response represents in the energetic expenditure rates of the whole organism, as follows: an average value of 10% for the energy fraction represented by byssus secretion (see Griffiths and King, 1979;Hawkins and Bayne, 1985), 40% for the energy fraction represented by the shell tissue thickening (Gardner and Thomas, 1987) and 50% for mean somatic and reserve tissues (Fuentes-Santos et al., 2019). PROFIT is expressed as a percentage, and each component of the index (TEN i , STI i , and CI i ) for a site i (OuN, OuS, MidW, and InC) was corrected to the mean value for the four experimental sites and the whole cultivation period, as follows: ANOVA was used to test the effects of both cultivation site and depth on the variability in the PROFIT index.
Normality and homogeneity of variances were examined by the Shapiro-Wilk W-test and Levene's test, respectively. All data were log-transformed when necessary, and rank transformation was applied when the heterogeneity persisted (Conover, 2012). Significant differences between experimental groups were identified by a posteriori tests (Tukey and Bonferroni). All analyses were performed using SPSS Statistics 23 (IBM España S.A., Madrid, Spain) and STATISTICA 7.0 software (TIBCO Software Inc., Palo Alto, CA, United States). All data are reported as means ± SD.

Environment
The sampling period, between March 2015 and August 2016, in the Ría de Arousa was characterized by alternating periods dominated by northerly and southerly winds. However, northerly winds prevailed during spring and summer (Table 1), while southerly winds clearly dominated in the winter of 2016, with a mean value of −832 m 3 km −1 h −1 and a daily maximum value of −3808 m 3 km −1 s −1 in January (Figure 2a). During the period dominated by southerly winds, the wave height increased greatly, especially in the outer areas of the Ría de Arousa, reaching 4 m at the end of December 2015 (Figure 2b). Propagation of the waves along the Ría de Arousa attenuated the wave height, which never exceeded 2 m, in the inner waters. The waves were lowest in summer, of average height less than 1 m ( Table 1).
The discharges from the Rivers Ulla and Umia were significant throughout the winter and spring of 2016 (Table 1) and occurred in three pulses that reached maximum continental inputs of around 400 m 3 s −1 in February (Figure 2a). However, accumulation of fresh water in the estuary, with minimum salinity values close to 18 psu in the surface waters, was observed in the interior of the Ría de Arousa in January 2016 (Figure 2e). The clear signal in the variability of salinity occurred during piling up of the offshore surface waters toward the coast as a result of the southern winds, as shown by the strong negative Iw values (Figure 2a). Similar continental inputs in February and April did not cause such low salinity values as they occurred during null or positive values of upwelling index and thus with a net offshore transport of surface waters. Despite these differences, the large freshwater inputs caused notable salinity gradients in the surface waters, with lowest values at InC and OuN sites and maximum salinity at OuS (Figure 2e). The variability in salinity observed in the upper layer at the Ría de Arousa was greatly reduced at depths between 6 and 12 m, in which the salinity ranged between 31 and 36 psu during the study period (Figure 2f). Both surface and subsurface waters were very homogeneous by the end of summer, with values around 35.5 psu (Figures 2e,f). Similarly, the temperature of the surface waters also varied greatly, ranging from 11 to 20 • C, in comparison to the temperature span of subsurface waters, which ranged between 12 and 19 • C (Figures 2c,d). The temperatures were characterized by an obvious seasonal cycle repeatedly broken by the entry of deeper colder waters during upwelling at the mouth of the estuary and less often, the presence of warmer waters coming from offshore under southerly wind domain. The arrival of cold water to the exterior sites during the upwelling episodes warmed the interior waters, which are also shallower in the area of the Ría de Arousa, especially in summer.
The food availability (FA), a proxy for available food supply for the mussels, was estimated as the product of the chlorophyll a (Chl-a) concentration and current velocity. The FA indicated strong short-term variability, with intense phytoplankton proliferation pulses at all experimental sites and a period where there was almost no food available for mussels, between December 2015 and January 2016 (Figures 2g,h and Table 1). The FA values were generally higher in the surface waters at the outer exposed northern OuN site where high current velocity and maximum Chl-a values coincided. The FA values were lowest in the southern outer site (OuS) in the Ría de Arousa. The values estimated for the area along the estuary between 6 and 12 m depth were lower than those found at the surface, but were more homogeneously distributed due to the higher concentration of chlorophyll in inland waters, as observed at the MidW and InC sites ( Table 1).

Mussel Tenacity (TEN)
The variability in mussel tenacity in relation to individual length in the selected experimental sites and for different sections along the ropes is shown in Figures 3A,B. Byssus tenacity was inversely related to the increase in mussel size according to negative allometric functions (see equations in Figure 3). This relationship is a consequence of the relatively lower magnitude of the increase in byssus attachment force relative to shell area during growth. Therefore, older and larger mussels are relatively weaker in terms of tenacity than juveniles (Figures 3A,B). ANCOVA applied to TEN variability reflected the size effect as a co-variable (p < 0.001) and also the impact of cultivation site (p < 0.001) and depth (p < 0.05) (Tables 2, 3). The TEN values were higher for the northern outer exposed site (OuN) followed by the southern outer site (OuS), and finally the weakest tenacity was observed for mussels in the middle and most inner sheltered sites (MidW-InC), with no differences between these sites (Figures 3A,B and Tables 2, 3). Cultivation depth had a minor (residual) impact on TEN, which was slightly lower (4.5%) in the deeper section (6-12 m) than in the upper section (0-6 m) of the water column at all sites; no effect was noted for the site × depth interaction term (Figures 3A,B and Table 2).
Spatial and temporal analyses of the standardized TEN (60 mm-SL) confirmed the main effect of the cultivation site (see Figure 3C with average TEN values for both depth intervals). This finding highlights the stronger byssus tenacity in mussels from the northern and southern outer exposed sites (OuN followed by OuS) than in mussels from the innermost locations, whereas seasonality was represented by greater TEN through the autumn-winter period than in the spring-summer period ( Figure 3C).

Shell Thickness Index (STI)
The variability in STI for the cultivated mussels is represented in Figure 4. Cultivation site was the only factor that significantly affected STI (Tables 2, 3), apart from the individual size effect, as illustrated by the allometric functions obtained (Figures 4A,B). Standardized STI (i.e., 60 mm-SL) is shown to represent the TABLE 1 | Variability of the environmental parameters temperature ( • C), salinity, chl-a (mg m −3 ), food availability (FA, mg m −2 s −1 ), wave height (m), upwelling index (10 3 m 3 km −1 s −1 ) and river input (m 3 s −1 ) in the different sites of cultivation and as a function of experimental time. 4.9 ± 3.5 6.6 ± 5.2 0.6 ± 0.5 0.8 ± 0.7 0.5 ± 0.1 Bold and italic numbers represent values for the water sections 0-6 and 6-12 m, respectively. Iw (bold) and river inputs are shown as single value for the whole experimental area (Ría de Arousa). temporal variability in the index ( Figure 4C) after removing the impact of mussel size. Two different aspects can be highlighted: the seasonal effect of the standardized STI, which was significantly higher (by approximately 15%) between autumn 2015 and winter 2016 (September-February) than in the summer period (May-August) for all cultivation sites ( Figure 4C). On the other hand, inter-site variability showed a significant effect of culture location on STI variability (p < 0.001) (Figure 4 and Tables 2, 3). Mussels from the most inner and sheltered site (InC) had much thinner shells than mussels from the other locations with the thickest shells being obtained for OuS mussels followed by OuN and MidW mussels (Figures 4A,B). No differences were observed in relation to depth of water ( Table 2).

Condition Index (CI) and Mussel Growth
The values of the condition index (CI), i.e., the ratio between soft and shell tissue weight, for the different cultivation sites and depths are shown in Figures 5A,B. The ANCOVA showed that CI was significantly influenced by both site and depth of cultivation (p < 0.01), together with the co-variable mussel size (p < 0.001) (Tables 2, 3). The main difference for the inter-site comparison corresponded to both outer sites (p < 0.01), the northern (OuN) and the southern (OuS) sites, with highest and lowest CI values, respectively (Figures 5A,B). All other comparisons between sites were not significant with the exception of the minor effect noted for the comparison OuN versus MidW (Table 3). Interestingly, the CI values were similar in very different sites in terms of hydrodynamics, food availability and salinity regime, like OuN and InC (i.e., in the outermost exposed and inner sheltered areas) (Figures 5A,B). Regarding cultivation depth, CI was significantly lower in deeper waters (6-12 m) than in the surface waters (0-6 m), approximately 10% higher on average for all sites, as the interaction term site x depth was not significant (Figures 5A,B and Table 2). Shell length growth curves for the mussels cultivated in the different experimental sites are shown in Figure 6. Values of the growth-related parameters of the Gompertz models used to analyze these curves are shown in Table 4. Significant differences between mussel populations referred to the maximum shell length obtained, i.e., L ∞ but also for the slope of growth vs. time regressions (-k) and the inflecion point of the curve (t ) ( Table 4). Asymptotic length (84-88 mm SL: considering the entire length of the ropes as the depth factor was not significant) was higher in mussels cultivated in the outer exposed northern site (OuN) FIGURE 3 | Byssus tenacity in the mussels cultivated at the four distinct locations in the Ría de Arousa, at depths of (A) 0-6 and (B) 6-12 m in relation to mussel size (SL, mm). Variability in byssus tenacity (as mean values for the entire length of the cultivation rope; moving average) in relation to (C) cultivation time once standardized to a common mussel size of 60 mm shell length. The dashed line represents the timing of the thinning out process usually carried out between seeding and harvest. than in mussels from the other sites, which did not differ significantly among them (67-73 mm SL) (Figure 6 and Table 4). This represented mean values of 14-22% greater shell length growth for OuN mussels relative to the other cultivation sites (Figure 6). Differences in the other two regression parameters (k and t ) are a consequence of the lengths obtained throughout the experimental period, highlighting the greater growth of mussels in the outer exposed OuN site (Figure 6). As water depth did Shell length (SL) of the mussels was used as co-variable. Values in bold are statistically significant (p < 0.05); ns: not significant. For analysis assumptions, see section "Materials and Methods."  not significantly affect the shell growth parameters, the plotted values in Figure 6 represent mean values for the entire length of the rope (0-12 m).

Integrating PROtective vs. FITness Responses: PROFIT Index
The PROFIT index was developed to illustrate the impact of local environmental conditions and how mussels respond with actions that are important for both survival (byssus tenacity and shell thickening) and yield (i.e., CI) with a simple ratio. The variability in this ecological indicator as a function of position along the Ría de Arousa is shown in Figure 7. ANOVA applied to the variability in PROFIT revealed a significant impact of the cultivation site, but no effect of depth or the interaction term (site × depth) ( Table 5). Then, the PROFIT values shown in Figure 7 corresponded to the entire length of the rope (0-12 m). In general, the PROFIT values were highest (p < 0.01) for mussels from the outer exposed northern site (OuN) throughout the cultivation period (Figure 7). Differences between mussels from the other locations (OuS-MidW-InC) and differences in relation to depth were not significant ( Table 5). Regarding the timing of harvest, for both outer and more exposed sites (OuN and OuS), PROFIT values were highest after 8 months of cultivation when the shell length reached approximately 60 mm (maximum PROFIT values of 33 and 12% for OuN and OuS, respetively). For the OuN site, where the mussels were largest (86 mm SL), maximun PROFIT values were only comparable to that reported at the end of the experiment (31%; after 18 months). For the OuS site, PROFIT values at the end of the experiment were still lower than those observed at 8 months (Figure 7). For the innermost sites, MidW and InC, the PROFIT values were highest (14-18%) at the end of the experimental period (18 months) (Figure 7), revealing a clear difference between outer and inner estuarine areas.

DISCUSSION
The aim of this study was to relate a number of field responses in the mussel Mytilus galloprovincialis, linking individual performance, understood as the capacity to withstand natural threats (flow currents, waves, physical disturbance by competition etc.), to condition index and the associated market value. For this purpose, byssus tenacity and shell thickening responses were considered proxies for the ability of mussels to attach to the substrate and resist hydrodynamic forces and predatory or competitive actions, respectively though the latter factor (biological pressure) would be of lower magnitude in the cultivation system as compared to rocky shore assemblages. As such actions may affect the commercial value by making the shell more fragile, the condition index was considered as proxy for the market value. The complementary analysis of the strategies used by the mussels during their life cycle enabled us to develop FIGURE 6 | Shell length growth curves for the cultivated mussels at the distinct experimental sites of the Ría de Arousa and considering the whole length of the cultivation rope (0-12 m). The dashed line represents the timing of the thinning out process usually carried out between seeding and harvest.
a novel, intuitive index integrating both energetic expenditure routes balancing protection/survival and fitness.

Protection and Fitness Along an Environmental Gradient
Differences in shell morphology, shape or thickness of bivalve mollusks are influenced by key environmental parameters such as food competition, substrate type, crowding, temperature, wave impact and predatory actions, among others (see Alunno-Bruscia et al., 2001;Beadman et al., 2003;Steffani and Branch, 2003;Valladares et al., 2010). In the present study, considering the spatial gradient between outer exposed and inner sheltered waters, we observed a significant decrease in shell thickness of mussels cultivated in the sheltered site (thinner and weaker shells). The same pattern was also observed for other protective tissue like byssus (weaker attachment). As well as being thinner, the shells of the mussels in the most inner and sheltered site (InC) were also more elongated, although shell height did not differ (see Supplementary Appendix A). Regarding the environment, the different patterns in shell shape were closely associated with the lower (two-fold) hydrodynamic values (wave height) in the most inner-sheltered site, more than any other parameter. Temperature and salinity can affect shell development  All parameters are statistically significant (p < 0.001). Distinct letters and numbers identify significant differences for growth parameters.
Frontiers in Marine Science | www.frontiersin.org FIGURE 7 | PROFIT index variability throughout the whole cultivation cycle (moving average) of the mussels at the four experimental locations of the Ría de Arousa. Values shown correspond to the entire length of the rope (0-12 m) because of the non-significant differences by water depth. The dashed line represents the timing of the thinning out process usually carried out between seeding and harvest. Since each component of the PROFIT (expressed as a percentage) was corrected to the mean value for the four experimental sites and the whole cultivation period, positive and higher values represent better performance of the mussels. (Kautsky et al., 1990;Nagarajan et al., 2006;Martinez et al., 2018;Telesca et al., 2018), but their impact, especially for temperature, can be considered minor in this study as the greatest differences between the most separated sites (inner InC and outer OuN-OuS) were reached in summer and winter periods respectively, with differences of 0.57 • C and 2.67 psu, for temperature and salinity, respectively. For the case of salinity, it is known that low values of this factor are accompanied by lower pH and CaCO 3 saturation states, reducing the availability of calcium ions and dissolved inorganic carbon for calcification (Miller et al., 2009). However, the manuscript also focused on the need and importance of integrated approaches to the new ecological PROFIT index reported here which integrates different energetic aspect that are important for the survival and fitness of organism. Since salinity variation is relatively high at the most inner and sheltered InC site which in turn must be accompanied by supplementary energetic cost for intracellular osmotic pressure adjustments (Neufeld and Wright, 1996), this factor must be considered further for testing the effect of these punctual salinity decreases. Furthermore, the hydrodynamics appear to be more important for explaining the differences in the protective structures since salinity and temperature were similar in the middle site (MidW) and the most exposed sites (OuN-OuS). Similar results were obtained for shell thickness, byssus tenacity and condition index in mussels from the MidW and inner InC sites, with comparable hydrodynamics. For most mussels in inner sheltered sites, thinner and weaker shells will increase shell fragility (see Penney et al., 2007). Both individual size (e.g., refugee size) and shell thickness are key factors regarding the vulnerability of the mussels as targets for predators (Leonard et al., 1999;Smith and Jennings, 2000;Zuschin and Stanton, 2001;Zuschin et al., 2003). Shell thickening is a well-known protective strategy in bivalves and is probably influenced by the evolutionary history of interactions between species (Freeman and Byers, 2006). Shell fragility is also important in aquaculture as the mechanical process of breaking the mussel clusters apart to separate individuals can lead to a significant loss of biomass (broken shells). Furthermore, factors such as shell appearance and integrity are also important, especially for the fresh market as broken shells can lead to rejection of the product for commercial purposes. Between 5 and 15% of the whole annual production may be lost due to broken shells being considered undesirable products (G. Sarà, personal communication, see Martinez et al., 2018). In Galicia, it is commonly known that mussels from outer/exposed sites are mainly destined for the fresh market, in which the appearance of the (more robust) shell is very important; mussels from the innermost sheltered sites are used in the canning industry or in processed products for which only the soft tissues are used. Consistently, in the present study, both strength of attachment and shell thickness were greater in outer and exposed sites (especially OuN) than in the inner areas. Future studies should address the biomass losses associated with shell failure or fragility as the spatial differences in both environmental factors and shell/byssus responses in the mussels from the Ría de Arousa are well known. Shell length growth of the mussels differed significantly depending on the location. Growth was higher in the most outer and exposed (northern) site (OuN) than in the other sites, particularly the inner-sheltered InC, in which growth in length was lowest at the end of the experimental period, regardless water depth. This finding is consistent with those of previous studies in the same embayment regarding the spatial patterns of growth, i.e., greater mussel growth and productivity for the northern outer site (OuN) (Pérez-Camacho et al., 2014). Similarly to growth parameters, STI was less affected by water depth, and only byssus TEN (residually) and CI (10%) decreased in deeper waters probably in association with lower availability of natural resources and the lower level of movement at depth than in surface waters. In contrast to shell-related parameters (thickening and growth), the variations in CI were very different. Highest and lowest CI values were obtained in the outer exposed northern (OuN) and southern (OuS) sites, respectively. Beside, no differences were noted between mussels grown at the OuN site and those cultivated at the inner cand calmer InC site throughout most of the study period, except in autumn 2015 and in the surface waters (0-6 m). Mussels growing in the southern exposed site (OuS) had to develop high-standard protective tissues (TEN and STI) as a consequence of the high-energy environment, e.g., hydrodynamics. However, the availability of natural food resources was lower in the southern site than in the northern site (with comparable hydrodynamics) and similar to that in the other inner, calmer locations. This energetic trade-off may have occurred at OuS where limited food resources available needed to be allocated to stronger byssus attachment and shell at the cost of soft tissue development, e.g., CI (see also Carrington, 2002;Babarro and Carrington, 2011). Following previous reasoning in regard to energy allocation patterns, mussels in the innermost sheltered site (InC) may have taken advantage of the less energy-demanding environment and in consequence, the energy not required for strong attachment or shell thickening may have been used for soft tissue growth.
Although the possibility that differences encountered in this study with regard to byssus tenacity and shell thickness index can partly reflect a genetic component in terms of divergence or selection cannot be completely ruled out without experimentation supporting that, it can be proposed for further research with complementary experimental design. High density of mussels in the rías and their spawning events in the previous spring-summer period are responsible of the seed attached along the rocky shore forming a relatively mixed genetic pool. Mussel seed used in our survey was obtained from the same intertidal rocky shore location of the Galician Coastline; consequently same origin was chosen to start mussel cultivation at the four different sites along the Ría de Arousa, which allow us to assume that variability shown in the mussel responses for the present study were due to phenotypic plasticity and adaptation patterns to key environmental parameters like food resources available, hydrodynamic, salinity etc. that characterized each zone of the ría. A detailed phylogeographic study to investigate genetic diversity in the Galician Rías showed that inter-population differences were half of that from other Iberian Atlantic samples using the same loci (Diz and Presa, 2009). Both large-scale (among five different Galician Rías) and short-scale (comparing inner and outer zones of single ría as performed in our own study) analyses showed a relatively weak divergence in genes which can be interpreted as not clear genetic differentiation (Diz and Presa, 2009).
The new ecological PROFIT index reported here integrates different energetic aspects that are important for the survival and fitness of organisms. This non-dimensional ratio may provide accurate information about the magnitude of certain energy allocation patterns that depend on widely heterogeneous environments in a given geographical area. PROFIT is a good indicator of the specific habitat conditions in terms of food availability, considering both natural seston load and hydrodynamics. Indeed, the northern exposed zone (OuN), according to food availability resources monitored, was identified as an environment where mussels may develop key metabolic processes that matter for the life cycle, with enough energetic resources to get high-standard protection through a strong byssus attachment and shell thickening, and high growth performance in relation to shell length and soft tissue production. Differences in the PROFIT index between the other experimental sites (OuS-MidW-InC) were of minor interest as was the effect of water depth. The fact that PROFIT values for mussels growing in the latter sites did not reach similar values in relation to commercial size and timing of cultivation period was a consequence of limiting energetic resources that could not meet the requirements for both protective (byssus and shell) and soft tissues equally. PROFIT index reported here can be easily updated in the near future with extra information about health of sampled mussels as well as environmental status of interest in mussel farming areas such as water pollution (D'Agata et al., 2014;Cappello et al., 2018), using the role of multi-marker approach to monitor pollution exposure or anthropogenic effects on immune system of the mussels (Giannetto et al., 2017;Caricato et al., 2019;Parrino et al., 2019).

Implications for Aquaculture
The most important finding for the aquaculture industry is that irrespective of the energy distribution between protective tissues and yield for a given environment, the PROFIT index will reflect both the importance of the spatial analysis and the optimal timing for harvesting. The timing is determined by the size of the mussel and the soft tissue yield (CI) as well as by the capacity of the mussels to remain attached to the cultivation rope (TEN) and to have a shell that can support the handling process without breakages (STI). The variability in TEN and STI in the sub-tidal mussels showed that juveniles were more strongly attached than adults (exponential decrease in TEN with size) whereas adults had stronger shells than juveniles (exponential increase in STI with size), as previously reported for intertidal individuals attach on rocky shore (see Babarro and Carrington, 2013). Thus, a general scheme can be considered, in which juveniles are more prone to being predated (thinner and weaker shells) and adults are more prone to being dislodged (weaker byssus tenacity) (see Supplementary Appendix B). The use of TEN and STI to calculate PROFIT enables integration of these often-overlooked characteristics, which vary throughout the life cycle of the mussel and are critical for maximizing the quality of the product depending on local environment.
The current perception of Spanish mussel industry is that market interests have shifted to smaller sizes of the cultivated product, as in France and Italy in recent years (see fluctuations in the listed price change for the period 2004-2012, Monfort, 2014). Accordingly, the percentage of large mussels obtained at harvest decreased significantly, from 24 to 8%, in the last 15 years, and the opposite was observed for the proportion of small sized mussels harvested, which increased from 51 to 68% (Aquaculture YearBook 3 ). Indeed, the minimum commercial size for mussels was 70 mm shell length a few decades ago, whereas the corresponding value for Spanish mussel cultivation is now 50 mm (Pérez-Camacho et al., 2013), which is similar to the size demanded by the market. In 2018, the European Commission endorsed the inclusion of small sized mussels in the Protected Designation of Origin (PDO) scheme used in Galicia for mussel certification with the aim of satisfying the European demands, while always preserving high quality.
The study findings show that the PROFIT index, which standardizes all aspects that matter for the survival and yield of the cultivated product, reached the maximum value in the outer sites after 8 months (winter), although with significant differences between the northern (OuN) and southern (OuS) sites according to food availability resources. The PROFIT value then remained lower, only reaching values close to the maximum at the end of the experimental period (18 months). Regarding mussel size, the maximum PROFIT value was reached at 60 mm shell length after 8 months, which is lower than the size above which a negative effect of high density on size/weight growth was reported for typical farming practices in Spain, 66 mm (Cubillo et al., 2012). This information can be used to offer a set of recommendations regarding the shortest time within which individuals in these locations can be harvested. New expectations and market tendencies can be taken into account, to prevent loss of market quality (CI variability) and to decrease the risk of dislodgement (based on byssus tenacity) during an extended cultivation period. By contrast, for mussels cultivated in inner and calmer sites, the PROFIT index reached the maximum value at the end of the cultivation period (18 months), possibly suggesting that the best compromise between optimal meat yield and strength of attachment and shell thickness is obtained much later than in outer and more expose sites due to lower availability of food resources. Consequently, depending on the cultivation site, harvesting could be managed according to the market interests for the specific size of mussels demanded and the appropriate cultivation timing, according to PROFIT, in order to reduce both total cultivation time and labor costs. Shortening the cultivation time by aiming for smaller mussels at harvest would allow for high culture densities and the thinning process to be omitted, which would represent a significant improvement in cultivated biomass and yield (Pérez-Camacho et al., 2013). This strategy would maximize growth and minimize dislodgement, given the limited intraspecific competition (Cubillo et al., 2012) and high tenacity (this study) of small mussels. Here, we complement this innovation strategy (Pérez-Camacho et al., 2013) with the underlying eco-physiological mechanisms to support any potential shift in the classical development of mussel raft cultivation, as a novel and integral management tool to valorize different production zones or exploitation methods.
PROFIT index can be used to optimize the timing of harvest and minimize dislodgements and consequently also minimize the negative effects of biomass on the sea floor.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AUTHOR CONTRIBUTIONS
JB and XP conceptualized and performed the experiment. ML helped in the field sampling organising trips and monitoring mussel responses on the raft. RF supported with data analysis, discussing, writing and editing the manuscript. JB wrote the first draft of the article. All authors added valuable comments on the writing the first draft and contributed to manuscript revisions and read and approved the submitted version.