A molecular and ecological study of Grillotia (Cestoda: Trypanorhyncha) larval infection in small to mid‐sized benthonic sharks in the Gulf of Naples, Mediterranean Sea

Abstract Aim Trypanorhyncha cestodes comprise a wide range of heteroxenous parasites infecting elasmobranchs as definitive hosts. Limited data exist on the larval infection of these cestodes and the role of intermediate and paratenic hosts in the life cycle of these parasites. We investigated the factors that determine the occurrence and the level of infection of Grillotia plerocerci in the skeletal muscles of various benthonic sharks and analyzed the parasites through an integrative taxonomic approach. Location Mediterranean Sea. Methods Sharks obtained as bycatch of commercial trawling activities (i.e., Etmopterus spinax, Galeus melastomus, and Scyliorhinus canicula) were used in this study. Data from a limited number of Dalatias licha and Scyliorhinus stellaris were also included. Grillotia plerocerci were molecularly characterized using the partial 28S large subunit rDNA. Boosted regression trees were used to model the relationship between the abundance of infection with both morphological and physiological predictors in each host. Results Plerocerci of Grillotia were detected in all shark species except S. stellaris. Host species significantly differed in terms of parasite abundance, with the highest and lowest prevalence and abundance of infection detected in G. melastomus and E. spinax, respectively. The relative influence of the traits involved in explaining the parasite abundance was related to the host size in G. melastomus, while both morphology‐ and physiology‐related traits explained the patterns observed in E. spinax and S. canicula. The 28S rDNA sequences shared an identity of ∼99.40% with a Grillotia species previously found in the Mediterranean Sea. At intraspecific level, two different genotypes were found. A first type was retrieved only from D. licha, whereas a second type was found in G. melastomus, E. spinax, and S. canicula. Main conclusions Present results suggest that the two genotypes could be involved in different consumer‐resource systems and confirm most of the examined shark species as transport hosts of Grillotia species for unknown larger top predators.

Moreover, no molecular data exist from the Mediterranean Sea except for a recent report of a massive infection of plerocerci by an unidentified Grillotia species in the anglerfish Lophius piscatorius Linnaeus, 1758 (Santoro et al., 2018).
Understanding the factors that determine the occurrence or level of infection in marine fishes has interested scientists for decades and numerous variables about parasite, host, and environment have been identified as important in the dynamics of the host-parasite systems (Dallarés, 2016;Palm, 2004;Santoro, Iaccarino et al., 2020;Timi & Poulin, 2020;Tompkins et al., 2011). In particular, the probability to become infected with larval forms of trypanorhynchans, as well as their levels of infection, depends on a number of biotic and abiotic factors at different geographical scales (Palm, 2004;Palm et al., 2017;. The distribution of larval and adult parasites is generally shaped through biotic factors involved in transmission pathways, such as trophic interactions between definitive, intermediate, and transport hosts. For parasites with a life cycle embedded within the marine food webs, biotic factors are also important to evaluate the physiological condition of the fish and the way the fish get resources from the environment, and also the probability for a host to acquire parasites. Therefore, the study of relationships between biotic factors and number of parasites within a host might also provide information about the biology and ecology of both hosts and parasites (Dallares, 2016;Dallares et al., 2016Dallares et al., , 2017Dallares, Padros et al., 2017;Palm, 2004;Santoro, Iaccarino et al., 2020;Timi & Poulin, 2020).
The aim of this study was to characterize through an integrative taxonomic approach the Grillotia plerocerci that infect small to mid-sized benthonic sharks in the Gulf of Naples (Mediterranean Sea) and investigate the features of the infection evaluating how the biological traits of the hosts can influence the occurrence and abundance of the parasites. Sharks used in the present study constituted the bycatch of scientific and commercial trawling operations (red shrimps' and pink shrimps' fishery activities) held with a commercial fishing vessel equipped with bottom trawl nets (mouth of 3 × 4 m in height and width, respectively; 40 mm mesh size), towed at ~2-2.5 kn on muddy bottoms . Samplings were performed in accordance with the permit n. 0008453 (issued May 15, 2020) by the Italian Ministry of Agricultural, Food and Forestry Policies, the guide for the care and use of animals by the Italian Ministry of Health and the ARRIVE guidelines.

| Collection and shark examination
After thawing, the sharks were weighed to the nearest 0.1 g and measured (fork length-FL) to nearest 0.1 cm; sex was determined by gonadal examination. A macroscopic gonadal maturity score (GMS) was recorded to investigate the phase of the reproductive cycle (1 = immature, 2 = maturing, 3 = mature, 4 = resting/regressing) as described in Follesa and Carbonara (2019 were calculated as suggested by Mouine et al. (2007).
The musculature of each fish was first cut in thin slices (about 0.5 cm in thickness) and examined under a dissecting microscope for trypanorhynchan larvae. Then, each muscle slide was shredded in a Petri dish with purified sea water and all parasites were collected.
Parasites attached to the musculature were extracted using scissors and tweezers. Subsequently, all larvae were counted, washed in physiological saline solution, and preserved in 70% ethanol or frozen (−20°C) for morphological and genetic analyses, respectively.
Larval forms were clarified in Amman's lactophenol, and the main morphological characters of scolex and bothria were used to identify the larvae at the lowest possible taxonomic level, according to the identification keys proposed by Palm (2004).
For SEM analyses, five plerocerci were also fixed overnight in 2.5% glutaraldehyde, then transferred to 40% ethanol (10 min), rinsed in 0.1 M cacodylate buffer, postfixed in 1% OsO 4 for 2 hr, and dehydrated in ethanol series, critical point dried, and sputter-coated with platinum. Observations were made using a JEOL JSM 6700F scanning electron microscope operating at 5.0 kV (JEOL, Basiglio, Italy).
Selected samples of skeletal muscles of infected fish were fixed in 10% neutral phosphate-buffered formalin for histological examination and processed by routine methods into paraffin blocks which were cut into 3μm-thick sections and stained with hematoxylin and eosin.

| Molecular analyses for Grillotia identification
Due to the high number of larvae retrieved, a representative subsample of 80 Grillotia plerocerci, 20 from each positive fish species (S. canicula, G. melastomus, E. spinax, and D. licha), was used for molecular analyses. Total genomic DNA (gDNA) from each larva was extracted using the Quick-gDNA Miniprep Kit (ZYMO RESEARCH), following the standard manufacturer-recommended protocol.
Sequence identity was checked using the Nucleotide Basic Local Alignment Search Tool (BLASTn) (Morgulis et al., 2008). The obtained 28S sequences were aligned with the published sequences of Lacistorhynchidae using ClustalX v. 2.1 (Larkin et al., 2007). JModelTest v. 2.1.10 was used to select the best fit model using the Akaike information criterion (Posada, 2008;Posada & Buckley, 2004). Representative sequences obtained in the present study were deposited in GenBank under the accession numbers MW838227-MW838238.
Bayesian inference (BI) was constructed using MrBayes, v. 3.2.7 (Huelsenbeck & Ronquist, 2001). The Bayesian posterior probability analysis was performed using the MCMC algorithm, with four chains, 0.2 as the temperature of heated chains, 5,000,000 generations, with a subsampling frequency of 500 and a burn-in fraction of 0.25. Posterior probabilities were estimated and used to assess support for each node. Values with a 0.90 posterior probability were considered well supported. Maximum likelihood (ML) analysis was performed using IQ-TREE (Nguyen et al., 2015) with 1,000 ultrafast bootstrap replicates. Clades were considered to have high nodal support if the ML bootstrap resampling ≥70%. Due to the phylogenetic position of the family Lacistorhynchidae (see Palm et al., 2009), the phylogenetic tree was rooted using Nybelinia surmenicola Okada in Dollfus, 1929 as outgroup. The 28S sequences from GenBank included in the phylogenetic tree are listed in Table 1.

| Statistical analyses
Prevalence and abundance of infection used in the present study follow Bush et al. (1997). Only host species with a sufficient number of individuals including E. spinax, G. melastomus, and S. canicula were used for statistical analyses. Differences in the abundance of Grillotia in the three shark species were measured by means of the nonparametric rank-based Mann-Whitney U test.
Boosted regression trees (BRTs) were then used to model the relationship between the number of Grillotia specimens found in each host with both morphological (FL, weight, and BCI) and physiological predictors (age, sex, SMS, GSI, and HSI) (Santoro, Iaccarino et al., 2020). BRTs are a class of nonparametric machine learning method where regression trees are combined in an ensemble boosting algorithm to improve predictive performance. Cross-validation (with max trees = 10,000) was used to find the optimal settings in terms of model parametrization and to avoid overfitting due to the relatively low number of specimens in our samples (De'Ath, 2007;Elith et al., 2008;Leathwick et al., 2006). Models were therefore trained by using values of learning rate (lr, the contribution of each regression tree to the model) of 0.01 and tree complexity (tr, model complexity in terms of allowed interactions among predictors) of 5, with bag fraction values (i.e., the proportion of observations used in selecting variables) of 0.5 and assuming a Gaussian error distribution.
Contributions of predictors were expressed by their relative influence (Friedman & Meulman, 2003), scaled in a 0-100% scale where higher percentage indicates a stronger influence on the response variable and can thus be considered more influential. Partial dependence plots were used to model the effect of predictors on the response variable after accounting for the average effects of all other variables in the model (Friedman, 2001;Friedman & Meulman, 2003).
Model performance was measured by using the explained deviance, computed as 1-(mean residual deviance/mean total deviance) and with values of 1 indicating a perfect fit. Model accuracy was assessed by means of the Spearman's rank correlation values between fitted and observed values. Finally, the H-statistic was used to measure the strength of interaction effects between predictors in BRT models, ranging between 0 and 1 and with higher values corresponding to larger interaction effects. Analyses were performed in R (R Development Core Team, 2018) using the packages gbm and dismo (Hijmans et al., 2013;Ridgeway et al., 2017) and following the prescriptions in Elith et al. (2008) and Elith and Leathwick (2017).

| General data
Numbers of examined hosts, biometrical data, physiological indices, and basic values of infection for all host species examined are listed in Table 2

| Molecular identification of Grillotia plerocerci
The The BI and ML tree topologies as inferred from the phylogenetic analysis showed that the 28S sequences clustered in a single clade, which also includes the sequence MH664054 deposited by Santoro et al. (2018). However, the two genotypes obtained and the sequence of Grillotia species already available in GenBank clustered in two separate branches, too ( Figure 5). The distance value between the sequences obtained here and the MH664054 sequence was 0.19%. The sequences included in this clade showed a close relationship with the sequence of Bathygrillotia rowei (Campbell, 1977) Beveridge & Campbell, 2012 (syn. G. rowei) retrievable in GenBank (DQ642765). Indeed, the distance values between the present Grillotia specimens and B. rowei resulted to be 3% and 3.1%, thus not guarantying conspecificity.

| D ISCUSS I ON
According to the GenBank database, 28S rDNA molecular data exist only from three (G. yuniariae, G. erinaceus, and G. pristiophori) out of the 18 species of Grillotia known at present, with a single 28S rDNA sequence available for an unidentified Grillotia species from the Mediterranean basin (Santoro et al., 2018). Molecular comparisons of the present plerocerci with data retrieved from GenBank confirmed that the present material did not match with the three above mentioned species. Moreover, sequence comparisons between plerocerci from D. licha and the other three shark species (E. spinax, G. melastomus, and S. canicula) documented differences which provided support for the existence of two potential genotypes, each of which could have a different definitive host and could be involved in different food web systems. In fact, whether D. licha preyed on small F I G U R E 6 Results of the boosted regression tree models showing the relative influence (in percentage) of morphological (FL, weight, and BCI) and physiological predictors (age, sex, SMS, GSI, and HSI) on the number of Grillotia specimens found in each host species. BCI, body condition index; HIS, hepatosomatic index; FL, fork length; GSI, gonadosomatic index; SMS, sex maturity score. Sharks in the figure have been redrawn from Compagno (1984aCompagno ( , 1984b and downloaded from https://fishb ase.org under a Creative Commons License-CC BY-NC 3.0. Dotted boxes on shark tails show approximately the body sites where Grillotia specimens have been collected sharks (e.g., G. melastomus and S. canicula) as suggested by previous studies in Spanish waters (Matallanas, 1982;Navarro et al., 2014), it would show at least a mixed infection with both genotypes of Grillotia sp. found in the present study. In support of our hypothesis, it has been suggested that parasites, in the process of host switch adaptation, may produce different genotypes that involve substitutions on a small number of nucleotides. Such minor genetic modifications can then become associated with different host species (Poulin et al., 2011;Woolhouse et al., 2001).
Life cycle proposed for members of Grillotia includes three/ four hosts, with copepods as first, schooling fishes as second, and larger predatory fishes as third intermediate hosts, with species of elasmobranchs acting as their definitive hosts (Beveridge & Campbell, 2007;Lubieniecki, 1976;Palm, 2004 (Lubieniecki, 1976;Palm, 2004). An association be- which in turn are known to prey on copepods (Dallarés, 2016). In the present study, we did not focus the attention on the quantification and identification to the lowest taxonomic level of the gastrointestinal contents. However, gross analysis of gastrointestinal material during the parasitological analysis revealed that remains of myctophid fishes and cephalopods were among the most common food items found in E. spinax, G. melastomus, and S. canicula, supporting a role of these prey as possible transport hosts for Grillotia specimens.
In Mediterranean Sea, recent records of larval forms of Grillotia infecting sharks include two reports in G. melastomus, in which the larvae were originally reported as Grillotia sp. (Dallares, Padros et al., 2017) and then as G. adenoplusia (Dallares et al., 2017), with a prevalence of infection ranging from 7.6% to 28.6% and a maximum mean abundance of 0.48 depending from the sampling localities, whereas no Grillotia specimens were detected in 41 individuals of S. canicula from the same locality (Dallares, Padros et al., 2017;Dallares et al., 2017). The difference in infection values between the North-Western Mediterranean studies (Dallares, Padros et al., 2017;Dallares et al., 2017) and the present one might be related to the different geographical features and high productivity in terms of food webs at the time of the sampling, two variables that may affect the number of intermediate and definitive hosts present in the area (Collins et al., 1984;Dallares, Padros et al., 2017;Hassan et al., 2002;Palm & Overstreet, 2000;Santoro et al., 2019;Santoro, Iaccarino et al., 2020).
The habitat use and the feeding ecology of the hosts might explain some results obtained in the present study (Jakob & Palm, 2006;Palm & Caira, 2008;Palm et al., 2007), among which: (i) the absence of Grillotia plerocerci in S. stellaris, which mostly use different habitat and prey than other shark species examined in the present study; (ii) the occurrence of the same Grillotia genotype in the three sympatric shark species with similar feeding habit (E. spinax, G. melastomus, and S. canicula); and finally, (iii) the occurrence of a different genotype in D. licha which has different feeding habit than the previous three species. Dallares, Padros et al. (2017) (Palm, 2004;Palm et al., 2007Palm et al., , 2009Palm et al., , 2017. The hunting strategy followed by the bluntnose sixgill shark is a tail-on approach preying from behind, which would presumably enhance the transmission of parasites concentrated in the tail region (Dallarés, 2016;Santoro et al., 2018;Seamone et al., 2014).
Several authors included either the length or the weight of the fish as an explanatory factor for trypanorhynch distributions with strong positive correlation between fish size and trypanorhynch abundance or intensity (see Palm, 2004). In the present study, the weight of the host was the most important predictor for the overall parasite abundance in G. melastomus and E. spinax and the second most important (after GSI) in S. canicula, while the length of the host alone showed negligible importance in all host species, with BRT values ranging from 4.8% to 12.9% ( Figure 6). Moreover, in G. melastomus, also the BCI (based on length and weight relationships) was a discrete predictor of abundance (18.8%) (Figure 6b), suggesting that in general body size is an important predictor for the overall parasite abundance.
Assuming that fish size increases with age, larger and older fishes are expected to accumulate in muscles a larger number of larval forms of parasites with long life span. MacKenzie (1985MacKenzie ( , 1987  and capable to accumulate in the host musculature over time (Collins et al., 1984). Results described in the present study agree with data from previous studies also regarding the scarce immune response in presence of large amounts of larvae (the maximum number of larvae found was 1421 in an adult male G. melastomus), suggesting that the present lesions are probably the result of direct tissue damage rather than an immune reaction targeted toward the parasitic antigens.
Present results support the hypothesis that larger fishes, having in general a higher food intake, consume more infected preys and have greater prevalence and abundance of food-transmitted parasites, and that the influence of host size on parasite abundance may apply especially to long-living parasite larvae. Moreover, larger hosts might better tolerate high parasite loads in terms of body mass, offering more space for parasite distribution, and this may be especially relevant for larval forms of trypanorhynchans which encyst in skeletal musculature (Collins et al., 1984;Hassan et al., 2002;Palm, 2004;Palm & Overstreet, 2000).
All these considerations, and in particular the smaller body mass, might also explain why E. spinax has lower abundance of infection.
Among the indices for evaluating the physiological condition of the fishes (GSI, HIS, and SMS), only GSI was a discrete predictor in S. canicula, with a value of 29.9%, while HSI (14.2%) showed lower relative importance ( Figure 6c). GSI is a good indicator of reproductive activity of fish, while the HSI is strictly related to the gonadal maturation and in general increases at early ripening and decreases at late ripening (Rizzo & Bazzoli, 2020). In the present study, most of the specimens of S. canicula were in reproductive activity (70.2%) suggesting that this physiological condition can modify the immune responsiveness increasing susceptibility to parasitic infections due to an increase in circulating levels of sex steroid hormones. For example, at time of reproductive activity, salmonid fish demonstrate immune deficiencies that include an inability to produce isohaemagglutinins that thus make fishes more susceptible to the ectoparasitic infestations (Pickering & Christie, 1980). Alternatively, it is also plausible that shark individuals in reproductive activity feed more acquiring more parasite individuals.

| CON CLUS ION
In conclusion, this study describes the first molecular characterization of larval forms of Grillotia in sharks from the Mediterranean Sea, infecting multiple hosts with potential evidence of the development of host specific lineages. According to the host-parasite list of the Shark-References Database (https://shark -refer ences.com/speci es/host-paras ites-list) and formerly published studies, this represents the first report of Grillotia specimens in E. spinax and S. canicula, and the first time that this larval form is found in D. licha. This is also the first contribution to the study of factors which may influence the occurrence and abundance of Grillotia specimens in sharks. Finally, the common occurrence of Grillotia, showing high prevalence and abundance among four species of small to mid-sized deep sea benthonic sharks in the Gulf of Naples, enlarges its range of distribution, confirming its wide presence in the deep bottoms of the Tyrrhenian Sea, and indicates that these shark species act as transport hosts of Grillotia species for larger top predators.

ACK N OWLED G M ENTS
Sampling in the Gulf of Naples was supported by the project ADViSE (PG/2018/0494374). The Maglione family (Fabrizio, Francesco, Giuseppe, Salvatore, and Vincenzo) (Giovanni Padre fishing vessel, Naples, Italy) offered the highest possible support during trawling activities.

CO N FLI C T O F I NTE R E S T
The authors report no conflicts of interest.