Hooked on you: shape of attachment structures in cymothoid isopods reflects parasitic strategy

Background Parasite attachment structures are critical traits that influence effective host exploitation and survival. Morphology of attachment structures can reinforce host specificity and niche specialisation, or even enable host switching. Therefore, it is important to understand the determinants of variation in attachment structures. Cymothoid isopods are striking ectoparasites of fishes that include the infamous ‘tongue-biters.’ They are known to parasitise hosts in one of four qualitatively distinct anatomical regions. Here, we quantify variation in cymothoid attachment structures — hook-like appendages called dactyli — and test whether differences in dactylus shape are correlated with parasite mode (where they attach), allometry, or both, using multivariate ordinary least squares regression. We also assess the influence of shared ancestry on shape using a molecular phylogeny to weight our models using phylogenetic generalised least squares regression. Results We find clear differences in shape between externally-attaching and internally-attaching cymothoids but also between anterior and posterior dactyli across various species with the same attachment mode. Allometric effects are significant for anterior but not posterior dactyli. Mouth-attaching species show greater shape variability than gill- and mouth-attaching species. We find no evidence that there are clade-specific patterns of association between parasite mode and dactylus shape. Conclusions Parasite mode appears to be the main driver of attachment morphology. This likely reflects several components of parasite ecology including feeding and functional demands of attachment in different microhabitats. Geometric morphometric approaches to the quantification of shape variation of simple structures is an effective tool that provides new insights into the evolvability of parasite attachment.

Parasitic strategy is largely conserved within cymothoid genera [5,8], but there is substantial variation in microhabitat use between species with the same parasitic mode, because distinct locations or orientations are used by different species. For example, at the genus level, Anilocra and Nerocila comprise exclusively skin-attaching species, but are found on the anterior and posterior regions of their hosts, respectively [9,10]. Within Anilocra, species are often site specific across the anterior region, for example A. haemuli is always found attached near the eye, while A. acanthuri attaches under the mouth [10]. Very few species are known to use more than one attachment mode, with all such examples recorded from atypical host associations. For example, the most common Brazilian freshwater species, Braga patagonica, is a branchial parasite of several fish species but on cultured Colossoma macroponum it is regularly recorded externally attached behind the dorsal fin [11]. The evolution of parasitic mode within Cymothoidae remains unclear, but there is some consensus that each mode has evolved more than once, and that skin-attachment is unlikely to be ancestral for the group as a whole [12,13].
Cymothoids are well adapted for ectoparasitism on their mobile fish hosts; a thickened cuticle affords protection; increased surface area of gill-bearing pleopods facilitates oxygen transfer in gill-and mouth-attaching species, and modified mouthparts enable the acquisition of blood meals in skin-attachers [7,14]. Externally-attaching species, relative to other parasitic modes, exhibit dorsoventral flattening, which reduces drag and minimises the energy expenditure of their hosts [14]. Crypsis is also displayed by some externally-attached species as a strategy to avoid predators such as cleaner fish [15]. The appendages cymothoids use for attachment are particularly characteristic: each of their prehensile pereopods, 'walking' limbs, terminate in a recurved dactylus, presence of which is a synapomorphy of cymothoid species [8]. Cymothoids and their close relatives exhibit a wide spectrum of trophic dependency from free living species, through temporary to obligate parasites, and cymothoids are thought to have evolved from either a cirolanid-like or an aegid-like ancestor [7,8,16]. Species in the families Cirolanidae and Aegidae do not possess recurved dactyli on their posterior pereopods, and adults retain their ability to swim [17]. In contrast, cymothoids lose the ability to swim after they have infested a suitable host, which drastically reduces the probability of finding another host in the event of being displaced [1]. Loss of swimming appendages may have evolved in concert with the origination of cymothoid dactylus morphology as a trade-off between an increasing reliance on host resources and maintenance of traits for the acquisition of new host individuals [14].
As an important trait for facilitating obligate ectoparasitism in cymothoids, we hypothesised that variation in dactylus shape would reflect differences in the functional demands of parasitising hosts in different locations. Gill-and mouth-attaching cymothoid species use their dactyli to penetrate host tissue, but also as hooks to clasp gill-rakers, tongues, or the upper palate [5,[18][19][20]. In contrast, externally-attaching cymothoids use dactyli to anchor themselves to host musculature and dermal tissues and are subject to greater hydrodynamic forces. We predicted that the externally-attaching species would have dactyli that are relatively longer, thinner, and 'needlelike' adapted for piercing flesh, while those of gill-and mouth-attaching species will be stouter, more recurved, and strengthened for 'gripping' . To test these predictions we used a geometric morphometric approach to quantify dactylus shape and assessed the influence of parasitic mode, size allometry, and phylogeny on shape variation.

Specimens
Cymothoid specimens used in this study are from collections at the Water Research Group, North-West University, Potchefstroom, South Africa, and the University of Salford, Manchester, UK. We took images of the first pereopod (P1) from 124 individuals across 18 species, and from 135 individuals of 19 species for the seventh pereopod (P7). Only P1 and P7 were measured, since these are recognised as the most useful for taxonomic studies because of considerable morphological variation between species; P2 − 6 show much less shape variation between species [20][21][22]. All specimens were adult females, each species was represented by at least three individuals, and there was a minimum of 26 individuals for each of three parasitic modes (Table 1). Flesh-burrowing specimens were not included due to insufficient sample numbers and because the seventh pereopods do not possess a recurved dactylus, only a simple stub [23]. All sampled species are known to occur in the ocean around southern Africa except for Anilocra chromis and Anilocra physodes which are found in the Caribbean and Mediterranean, respectively.

Image and landmark acquisitions
We captured high resolution digital images of the P1 and P7 dactyli for each individual using a Nikon DS-Fi1 camera fitted to a Nikon SMZ1500 stereoscopic microscope. For each of P1 and P7, we plotted 39 semi-landmarks to describe two curves between three fixed landmarks ( Fig. 1) with tpsDig2 [24]. The first fixed landmark was located at the medial junction with the propodus, and the second landmark was placed at the distal tip of the dactylus. Due to differences between individuals in how the propodus overlaps and obscures the dactylus, we drew a line between these first two landmarks, and another at a 5 o angle from this. The third fixed landmark was placed on the lateral edge of the dactylus at the intersection of the 5 o line, thus removing joint shape information from the same relative point in each specimen. The first curve was plotted between the first and second landmarks, along the medial edge of the dactylus, re-scaled by length with 13 semi-landmarks. We used the same method for the second curve, between the second and third landmarks, but with 26 semi-landmarks as this edge is between 150-200% the length of the medial curve. Thirty images were plotted a second time to calculate landmarking error.

Phylogeny reconstruction
We generated a new phylogenetic tree using molecular sequence data from the barcode region of mitochondrial Cytochrome Oxidase subunit I (COI: [25]). A sequence was obtained from one representative of each species, for which an image was captured, except Ceratothoa famosa, Elthusa sp., Nerocila sigani, and gen. nov. et sp. nov. There are no publicly available sequences for these species, and the preservation condition of our specimens did not produce DNA sufficient for PCR. For better preserved specimens, we extracted DNA from a single pereopod using a Machery-Nagel spin column kit, before amplification with Polymerase Chain Reaction (PCR) following the protocol in [22]. Each PCR product was purified, then sequenced in both directions on an ABI 3630 Genetic Analyzer, and we generated consensus sequences with Geneious R10 (https://www.geneious.com). We added sequences for Nerocila depressa (MH425627) and Anilocra chromis (KY562736), and used the aegiid, Aega psora (FJ581463), as an outgroup. We aligned nucleotide sequences with TranslatorX [26] and used Gblocks v0.91b [27] to remove ambiguously aligned sites resulting in a trimmed alignment of 581bp.
Due to the length of alignments and our concern that phylogenetic signal might be limited, we applied a backbone constraint to our tree searches based on the latest available molecular phylogeny for cymothoids [12]. This allowed us to better place taxa for which we have no phylogenetic information and to calculate branch lengths. Specifically, these included hard constraints on genera represented by more than two species, while restricting the potential placement of other taxa within those clades. The best fitting substitution model for the alignment (GTR + ) was selected using corrected Akaike's Information Criterion (AICc) in jModelTest 2.1.10 [28,29]. Tree searches were performed under a Maximum Likelihood (ML) framework with RAxML v8 [30], and topological support was assessed with 500 rapid bootstrap replicates. The topology of the constrained tree was then compared against the unconstrained topology using the Shimodaira-Hasegawa test [31]. We ultrametricised the tree using the function 'chronos' of the ape package in R under a correlated trait model, and an optimised value for the smoothing parameter (from 5,000 starting values) was selected as that which produced the tree with the highest penalised log-likelihood [32].

Statistical analyses
All subsequent analyses used R v3.4.4 [33] and the packages ape 5.2 [34], geomorph 3.0.7 [35], phytools 0.6 [36], nlme 3.1 [37] and stats [33]. We obtained individual shape variables from our raw landmark coordinates by Generalised Procrustes Analysis (GPA; [38], using Procrustes distance to optimise locations of semi-landmarks [39]. This superimposition produces a set of scaled and aligned Procrustes coordinates that minimises location, orientation, and size differences between samples, thus retaining only information related to geometric shape. The mean consensus shapes for each species in both the P1 and P7 configurations were aligned in a second GPA to produce species level Procrustes coordinates [40]. We then performed Principal Components Analysis (PCA) on the aligned coordinates to visualise shape differences plotted as morphospaces and PCA backtransformations [40]. The PCA scores were also used as shape variables in regression analyses, for which we retained the first twelve principal components with non-zero eigenvalues. Together these explained over 99% of variation in both our P1 and P7 datasets. For phylogenetic analyses we generated sets of shape variables for P1 and P7 with the same method including only the species present on our phylogeny, and using phylogenetic PCA (pPCA) -P1 phy and P7 phy datasets.
At this point, we quality-checked our full P1 and P7 datasets to identify landmark or analytical problems. First, for each dataset, we calculated the mean Procrustes distance of each sample to the P1 or P7 global consensus shape, where outlying data points might indicate landmark error. Outliers were retained if they were consistent within a species. For example, in the P1 dataset individuals of Norileca indica all appear above the upper quantile, which reflects genuine shape information, rather than error. Error in our data acquisition steps was assessed by nested ANOVA between our 30 repeated landmark sets to calculate the ratio between total Mean Squared Error (MSE) and that contributed by our replicates. We found our digitisation to be over 97% repeatable, and, therefore, that measurement error did not significantly influence variation in shape.
For both our P1 and P7 datasets, we first fitted a 'full' model including all covariates and their interactions, with single factor non-parametric Multivariate Analysis of Covariance (Procrustes npMANCOVA). The twelve principal components describing shape formed our response and as predictors we modelled size as a continuous covariate, parasitic mode as a three-level factor, and the interaction between size and mode. For size, we used the mean of log-transformed centroid sizes for each species, which are the square root of the sum of squared distances between the landmarks of each specimen and their centroid [41]. We then proceeded with pairwise tests to assess homogeneity of slopes (whether patterns of shapeallometry are common across parasitic modes), and group means (testing shape differences between groups after accounting for variation in size). Finally, for both datasets we also used the residuals from a regression of shape on size (equivalent to allometry-free shapes) to test shape differences between parasitic modes. This latter analysis is only appropriate where there is a common allometry between groups [42]. For each analysis we used 10,000 permutations of the Residual Randomisation Permutation Procedure (RRPP [43]) to generate empirical sampling distributions for significance testing, and from which effect sizes were estimated as standard deviates (z-scores).
We applied the same modelling approach to our P1 phy and P7 phy datasets but using Phylogenetic Generalised Least Squares (PGLS) rather than Ordinary Least Squares regression (OLS). Data collated from related species violate the assumption of OLS that residual error is independent between observations [44]. It is possible to account for this autocorrelation by weighting the error structure of the regression assuming a model of trait evolution to calculate covariance of traits among species [45,46]. The default method for PGLS in 'geomorph' (function procD.pgls) assumes traits evolve under Brownian motion, which may not be realistic for all datasets. Therefore, we estimated Pagel's lambda as a measure of phylogenetic signal in the residuals of our 'full' models using Restricted Maximum Likelihood (REML) implemented in the function 'gls' . Using these estimated values of lambda (P1 phy = 0.18 and P7 phy = 0.13 ) we then calculated phylogenetic variance-covariance correction matrices and conducted PGLS with 'procD.lm' . Calculating separate variance-covariance matrices in this way allowed us to complete all regression analyses with 'geomorph' functions and the RRPP method, while using a better fitting model. To infer patterns of morphospace occupation across the evolution of parasitic modes we projected our phylogeny into morphospaces to create phylomorphospaces for P1 phy and P7 phy .

Dactylus shape variation
Principal Component Analysis (PCA) clearly separates dactylus shape of externally-attached cymothoids from gill-and mouth-attaching species for both P1 and P7, where they occupy a distinct region of morphospace (Fig. 2). The vast majority of shape variation is accounted for by the first two principal components: 77.3% and 79.3% for P1 and P1 phy , 88.5% and 89.7% for P7 and P7 phy . For P1, higher PC1 values indicate increased curvature, while PC2 reflects changes in the width of the dactylus, especially the proximal width. Interestingly, all externallyattaching species have positive PC1 values for P1, and gillattaching species share similar curvatures with a narrow PC1 range. P7 dactyli are thinner with margins subparallel at lower values of PC1, and are flatter for negative values of PC2 (Fig. 2). Similar to P1, gill-attaching species occupy a small area of P7 morphospace, while mouth-attaching species cover the widest area, indicating a broad range of dactylus shapes.

Relationship between dactylus shape, size and parasitic mode
Dactylus shape is significantly and strongly correlated with parasitic mode, as evidenced by the positioning of species in morphospace and from the results of Procrustes npMANCOVA ( Table 2). The interaction term between size and mode has a weak and non-significant effect for P1 (z-score = 1.09, p-value = 0.86), whereas both size and mode exhibit large and significant effects (z-score = 1.92 and 2.98, respectively). For P7 we find no significant effect of size on shape (z-score = 0.81, p-value = 0.23), but mode (z-score = 2.82, p-value <0.001) and the interaction of size and mode (z-score = 1.66, p-value = 0.03) are both significant. Despite the significant interaction term in P7 it accounts for little of the total variance in the model (R2 = 0.09). Pairwise comparisons of angles between group allometric slopes show that the size-mode interaction in P7 is driven entirely by differences between external-and mouth-attaching species (z-score = 1.85, p-value = 0.01). As expected from results of the full model, P1 comparisons of group slopes show no pairwise differences. For both appendages there are significant differences between the mean shapes of externally-attaching cymothoids and the other two modes, but not between gill-and mouth-attaching species.

Phylogenetic context
We did not find any significant differences between the log-likelihoods of our constrained and unconstrained topologies using Shimodaira-Hasegawa test(diff = -0.03,   (Fig. 3a). As expected, the constrained topology is consistent with that of [12] except we do not recover a sister relationship between Cymothoa and Nerocila. Evidence from previous work continues to show that evolution of attachment mode in cymothoids is homoplastic, with gill and external attachment likely to have arisen independently at least twice or for there to have been secondary reversals [12,13]. The first two components of pPCA for P1 phy account for 48.4% and 30.8% of variation in shape respectively, while for P7 phy pPC1 reflects 47.7% and pPC2 42.2%. Mouth-attaching genera occupy a much greater portion of phylomorphospace than other modes and exhibit similarly striking patterns as to our full datasets: Cymothoa eremita, for example, exhibits a P7 morphology most similar to external species than to the congener C. sodwana, (Fig. 3c). Gill-attaching species have highly similar P7 morphologies but Norileca indica has a P1 morphology distinct from the other gill-attaching species, perhaps due to its unusual attachment orientation on the ventral side of the operculum, facing anteriorly [47]. Most gill-attaching species are found anchored to the gill filaments with their dorsal side facing the operculum. Accounting for phylogeny with PGLS results in large and significant shape differences between parasitic modes for both P1 phy and P7 phy (Table 3). There is no evidence of separate allometric slopes between groups for either dataset using PGLS but, as with P1, allometry remains an important influence on shape for P1 phy . We find the same pattern in pairwise comparisons as with OLS regression from our full datasets, where the mean shapes of gill-and mouth-attaching species are significantly different from externally-attaching species, but not from each other.

Discussion
Our results indicate that the shape of cymothoid dactyli are strongly influenced by parasitic mode corroborating theory that adaptation to hosts and local environments is a driver of parasite phenotypes [48]. Other empirical studies also link the morphology of parasite attachment to various ecological factors: in relation to host size in feather lice [49]; host specialisation in platyhelminth fish parasites [50]; host thermal regulation strategy in acanthocephalans [51]; host biogeographic plasticity [52]; and, as we also find here, microhabitat [53,54]. However, phylogeny may also explain morphological variation where more closely related species are expected to share similar traits. For example, [55] and [56] both suggest that integration between parts of monogenean attachment organs results in phylogenetic constraint of shape variation. Both of these studies find that divergent parasite species that infest the same host species possess differently shaped attachment organs (called haptors). We have found no evidence that there are clade-specific patterns of association between parasite mode and dactylus shape. Our results imply that different clades can converge on the same dactylus shape, which presumably is well-adapted for that particular parasite mode. Large effect sizes are consistently observed for parasitic mode in both OLS and PGLS analyses of P1 and P7 datasets. Underlying this pattern are significant shape differences between externally-attaching species, and both gill-and mouth-attachers. Interestingly, however, the shape dissimilarities of external and internal species in morphospace are also different between P1 and P7 dactyli. We also find that the shape of P1 dactyli is influenced by size, but this is not the case for P7. Altogether, this suggests that anterior and posterior dactyli may function differently, and in externally-attaching species the differences are particularly acute. Certainly, pereopods are broadly arranged in two opposing angles between pereopods 1-3 and 4-7, which is suggested to enhance  [8,14]. P1 dactylus shapes are highly recurved, whereas for P7 the shape is flatter and more slender than those of gill-and mouth-attaching species. Gill-and mouth-attaching species, on the other hand, possess more similar P1 and P7 morphologies. All isopods exhibit biphasic moulting [57]; a particularly useful trait for cymothoids in that they do not need to completely detach from the host during ecdysis. It has previously been suggested that this preadaptation may have been key to isopods evolving ectoparasitic lifestyles [14]. The risk of detachment during ecdysis is presumably far greater for external-attaching species and posterior dactyli could provide secure anchorage if they function to pierce host tissue, as P7 morphology suggests they might. Inhabiting the external environment seems to result in convergence of P7 dactylus morphology, where relatively less related species share similar morphologies. Externalattaching species have a narrow range of PC1 scores, and Anilocra capensis shares a more similar morphology with A. chromis and Nerocila depressa than its closest relative in our phylogeny, A. physodes. Similarly, gill-attaching species share very similar P7 morphologies. Welicky et al. [58] and Pawluk et al. [59] have shown that cymothoid and host size are tightly correlated except for gill-attaching species. This suggests that available space in the opercular cavity does not scale linearly with host size much in the same way that, to draw parallels with eye anatomy, the upper and lower tarsi are near flush to the sclera or cornea regardless of size. Such a restriction could limit growth patterns of gill-attaching species, including their dactyli. Morphology of gill-attaching species appears to neatly fit their location and orientation within the gill cavity, where they are typically asymmetrical, twisted to one side, and either completely flat or extremely concave. Therefore, there could be two constraints operating on P7 morphology: exposure and detachment risk related to external parasitism, and for gill-attaching species the volume of the opercular cavity.
We observe the opposite pattern in mouth-attaching species, which exhibit highly divergent P7 forms. Free from hydrodynamic forces and the risks of being brushed or cleaned off, a fish mouth is, perhaps counter-intuitively, a safe environment in which to reside as an ectoparasite. Mouth-attaching species have adapted to 'bite tongues' but they also secure themselves to the upper palate, as in Cinusa tetrodontis [18], wrap their dactyli around the tongue, or attach to the inner cheek and gill arches. Most mouth-attaching species are oriented anteriorly, but a few species face posteriorly into the throat, like all species in the genus Isonebula [60]. Therefore, mouth-attachment as a discrete parasitic mode actually conceals considerable microhabitat variation that could explain the diversity of P7 dactylus shapes. In addition, host specificity and the manner of parasite speciation might also be important. Except for Cymothoa eremita, each mouth-attaching species in our dataset is known from fewer than three host species, and even C. eremita shows strong preference toward Parastromateus niger [5,20]. Host switching is a common model explaining parasite speciation, but only a subset of potential new hosts may be a suitable match for parasite phenotypes and it is likely that new host environments are initially suboptimal. After a host-switch event, subsequent co-adaptation between the parasite and the new (suboptimal) host environment would then refine morphological features, potentially driving shape differences between parasite species.
Another possible determinant of cymothoid attachment morphology is feeding ecology. It has previously been suggested that while externally-attaching species depredate hosts, feeding on blood and tissue, gill-and mouthattaching species may additionally take up energy from prey their hosts consume [61][62][63][64]. Certainly, the relatively benign effects of many gill-and mouth-attaching species could support this view [65,66]. Differences in feeding ecology are also consistent with mouthpart morphology, where gill-and mouth-attaching cymothoids possess molar processes and incisors adapted to grind or slice tissues. Externally-attaching species, however, have mouthparts that interlink to form a functional sucking cone, adapted for fluid intake [14]. In externally-attaching species, the incisor process is narrowed to a point rather than a blade, perhaps for piercing host flesh, while the mouthparts are angled downward toward the host. Furthermore, the primary mouthparts are restricted in their movement, only able to motion in an inwardoutward direction. Such a movement requires a suitable counter-force to hold the mouth close to the host [14]; a role that might be performed by the anterior pereopods and, in part, determine the P1 shape of externallyattaching species.
Size is also a significant predictor of P1 shape in OLS and PGLS analyses. As stated previously, the external environment must impose unique challenges, not least of which is hydrodynamic force. We suggest that the directionality of this force could be key to understanding the different influences of size between P1 and P7 dactylus shapes. All externally-attaching species orientate anteriorly, parallel with their hosts [10,16]. In our dataset, they also possess the largest and most recurved P1 dactyli (mean centroid size = 7.54). The anterior dactyli in cymothoids are angled in a manner ideal for pulling against the direction of flow, which requires less effort than pushing due to reduction of frictional force. Such efficiency is likely complemented in externally-attaching species by increased strength derived from P1 dactylus shape and size. Drawing parallels from structural engineering, we can conceive that steeper arches with more symmetrical parabolic geometries are able to resist greater forces in compression [67]. Gill-and mouthattaching species, on the other, hand are not subject to the same demands of the external environment, and possess P1 dactyli that are not as recurved and are perhaps better suited to withstand other structural displacements like shear.

Conclusions
We used geometric morphometrics, multivariate analyses, and phylogenetic comparative methods to quantify shape variation in the attachment structures of cymothoid isopods to determine whether shape differences between species correspond to parasite mode. In addition, we assessed the relative influences of allometry and shared ancestry. We found that ecology is the primary driver of dactylus shape. Separate lineages appear to have independently evolved similar dactylus morphologies that are presumably optimal for particular parasitic modes. The clearest differences are between externally-attaching species, and those attached internally in the gill or mouth. Geometric morphometrics is a powerful method for uncovering complex patterns from simple outline shapes like dactyli. Of particular note are the shape differences we found between anterior and posterior pereopod shapes. These anterior-posterior differences are also characteristic of the different parasitic modes and likely reflect adaptations to obligate parasitism and differences in feeding ecology between externally-attaching, and gill-and mouth-attaching cymothoids.