Grouping behaviour impacts on the parasitic pressure and squamation of sharks.

The evolution of grouping behaviour involves a complex trade-off of benefits and costs. Among the latter, an increase in the risk of parasitic transmission is a well-documented phenomenon that has likely promoted the evolution of defensive mechanisms in aquatic vertebrates. Here, we explore the relationship between grouping behaviour, parasitic richness (∼parasitic pressure), and the evolution of potential defensive traits in the squamation of sharks through phylogenetic, standard and zero-inflation regression models. Our results demonstrate that sharks that frequently aggregate show increased parasitic pressure, which may constitute an agent of selection. Accordingly, their squamation is characterized by large-scale crown insertion angles and low-scale coverage, which are interpreted as traits that compromise parasite attachment and survival. These traits are less evident in regions of the body and ecological groups that are subjected to high abrasive stress or increased drag. Thus, the squamation of sharks responds to a compromise between various functions, where protective and hydrodynamic roles prevail over the rest (e.g. ectoparasitic defence and bioluminescence aiding). This work establishes a quantitative framework for inferring parasitic pressure and social interaction from squamation traits and provides an empirical basis from which to explore these phenomena through early vertebrate and chondrichthyan evolution.

Introduction

Grouping behaviour is a widespread phenomenon among aquatic vertebrates [1-3] that involves a complex trade-off of benefits and costs associated with intraspecific and interspecific interactions [4]. Possible drivers for the evolution of grouping behaviour in aquatic vertebrates include advantages on swimming energetics and migration [5], reduced vulnerability to predators due to dilution or confusion effects [6,7], and cooperative hunting [8], among others. However, the establishment of groups may also entail a series of disadvantages, such as a decrease in foraging efficiency [9] or an increased risk of parasitism [10-13]. The latter is particularly relevant in the case of monoxenic parasites (i.e. those with direct life cycles), which require only one host to complete their development and, therefore, the probability of successful transmission is enhanced by higher host densities [14].

Parasitic interactions are an essential component of aquatic communities [15,16]. They trigger the establishment of evolutionary ‘arms races’ between hosts and parasites and promote the evolution of defensive adaptations in the former [17], which may include physiological and immunological responses [18], behavioural strategies [10,19,20] and morphofunctional traits [21]. Most monoxenic parasites of primarily aquatic vertebrates are ectoparasites that predominantly anchor to the branchial or skin surfaces (i.e. Monogenea, Copepoda, Hirudinea and Thecostraca) [22,23]. These groups cause significant damage to the respiratory organs and epidermal structures, resulting in the exposure of internal tissues and facilitating the development of infections [22-24]. Hence, several authors have suggested that the squamation of sharks plays a primary role in preventing the settlement of ectoparasites [25,26].

The squamation of sharks is a complex anatomical structure consisting of minute dermal denticles (i.e. placoid scales) that display topological, ontogenetic, sexual and interspecific morphological variation highly correlated with different scale functions and ecological attributes. Consequently, a well-defined categorization of placoid scales into morphotypes and functional types has been established in the literature [21,26-28] (electronic supplementary material, figure S1): the gracile carinated scales present in active pelagic species are involved in drag reduction [29-31]; the robust densely packed scales with a rhomboidal or rounded shape that are common in all sharks in areas of the body that are subjected to high abrasive stress (i.e. the rostrum, mouth and leading edges of the fins), and covering most of the body in demersal sharks, are associated with protection against abrasion or stabilization of the leading edges of the fins [32]; the star-shaped scales with concave facets and needle-like scales accommodate photophores and permit the passage of the light emitted by those organs in bioluminescent species [33,34]; the hook- or thorn-like scales found on the skin of slow schooling sharks are presumed to play a role in preventing the settlement of ectoparasites and epibionts [21,26,27]; and the strongly carinated scales with well-developed lateral wings that cover the bodies of slow-moving species in open water, and the dorsolateral region of demersal species that live on soft substrates, are associated with more generalized functions [27] (see also [35,36] for more specific roles of placoid scales). Although hydrodynamic, abrasion resistance and bioluminescence-aiding functions are well proven in previous experimental and/or observational studies [26,33,37-41], empirical evidence to support a potential defensive role against ectoparasites is anecdotical [42-44] or even inconsistent ([21] and references therein).

Resolving the interrelations regarding grouping behaviour, (ecto)parasitic pressure and potential defensive morphological adaptations in the squamation of sharks is, therefore, critical to (1) achieve an improved understanding of the trade-offs and costs associated with the evolution of social behaviours and the drivers of parasitic interactions in sharks and, more generally, in aquatic vertebrates and (2) uncover potential morphofunctional signals in the squamation of sharks associated with grouping behaviour and the inherent risk of increased parasitic pressure, thus, providing a framework to indirectly infer and track these phenomena through chondrichthyan and early vertebrate evolution. Here, through the implementation phylogenetic, standard and zero-inflated regression models, we explore whether grouping behaviour in sharks is associated with higher parasitic pressure and assess whether this can be linked to specific squamation traits that may prevent the settlement and compromise the survival of ectoparasites, with direct implications in the design of antifouling biomimetic surfaces.

Material and methods

Database compilation

We generated a database on 213 shark species that included information on parasitic richness (PR) (i.e. the number of parasitic species in a given host species), the occurrence of grouping behaviour, and ecology (electronic supplementary material, table S1 and figures S2 and S3). Parasitic richness was calculated for each species based on records compiled by Pollerspöck & Straube [45], and only the main groups of monoxenic ectoparasites of sharks were considered (i.e. Monogenea, Copepoda, Hirudinea and Thecostraca). The shark species were categorized according to the frequency they present grouping behaviour, based on the information provided by Ebert et al. [46], the IUCN [47] and Froese & Pauly [48]. We classified the species into ‘taxa that show grouping behaviour infrequently’ (i.e. always found as solitary individuals or with sporadic records of aggregations) and ‘taxa that show grouping behaviour frequently’ (mostly recorded in aggregations) (GB1). In order to test the sensitivity of the model outcomes to uncertainty in this classification, we considered a second categorization in which the former category was split into ‘solitary taxa’ and ‘taxa that show grouping behaviour occasionally’ (GB2). The occurrence of grouping behaviour was considered to be unknown for species with poorly studied ethology and ecology (according to [46]), which were excluded from subsequent analyses. The shark species were further categorized into ecological groups (EG) following the classification system proposed by Ferrón & Botella [27] (i.e. strong-swimming pelagic species, schooling species of low to moderate speed, demersal species on rocky substrates and in caves, demersal species on sandy and muddy substrates, mesopelagic bioluminescent species, and slow species of the open water) and based on the ecological data reported by Ebert et al. [46], the IUCN [47], and Froese & Pauly [48]. Finally, we also recorded the number of published studies (NS) concerning each of the considered shark species based on the bibliographic references listed in Pollerspöck & Straube [45].

Shark squamation

Simultaneously, we studied the squamation of 63 shark specimens of 49 different species from the collections of the Museum für Naturkunde of Berlin (MfN, Germany) and the Australian Museum of Sydney (AM, Australia) (electronic supplementary material, table S2). This selection represented suitable coverage of the diversity of the lifestyles displayed by sharks and included representatives of all the ecological groups described by Ferrón & Botella [27]. We acquired skin samples from eight different topological positions (P) (figure 1) using dermatological punches that had a diameter of 5 mm and 10 mm. The samples were cleaned mechanically (using soft toothbrushes and distilled water at high pressure) and chemically (using sodium hypochlorite [NaOCl, 5–6.5%]), according to the procedure described by Reif [26]. The samples were then visualized under a Leica Z16/DFC500 macroscope and analysed using ImageJ software [49] at the Cavanilles Institute for Biodiversity and Evolutionary Biology at the University of Valencia (Spain). For each sample, we measured the scale crown insertion angle (CA) (the angle of the crown of the scale to the vertical axis), which was averaged from measurements taken of ten scales that were randomly selected from each sample, and the scale coverage (SC) (the percentage of the skin surface covered by scales), which was measured in an area of 2 × 2 mm located at the centre of each sample (figure 1). These two variables were considered as they were suggested by several authors as potential traits that affect the attachment and survival of ectoparasites [21,44] (see §4 for further details).

Figure 1

Topological positions of the placoid scale samples (upper), and the measurements considered in this study (lower). Topological positions: 1, dorsolateral region posterior to the pectoral girdle; 2, dorsolateral region posterior to the pelvic girdle; 3, ventral region posterior to the pectoral girdle; 4, ventral region posterior to the pelvic girdle; 5, tip of the rostrum; 6, most anterior part of the lower jaw; 7, gill slit between the first and second gills; 8, leading edge of the pectoral fin. Samples 1, 2, 7 and 8 were taken from the left side of the body. Measurements: CA, scale crown insertion angle; SC, scale coverage. The scale in the upper panel is presented in lateral view. The shark drawings are courtesy of Hugo Salais (HS Scientific Illustration). (Online version in colour.)

Phylogenetic comparative methods

The relationship between grouping behaviour, parasitic pressure and squamation traits was assessed in two successive steps, by implementing phylogenetic generalized least squares (pGLS) models in the package ‘caper’ v. 1.0.1 [50] using R [51]. This methodology allows for the fitting of a linear regression to explore the impact of one or several predictor variables on a single response variable while controlling for a potential phylogenetic signal in the response [52]. First, we tested whether parasitic pressure was increased by the frequency of grouping behaviour by performing pGLS regressions between the PR and all possible combinations of the predictors GB1/GB2 and NS. PR was considered to be an approximation of parasitic pressure according to previous studies ([53] and references therein). The NS was incorporated to account for a potential bias of PR due to uneven sampling effort [54]. Due to the possibility that some of the records where PR is equal to zero represent true absence of parasites (i.e. structural zeroes), we checked for zero-inflation in the models using the R package ‘performace‘ v. 0.8.0 [55] and ran zero-inflated count data regressions using the R package ‘pscl’ v. 1.5.5 [56]. Second, we assessed whether the CA and SC varied according to the location on the body, lifestyle and/or frequency of grouping behaviour by performing a second set of pGLS regressions between CA/SC and all possible combinations of the predictors P, EG and GB1/GB2.

In all cases, we included null models (i.e. without covariates) for comparison. The Akaike information criterion (AIC) and likelihood-ratio tests for non-nested models (LRT) were employed to compare the goodness of fit for the different models. LRT were performed using the R package ‘nonnest2’ v. 0.5.5 [57] based on Vuong's [58] theory. Model parameters were standardized using the R package ‘effectsize’ v. 0.6.0.1 [59]. Multicollinearity was checked by calculating variance-inflation factors (VIF) and generalized variance-inflation factors (GVIF) with the R package ‘car’ v. 3.0.12 [60] and considering a threshold of GVIF(1/(2*Df)) = 2 (equivalent to VIF = 4) [61]. All pGLS analyses were conducted with consideration of the phylogenetic tree provided by Vélez-Zuazo & Agnarsson [62] after trimming the taxa that were not included in the respective datasets using the R package ‘ape’ v. 5.5 [63]. Species for which the squamation of more than one specimen was studied were included as polytomies that were resolved randomly with zero branch length. Additionally, all the analyses were repeated using non-phylogenetic generalized least square models (GLS).

Results

Parasitic richness varies considerably among species with different frequencies of grouping behaviour and numbers of published studies (electronic supplementary material, figure S3). The best-supported pGLS model (PR ∼ NS + GB1) indicates that parasitic richness is mostly explained by these two predictors (electronic supplementary material, table S3). This model supports that the taxa that show infrequent grouping behaviour have lower PR than those that show frequent grouping behaviour and that PR increases with the number of published works (table 1). This is also supported by the models that include a finer categorization of grouping behaviour (i.e. PR ∼ NS + GB2), in which the PR increases successively in solitary species, species that show occasional grouping behaviour, and species that show frequent grouping behaviour (electronic supplementary material, data S1). Grouping behaviour affects PR even in models that do not account for the number of published studies, but these explain a lower amount of variation and have high AIC scores (electronic supplementary material, table S3 and data S1). Two additional candidate models (i.e. PR ∼ NS + GB2; PR ∼ NS × GB1) show a similar support on the basis of the AIC scores (i.e. ΔAIC < 2) and LRT results (p-value = 0.73 and 0.64, respectively) (electronic supplementary material, table S3).

Table 1

Details of the best-supported phylogenetic (pGLS) and standard (GLS) regression models, showing the associated coefficients and p-values. CA, scale crown insertion angle; PR, parasitic richness; SC, scale coverage. Predictors: EG, ecological group (School, schooling species of low to moderate speed; Biolu, mesopelagic bioluminescent species; DemRo, demersal species on rocky substrates and in caves; DemSa, demersal species on sandy and muddy substrates; Pelag, strong-swimming pelagic species; Open, slow species of the open water); GB1, grouping behaviour (Infrequent, Frequent); GB2, grouping behaviour (Solitary, Occasional, Frequent); NS, number of studies; P, position (the numbers refer to the topological positions in figure 2).

	pGLS						LM
	PR ∼ NS + GB1		CA ∼ P + EG + GB1		SC ∼ P + EG + GB2		PR ∼ NS + GB1		CA ∼ P + EG + GB1		SC ∼ P + EG + GB2
	estimate	p-value	estimate	p-value	estimate	p-value	estimate	p-value	estimate	p-value	estimate	p-value
intercept	0.54	6.51 × 10−1	134.65	2.20 × 10−16	41.36	5.35 × 10−13	1.81	4.39 × 10−2	129.71	2.00 × 10−16	50.01	2.00 × 10−16
NS	0.02	2.00 × 10−16	—	—	—	—	0.02	2.00 × 10−16	—	—	—	—
GB1	—	—
frequent	0.00	—	—	—	0.00	—	0.00	—	—	—	0.00	—
infrequent	−2.81	1.23 × 10−2	—	—	6.34	1.36 × 10−1	−3.23	8.55 × 10−4	—	—	4.65	6.21 × 10−2
GB2
frequent	—	—	0.00	—	—	—	—	—	0.00	—	—	—
ocassional	—	—	−2.80	5.80 × 10−1	—	—	—	—	−2.59	3.42 × 10−1	—	—
solitary	—	—	−15.46	2.80 × 10−2	—	—	—	—	−12.65	1.07 × 10−6	—	—
P
P1	—	—	0.00	—	0.00	—	—	—	0.00	—	0.00	—
P2	—	—	0.07	9.77 × 10−1	1.26	7.56 × 10−1	—	—	−0.08	9.80 × 10−1	0.90	8.06 × 10−1
P3	—	—	−0.98	6.99 × 10−1	1.84	6.49 × 10−1	—	—	−1.87	5.52 × 10−1	1.81	6.22 × 10−1
P4	—	—	−0.13	9.58 × 10−1	1.36	7.36 × 10−1	—	—	−1.00	7.51 × 10−1	2.15	5.58 × 10−1
P5	—	—	−22.92	2.20 × 10−16	31.46	1.78 × 10−13	—	—	−17.97	2.62 × 10−8	22.39	3.02 × 10−9
P6	—	—	−11.25	1.31 × 10−5	32.20	5.46 × 10−14	—	—	−11.10	4.78 × 10−3	22.12	4.58 × 10−9
P7	—	—	−0.97	7.02 × 10−1	26.49	3.05 × 10−10	—	—	−1.21	7.02 × 10−1	18.97	4.06 × 10−7
P8	—	—	−19.75	1.59 × 10−13	29.34	4.73 × 10−12	—	—	−16.13	5.21 × 10−7	21.33	1.48 × 10−8
EG
School	—	—	0.00	—	0.00	—	—	—	0.00	—	0.00	—
Biolu	—	—	13.94	1.32 × 10−1	−27.25	1.55 × 10−5	—	—	4.76	1.65 × 10−1	−28.52	6.25 × 10−13
DemRo	—	—	−17.42	1.32 × 10−1	36.94	3.70 × 10−6	—	—	−18.07	5.17 × 10−11	32.70	2.00 × 10−16
DemSa	—	—	−26.70	6.72 × 10−3	39.44	2.29 × 10−5	—	—	−18.99	1.76 × 10−8	31.64	2.31 × 10−15
Pelag	—	—	−20.86	1.96 × 10−2	35.39	5.78 × 10−7	—	—	−21.62	5.22 × 10−15	33.25	2.00 × 10−16
Open	—	—	−19.60	2.10 × 10−1	35.11	4.54 × 10−4	—	—	−15.91	1.38 × 10−5	33.02	2.91 × 10−14

Figure 2

Scatterplots representing scale crown insertion angle and scale coverage in various topological positions of the studied shark specimens. The numbers refer to the topological positions in figure 2 and the specimens are coloured by ecological group. (Online version in colour.)

Zero-inflated regressions also support that parasitic richness increases with the frequency of grouping behaviour and the number of published studies, according to the best-supported model (PR ∼ NS + GB2) (count model in electronic supplementary material, tables S4 and S5). No other candidate model provides a similar fit to the data on the basis of the AIC scores and LRT results. Zero-inflation model supports that the odds that a shark species has a parasitic richness of zero significantly increases in solitary species and decreases with the number of published studies (by 5.20 and 0.99, respectively; zero-inflation model in electronic supplementary material, table S5 and data S1).

The scale crown insertion angle and the scale coverage vary significantly among the different species and body locations (figure 2). The scale crown insertion angle is largely explained by grouping behaviour, the topological position on the body and the ecological group, according to the best-supported pGLS model (CA ∼ P + EG + GB2) (electronic supplementary material, table S3). This model supports that the CA increases successively in solitary species, species that show occasional grouping behaviour and species that show frequent grouping behaviour. A similar CA is observed among the scales on the gill slits and dorsolateral and ventral regions of the body, while smaller CA is observed in the scales of the mouth, rostrum and leading edge of the pectoral fins (table 1 and figure 2). Schooling species of low to moderate speed display larger CA than strong-swimming pelagic species, demersal species on rocky substrates and in caves, demersal species on sandy and muddy substrates, and slow species of the open water. However, they display smaller CA than mesopelagic bioluminescent species (table 1 and figure 2). Four additional candidate models (i.e. CA ∼ P + EG; CA ∼ P; CA ∼ P + GB2; CA ∼ P + EG + GB1) show similar AIC scores but LRTs support that the model with the lowest AIC score is significantly different from all of them (p-value = 0.09, 0.03, 0.04 and 0.06, respectively; electronic supplementary material, table S3).

Similarly, the scale coverage is largely explained by grouping behaviour, the topological position on the body and the ecological group, according to the best-supported pGLS model (SC ∼ P + EG + GB1) (electronic supplementary material, table S3). This model supports that species that show occasional grouping behaviour have greater SC compared with species that show frequent grouping behaviour. Similar SC is observed on the dorsolateral and ventral regions of the body, while a greater SC is observed on the gill slits, mouth, rostrum and leading edge of the pectoral fins (table 1 and figure 2). Schooling species of low to moderate speed display less SC than strong-swimming pelagic species, demersal species on rocky substrates and in caves, demersal species on sandy and muddy substrates, and slow species of the open water. However, they show a greater SC compared to mesopelagic bioluminescent species (table 1 and figure 2). Two additional candidate models (i.e. SC ∼ P + EG; SC ∼ P + EG + GB2) show a similar support on the basis of the AIC scores but the p-value from the LRT is comparatively low for the model excluding grouping behaviour (p-value = 0.11; electronic supplementary material, table S3).

In terms of model support, the predictors included in the best-supported models and the significance of the factor levels, the results derived from the GLS models are virtually the same as those obtained from the pGLS models (tables 1 and S3). Candidate GLS models excluding grouping behaviour clearly fit the data less well than those including this covariate, especially when considering the results from LRTs (electronic supplementary material, table S3). We detect no significant multicollinearity among the explanatory variables (electronic supplementary material, data S1).

Discussion

Grouping behaviour involves increased parasitic pressure in sharks

The results derived from the pGLS, GLS and zero-inflated regression models demonstrate that grouping behaviour affects the richness of monoxenic ectoparasites, even when various categorizations of the former are considered (electronic supplementary material, tables S1 and S4). Species that show frequent grouping behaviour have greater parasitic richness than those that display a solitary lifestyle or aggregate occasionally (table 1 and electronic supplementary material, table S5). This agrees with the fact that solitary species might be contributing to zero-inflation (electronic supplementary material, tables S4 and S5). Parasitic richness is a reliable estimate of parasitic pressure and the impact on hosts, as shown by independent evolutionary, genetic, immunological, physiological and behavioural evidence ([53] and references therein). According to this, our study supports that shark taxa that tend to aggregate more frequently are exposed to a stronger parasitic pressure, which might constitute a selective agent for the evolution of defensive mechanisms.

Besides parasitic richness, other metrics such as intensity (i.e. average number of individual parasites of a given type present in an infected host) or prevalence (i.e. percentage of the host population infected with a parasite) are commonly used in studies that aim to quantify ecological and evolutionary consequences of parasitism [64]. The integration of those parameters within the framework here presented could provide novel insights about their interrelationships with parasitic richness and pressure [65]. However, the design of such large-scale comparative studies is hindered by the fact that intensity and prevalence records in sharks are still scarce. Other variables of potential interest, which are equally scant in the literature, include the number and density of individuals in the aggregations or information on group composition (e.g. number of species, age or sex) [1].

The results of this study also show that parasitic richness is correlated with the number of studies published on each shark species, with the PR being greater for the taxa that were studied extensively (those that included a broader bibliography). This indicates a requirement for the incorporation of variables that account for potential bias due to uneven sampling effort when studying and modelling parasitic richness rather than solely considering raw counts of parasite species [54].

Parasite control mechanisms of shark squamation

The pGLS and GLS models reveal that the scale crown insertion angle and the scale coverage are affected by the frequency of grouping behaviour (electronic supplementary material, table S3). Thus, the taxa that frequently aggregate have squamation with a larger CA and less SC than those that are mostly solitary or aggregate occasionally (table 1). We propose that this phenomenon may be a defensive response for the prevention of parasite attachment and to compromise the survival of parasites on the skin of species with increased parasitic pressure associated with grouping behaviour.

In support of this hypothesis, previous reports on ectoparasite attachment locations indicated that there is a preference for certain areas that are characterized by the absence of denticles and/or extensive wear [66]. This seems particularly true for species that have scales with crowns that point upwards, for which parasites are mostly restricted to the mouth cavity, gills and spiracle (e.g. [67-71]). Conversely, research has shown that a dense coverage of scales eases the attachment and migration of the larvae of the main groups of shark ectoparasites [44]. This is also evident in the general occurrence of copepods on the body surface of pelagic sharks that have imbricated tilted denticles (e.g. [67,72,73]). This type of squamation may further provide a suitable microhabitat for the larvae of certain parasitic groups (e.g. Monogenea and Copepoda) that often occupy the free spaces between the lower surface of the crown and the skin [74,75], while a sparse coverage of scales may increase exposure to a more hostile external environment and predators [76]. Therefore, a combination of a large CA and low SC could reduce the effective surface area for the adhesion and migration of ectoparasites and provide an inadequate microhabitat for the larvae.

The squamation of sharks is shaped by a complex trade-off of functional drivers

Besides the frequency of grouping behaviour, the pGLS and GLS models show that the scale crown insertion angle and the scale coverage also vary between the various regions of the body and the lifestyle of each species (electronic supplementary material, table S3). The topological and ecological variation of squamation has been illustrated in several previous studies ([26,27] and references therein). Here, we demonstrate that the rostrum, mouth and leading edge of the pectoral fins have scales with a smaller CA and greater SC than the dorsolateral and ventral regions of the body. Schooling species of low to moderate speed show a larger CA and less SC than the strong-swimming pelagic species, both groups of demersal species, and slow species of the open water. However, they showed a smaller CA and greater SC than the mesopelagic bioluminescent species (table 1 and figure 2).

We propose that this complex variation of patterns regarding shark squamation may have arisen due to an adaptive trade-off [32,77]. When subjected to various selective pressures (i.e. high abrasion, increased drag, increased ectoparasitism, or the presence of photophores), the optimization of the CA and SC seems to follow divergent trajectories. Thus, while a large CA and less SC may provide defence against ectoparasites (see above), a small CA and greater SC may enable abrasion resistance and drag reduction [32]. Our results suggest that morphofunctional adaptations to abrasion resistance and drag reduction prevail over those associated with defence against ectoparasites when increased abrasion or drag occurs. Accordingly, the areas of the body that are subjected to high abrasive stress (i.e. the mouth, rostrum and leading edges of the pectoral fins in all sharks, and most of the body of demersal species) and increased flow velocity (i.e. most of the body of strong-swimming pelagic species) have scales with comparatively small CA and greater SC, despite the fact that many of these pelagic and demersal species often occur in aggregations (figure 2 and electronic supplementary material, table S1) [46-48]. This view is further reinforced by the differences in the standardized coefficient magnitudes of the models, which indicate that these particular lifestyles and body regions exert a stronger effect on CA and SC than the fact of belonging to a species with a high frequency of grouping behaviour (electronic supplementary material, data S1). The scales of these body regions and ecological groups present conserved morphologies, generally as a result of convergent evolution, and are associated with the abrasion-resistant and drag-reduction functional types described by Ferrón & Botella [27].

Conversely, the squamation of bioluminescent sharks is characterized by the sparse presence of scales, enabling the accommodation of photophores without intercepting the light emitted by them [26,33,78,79] (but see also [80]). Our results support this and show that the bioluminescent species have the lowest SC among all the sharks studied (table 1 and figure 2). By contrast to species in other ecological groups, bioluminescent species show significant diversity regarding scale crown morphology (electronic supplementary material, figure S1) [26,27,33]. This fact is reflected in our results in the wide range of CA values observed among sharks of this ecological group (figure 2), suggesting that this trait plays a minor role in bioluminescence aiding. Considering this, no trade-off is expected between the morphofunctional adaptations of squamation associated with bioluminescence aiding and defence against ectoparasites when photophores and high ectoparasitic pressure co-occur in the same species. Our results show that the largest and smallest CA correspond to scales of etmopterid and dalatiid taxa (electronic supplementary material, table S2), which are mostly regarded as schooling and solitary sharks, respectively [26,33]. This indicates that part of the morphological diversity of bioluminescent shark scales could be associated with the frequency of grouping behaviour and the associated increased parasitic pressure. However, this should be cautiously considered as the ecology and ethology of these deep-water sharks remain poorly understood.

Concluding remarks and further implications

Here, we constrain the relationship between grouping behaviour, ectoparasitic pressure and associated defensive mechanisms regarding squamation in sharks through phylogenetic, standard and zero-inflated regression models. We show that a higher frequency of grouping behaviour is associated with increased parasitic pressure, which may trigger the evolution of squamation with a large CA and low SC. These are interpreted as traits that prevent the settlement and compromise the survival of ectoparasites. This phenomenon is less evident in the body regions and ecological groups that experience high abrasive stress or increased drag, suggesting that the squamation of sharks responds to a complex trade-off of functions where protective and hydrodynamic roles prevail over the rest (e.g. ectoparasitic defence and bioluminescence aiding). The analysis of these data in the context of Pareto optimality theory might provide a quantitative framework for explicitly testing this hypothesis [81,82].

The results of our study shed light on the cost–benefit balance that determines the evolution of grouping behaviour in vertebrates. Moreover, we have detected and quantified several squamation traits associated with high ectoparasitic pressure. These results constitute a basis to indirectly infer this phenomenon, not only in living taxa but also in extinct groups and may also enable us to trace the evolution of parasitic and social interactions in the evolutionary history of vertebrates and chondrichthyans in the future. The comparatively good representation of vertebrate squamations in the fossil record ensures a wide applicability of this approach to extinct taxa. Data extracted from this novel methodological framework might be crucial to complement the records of host–parasite body fossils in association [83-89] and pathologies [90,91]. Still, this more anecdotal evidence has the unique potential to provide estimations of parasitic intensity or prevalence as well as information on the nature of ectoparasites and their attachment site preferences [92-97]. The scarcity of such records could be somewhat alleviated not only by dedicating more effort to the recovery of associated parasites but also by implementing less aggressive preparation techniques on potential host specimens [96,97]. Ultimately, the results of this study may serve as a source of inspiration for the development and fine-tuning of biomimetic surfaces that involve simple denticulated designs for antifouling purposes [98-102].