Monogenean anchor morphometry: systematic value, phylogenetic signal, and evolution

Background. Anchors are one of the important attachment appendages for monogenean parasites. Common descent and evolutionary processes have left their mark on anchor morphometry, in the form of patterns of shape and size variation useful for systematic and evolutionary studies. When combined with morphological and molecular data, analysis of anchor morphometry can potentially answer a wide range of biological questions. Materials and Methods. We used data from anchor morphometry, body size and morphology of 13 Ligophorus (Monogenea: Ancyrocephalidae) species infecting two marine mugilid (Teleostei: Mugilidae) fish hosts: Moolgarda buchanani (Bleeker) and Liza subviridis (Valenciennes) from Malaysia. Anchor shape and size data (n = 530) were generated using methods of geometric morphometrics. We used 28S rRNA, 18S rRNA, and ITS1 sequence data to infer a maximum likelihood phylogeny. We discriminated species using principal component and cluster analysis of shape data. Adams’s Kmult was used to detect phylogenetic signal in anchor shape. Phylogeny-correlated size and shape changes were investigated using continuous character mapping and directional statistics, respectively. We assessed morphological constraints in anchor morphometry using phylogenetic regression of anchor shape against body size and anchor size. Anchor morphological integration was studied using partial least squares method. The association between copulatory organ morphology and anchor shape and size in phylomorphospace was used to test the Rohde-Hobbs hypothesis. We created monogeneaGM, a new R package that integrates analyses of monogenean anchor geometric morphometric data with morphological and phylogenetic data. Results. We discriminated 12 of the 13 Ligophorus species using anchor shape data. Significant phylogenetic signal was detected in anchor shape. Thus, we discovered new morphological characters based on anchor shaft shape, the length between the inner root point and the outer root point, and the length between the inner root point and the dent point. The species on M. buchanani evolved larger, more robust anchors; those on L. subviridis evolved smaller, more delicate anchors. Anchor shape and size were significantly correlated, suggesting constraints in anchor evolution. Tight integration between the root and the point compartments within anchors confirms the anchor as a single, fully integrated module. The correlation between male copulatory organ morphology and size with anchor shape was consistent with predictions from the Rohde-Hobbs hypothesis. Conclusions. Monogenean anchors are tightly integrated structures, and their shape variation correlates strongly with phylogeny, thus underscoring their value for systematic and evolutionary biology studies. Our MonogeneaGM R package provides tools for researchers to mine biological insights from geometric morphometric data of speciose monogenean genera.


INTRODUCTION
The Monogenea is a class of flatworms (Platyhelminthes) that are primarily ectoparasites of fish (Whittington, 2005;Hayward, 2005). An adult monogenean parasite has well-developed attachment appendages located at its anterior (prohaptor) and posterior (opisthaptor) regions that help it to resist physical dislodgement from the host. The posterior attachment organs consist of sclerotized hard parts such as hooks, anchors and clamps. Ecologically, monogenean parasites are characterized by their strong host specificity (Whittington et al., 2000). The Monogenea has several desirable features that make it invaluable as a model system for studying evolutionary processes that resulted in its past diversification and Manuscript to be reviewed present diversity (Poulin, 2002). Primarily, many of its genera are speciose, morphologically diverse, show well-resolved phylogenies at the familial level Kritsky, 1997, 2001;Mollaret et al., 2000), and samples can be easily obtained in large numbers. It has been used as a model to shed light on ecological forces that shape species community and structure (Rohde, 1979;Mouillot et al., 2005;Raeymaekers et al., 2008), to investigate processes leading to speciation and its maintenance (Rohde and Hobbs, 1986;Rohde 1994;de Meeus et al., 1998;Šimková et al., 2002;Hahn et al., 2015;Vanhove and Huyse, 2015), to elucidate host-parasite evolutionary ecology Huyse et al., 2003;Huyse and Volckaert, 2005;Šimková et al., 2006;Šimková and Morand, 2008;Mendlová and Šimková, 2014;Grégoir et al., 2015), and to explore the extent of correlation between phenotype variation in attachment organs and factors such as phylogeny, host specificity and geographical location .
Morphometric variation in anatomical structures of interest can be studied using two approaches. Traditional morphometrics (Reyment et al., 1984;Marcus, 1990) is characterized by the use of lengths of defined positions on anatomical structures of interest (or their ratios) as input data for multivariate statistical analyses. While such variables may measure size adequately, they are generally not effective for capturing shape information present in the geometry of a set of defined points of an object (Rohlf and Marcus, 1993). A large proportion of biological variation due to shape differences is therefore missed when an analysis uses only information from variation in length variables.
With the development of geometric morphometrics over the past three decades, researchers now have, at their disposal, a powerful method for extracting, visualizing and combining shape data with other data types such as molecular phylogenies to attain an integrative evolutionary analysis (Rohlf and Marcus, 1993;Adams et al., 2004;. Digitization of the anatomical structure of interest provides the key to the acquisition and use of a new type of data -landmark coordinates, from which shape information can be effectively extracted, and then analyzed, using new tools such as Procrustes superimposition, thin plate splines, relative warp analysis and elliptic Fourier analysis. Geometric morphometrics is now commonly used in systematics and evolutionary biology research where analysis of shape can be expected to provide new insights to complement traditional morphometric, phylogenetic or biogeographic analyses. A cursory search in major biological journal databases for recent publications having "geometric morphometrics" in their keywords revealed that geometric morphometrics is widely used to study various biological aspects, in diverse phyla, such as fish taxonomy (Sidlauskas et al., 2011), plant taxonomy (Conesa et al., 2012), gastropod shell shape variation (Smith and Hendricks, 2013;Cruz et al., 2012), morphological adaptation in birds (Sievwright and MacLeod, 2012), fly wing evolution (Pepinelli et al., 2013), turtle neck shape evolution (Werneburg et al., 2015), beetle speciation (Pizzo et al., 2013) and species boundary problems in butterflies (Barão et al., 2014). Because of the inherently digital nature of geometric morphometric data, its increasing prominence in morphological studies accentuates the role of informatics in modern taxonomy (Wheeler, 2007).
In morphological analyses of monogeneans, taxonomists often prioritize prominent sclerotized parts such as the copulatory organ, because qualitative variation in the latter is frequently sharp and easy to describe. Nonetheless, morphometric variation in all sclerotized parts of monogeneans has been studied for a long time from the perspective of systematics (e.g. Shinn et al., 2001) and evolutionary ecology (e.g. Poisot and Desdevises, 2010;Mendlová and Šimková, 2014). Hard parts such as the anchors are ideal for geometric morphometric analysis because they are not easily deformed by compression when mounted onto slides (Lim and Gibson, 2009). Anchor shape and size are taxonomically informative. Typically, size information is captured quantitatively in the form of distances between two defined -points on the anchors, and shape information is captured qualitatively in the form of character states. The analysis of monogenean morphometric data has been, and continues to be, dominated by the application of Manuscript to be reviewed traditional morphometrics (e.g. Mariniello et al., 2004;Shinn et al., 2004;Tan et al., 2010;Hahn et al., 2011;Soo and Lim, 2012). To date, there are only few examples (Vignon and Sasal, 2010;Vignon, 2011;Vignon et al., 2011;Rodríguez-González et al., 2015a) of applying geometric morphometrics to analyze monogenean anchor shape variation to overcome the limitations in data resolution inherent in standard morphometric and qualitative analyses. The paucity of geometric morphometric studies, however, belies the importance of this approach in uncovering intraspecific shape variation in anchors that can be invaluable for species delimitation, particularly in resolving synonymyiess (e.g. Pérez Ben et al., 2014), as well as for testing hypotheses of morphological integration (Olson and Miller, 1958) and evaluating levels of phenotypic plasticity (Pfennig et al., 2010).
As anchors serve a functional purpose, a priori, it is unclear whether phenotypic similarity of anchors among species is an outcome of adaptive processes related to the ecology or the morphology of their fish host, or simply a reflection of their phylogenetic constraint . If the presence of phylogenetic signal in anchors can be statistically established, evolutionary analysis of shape and size change can then be used to elucidate trends in particular clades. The results are expected to be useful for guiding the selection of appropriate anchor morphometric variables for conversion into morphological characters that have lower levels of homoplasy, thus overcoming the problem of unnecessary homoplasy of a morphological character arising from poor quality and insufficient number of character states (Perkins et al., 2009).
In this paper, we developed an integrative analysis that uses data from anchor morphometry and morphology, as well as DNA sequences, that allows the investigation of broad aspects in the systematic biology of monogeneans, such as species discrimination, evolutionary ecology, phylogenetic signal, and morphological integration. For illustration, we used data obtained from 13 recently described species belonging to the Ligophorus (Monogenea: Ancyrocephalidae) genus, a particularly speciose genus with 63 species known to date (Table   S1). Thus, our study covered all the Ligophorus species (approximately 20% of global Ligophorus diversity) currently found on two mullet host species in Southeast Asia. Mullets are eaten by people in this region.
Specimens were initially mounted in modified ammonium picrate glycerine, and subsequently converted into unstained, permanent mounts in Canada balsam. The opisthaptorial sclerotized hard parts of Ligophorus consist of a pair (left and right) of dorsal and ventral anchors, bars, and marginal hooks. Digital images of these hard parts were taken from labelled mounted slides using a light microscope with Leica digital camera (DFC 320) connected to the QWin plus image analysis software (Leica Microsystems, Germany) under 40x magnification, saved as jpeg files and organized into folders. Three species (L. fenestrum, L. liewi and L. kederai) showed a probable fixed character state -that of the presence of fenestrated structures on anchors of all examined specimens, since no variation in this character state was observed in all examined specimens of these three species.

Data Acquisition
Where possible, we measured a specimen's body length from its anterior to posterior end, and body width at the midpoint of its body. To obtain landmark coordinate data from the anchors, we used TPSDIG2 (Rohlf, 2013;Rohlf, 2015). Eleven landmarks (LM; six Type I, five Type III) were placed sequentially on the right and left ventral and dorsal anchors of each specimen ( Fig.   1). The set of all 11 landmark coordinates makes up a specimen's landmark configuration.
Six of the landmarks are Type I (LM1,LM2,LM3,LM5,LM7,LM8), while the remainder (LM4,LM6,LM9,LM10,LM11) are Type III (i.e. semi-landmarks). LM1and LM3 are the inner root point and the outer root point, respectively. Sandwiched between them is LM2, the groove point. LM5 is the dent point, while LM7 is the curve point. The tip point is represented by LM8.
The semi-landmarks were defined relative to Type I landmarks. The horizontal (towards the outer root point) and vertical projections (towards the curve point) from LM2 intersect with the anchor outline to give LM4 and LM6, respectively. LM9 and LM10 are the intersection points between the vertical projection from LM7 and LM1 with anchor outline, respectively. The projection from LM2 perpendicular to the vertical projection from LM1 touches the anchor outline to define LM11. We used the set of landmarks LM1 to LM4 and LM11 to represent the shape of the root compartment, and the set LM5 to LM10 to represent the point compartment. For geometric morphometric analysis, semi-landmarks were not specially treated (e.g. employing sliding landmark analysis), following MacLeod (2013) that such treatments may introduce distortions to the original geometrical relationship that lead to complicated interpretations of the result. The anchor images and their corresponding landmark coordinate data are given in Supplemental Material D1 and D2, respectively.
The monogeneaGM package contains a suite of R functions for three primary analyses: anchor landmark coordinate data quality control, visual checks and Generalized Procrustes Analysis; hierarchical clustering and principal component analysis; and analysis of anchor shape change using directional statistics. Some of the functions are suitable for general use. For example, the phylomorphospace visualization functions -tpColorPlot2d for twodimensional data, and tpColorPlot3d for three-dimensional data, can take in two additional arguments: a color transparency control, which is useful for improving graphical presentation when there are substantial overlaps of data points, and an option to superimpose a user-supplied phylogeny onto defined centroids in the scatter plot, if this information is available.

Manuscript to be reviewed
Despite careful slide preparation, it is inevitable that anchor images of some specimens would contain substantial amount of non-biological shape variation caused by incongruent image and object planes (Arnqvist and Mårtensson, 1998). The inclusion of these poor quality data in downstream analyses is undesirable, as they introduce noise into an analysis that can potentially complicate the interpretation of results. To mitigate this problem, we developed a quality control procedure to filter out poor quality images. In this procedure, we first computed all pairwise Euclidean distances between landmarks for the left and right forms of dorsal and ventral anchors. . The magnitude of this measure is straightforward to interpret -it is high (maximum 100) for good quality specimens and low (minimum 0) for poor ones. Figure 2 shows examples of poor and good quality specimens together with their Tukey Mean-Difference (TMD) plots, respectively. Specimens with Q of 10 or more (n=437 ; Table S3

Converting Pairwise Euclidean Distances in Arbitrary Units to Physical Units
We used a subset (n=96) of the total specimens with quality score above 10 (n=437) and measured the physical distances from LM1 to LM3 and from LM1 to LM5 in these samples using QWin plus image analysis software (Leica Microsystems, Germany). We then regressed the Manuscript to be reviewed physical distances against the computed pairwise Euclidean distances to determine the linear equation for converting arbitrary distance units into their physical units (in μm). Thus, all pairwise Euclidean distances computed from raw landmark coordinates could be converted to physical distances by multiplication with a factor of 0.2 followed by addition of 0.9 (Fig. S3).

Geometric Morphometric Analysis
For each species, we performed Generalized Procrustes Analysis (GPA; Gower, 1975;Rohlf and Slice, 1990) to align the sample landmark configurations for both ventral and dorsal anchors, using the gpagen function in the geomorph package (Version 2.1.1; Adams and Otarola-Castillo, 2014). The resulting GPA coordinates of the left and right forms were then averaged. GPA removes the effects of translation, rotation and scaling so that the resulting landmark configurations have minimum sum of squared distances with respect to the mean landmark configuration (Adams et al., 2004). Nevertheless, even after GPA, comparison of anchor shape variation can still be potentially confounded by the presence of non-biological variation in the landmark configuration. Specifically, if many samples of a species have anchors lying in one particular position, it would not be clear whether variation between its members' mean GPA landmark configuration and those of other species constitutes genuine biological variation or mathematical artifact. Typical application of geometric morphometrics in nonmicroscopic objects (fly wings, skulls, etc.) does not usually suffer from this problem, since specimens from species to be compared can be manipulated into standardized positions before imaging.
To ensure that the landmark configurations of all 13 species were comparable, we determined the angular deviation of LM7 from the x=0 line, and rotated all landmark coordinates by this amount with the origin as pivot. This has the effect of creating standardized landmark configurations for specimens across all species since the x-coordinate of LM7 is always zero after adjustment. The GPA coordinate data of all specimens thus obtained were subjected to another iteration of GPA to produce the final shape alignment, which was organized into a data matrix with rows representing specimens and columns representing the 44 GPA landmark coordinates.

Molecular Phylogenetic Analysis
We used DNA sequence data from three nuclear markers: 28S rRNA, 18S rRNA and ITS1 from the 13 Ligophorus species to infer their phylogenetic tree. These markers are a mainstay of parasitic platyhelminth molecular phylogenetics (Littlewood, 2008), and the combination of fast (ITS-1) and slow (28S rRNA, 18S rRNA) evolving sequences should provide sufficient resolution for the inferred phylogenetic tree (Lockyer et al., 2003;Blair, 2006;Waeschenbach et al., 2007), allowing the relationship of even closely related monogenean species to be resolved (e.g. in Gyrodactylus (Gilmore et al., 2012)). Ideally, the inclusion of mitochondrial markers such as cytochrome oxidase I (COI) would provide a more representative sampling of the genome, and hence, a more reliable phylogeny. However, the absence of mitochondrial data may not impact the quality of the inferred phylogeny, since, for monogeneans, it has been shown in Diplectanidae and Gyrodactylus that both rRNA and mitochondrial markers are effective for species identification . Regardless of whether nuclear or mitochondrial genes are used, it is important to keep in mind that gene trees do not always reflect the species tree (Maddison, 1997).
Partial 28S rRNA (~800bp) and ITS1 (~750bp) sequence data were obtained from , whereas 18S rRNA sequence data were generated in the present study. Briefly, the Ligophorus specimens were removed from host gills, identified morphologically and then preserved in 75% ethanol. Genomic DNA was extracted from samples using DNAEasy extraction kit (QIAGEN, Hilden, Germany). About 5µl of the extracted DNA was used as template in the PCR reaction to amplify the partial 18S rRNA sequence using two primers: WormA (5'- GCGAATGGCTCATTAAATCAG-3') (Littlewood and Olson, 2001) and new930F (5'-CCTATTCCATTATTCCATGC-3') (modified from Littlewood and Olson, 2001). The PCR reaction (50µl) was carried out in a solution containing 1.5mM MgCl2, PCR buffer (Fermentas), 200µM of each deoxyribonucleotide triphosphate, 1.0 µM of each PCR primer and 1U of Taq polymerase (Fermentas), in a thermocycler (Eppendorf Mastercycler) using the following conditions: initial denaturation at 95°C for 4 minutes, followed by 35 cycles of 95°C, 52°C and 72°C for one minute each, with final extension at 72°C for 10 minutes. An aliquot (10µl) from the amplicons were electrophoresed in 1.3% agarose gel, stained with ethidium bromide and viewed under an ultraviolet illuminator. The remaining 40µl of each amplicon was purified using a DNA purification kit (QIAGEN, Hilden, Germany) and subjected to automated DNA sequencing (ABI 3730 DNA Sequencer, First Base Laboratories, Kuala Lumpur) using the same primers used for PCR amplification. Approximately 750 bp of the 18S rRNA sequence were amplified and sequenced for the 13 Ligophorus species (Table 1).

Ligophorus species
Host species Locality (Malaysia) GenBank Accession no.  For phylogenetic analysis, we first aligned sequences for each marker using MAFFT (Version 7 at http://mafft.cbrc.jp/alignment/server/; Katoh and Standley, 2013). The alignment parameters used were the Q-INS-i iterative refinement method, the 1PAM/=2 nucleotide scoring matrix with a gap opening penalty of 1.53. We then concatenated multiple sequence alignments of the three nuclear markers (28S rRNA-ITS1-18S rRNA). The resulting supermatrix was used as input in IQ-TREE (Nguyen et al. 2015) to construct the maximum likelihood (ML) phylogenetic tree (Felsenstein, 1981;Felsenstein, 2003). IQ-TREE is a state-of-the-art ML tree construction pipeline that integrates DNA model selection, maximum likelihood tree search, and rapid bootstrap analysis (Minh et al., 2013). For model selection, we used the Bayesian Information Criteria, and did not consider "G+I" models, following Yang (2006) that modelling proportion of invariable sites in the presence of gamma rate variation across sites creates model identifiability issues. MEGA (Version 6; Tamura et al., 2013) was used to edit the resulting phylogenetic tree.
We annotated the tree with the morphology of anchors, bars and male copulatory organ to allow visual assessment of overall phylogenetic and phenotypic correlation.

Species Discrimination
Discriminating a monogenean species is a complex art that involves the comparison of qualitative features of numerous anatomical structures: the male copulatory organ, female reproductive organ, anchors, bars and marginal hooks. Among the sclerotized hard parts, multivariate morphometric analyses of shape and size variables of suitable anatomical structures provide a quantitative means for species discrimination, which is invaluable for complementing the results from qualitative morphological analyses.
To visualize species clustering in low dimension morphospace, we applied Principal Component Analysis (PCA) separately for the ventral and dorsal anchors using their GPA coordinate data. The trade-off between loss of information through dimensional reduction and gain of interpretation via visualization in PCA can, however, make it difficult to judge how well members of the same species cluster together in the PCA scatter plots, especially when there are overlaps between different species clusters. To overcome this problem, we complemented PCA results with the cluster heat map (Wilkinson and Friendly, 2008), a powerful method for organizing high-dimensional multivariate data that allows visual detection of patterns of variation. The cluster heat map first maps numerical information in the cells of the input data matrix to corresponding color codes. Then, a hierarchical clustering algorithm is applied to cluster the samples by similarity, in such a way that within cluster variation is always smaller than between cluster variation. For the current analysis, we estimated similarity between each pair of samples using the Manhattan distance metric. The resulting distance matrix was then used as input for hierarchical clustering of samples using the Ward algorithm.
To assess the impact of applying the quality control procedure, we compared cluster heat maps generated using all samples, and using only samples that passed data quality control. Heat map construction was done using the heatmap.2 function in the gplots package (Version Manuscript to be reviewed 2.13.0; Warnes et al., 2014). We found the simple heat map a good alternative to inspection of the PC loadings table when trying to interpret the first few PC axes biologically.

Testing for the Presence of Phylogenetic Signal in Anchor Shape
Species with different shapes are localized in particular regions of the morphospace.
When a phylogeny is superimposed onto this morphospace, a phylomorphospace is induced, and it becomes possible to evaluate whether common descent or convergent evolution is likely to have shaped phenotypic similarity (Klingenberg and Ekau, 1996;Sidlauskas, 2008;Revell, 2014). If anchor shape contains substantial phylogenetic signal, then we expect the phylogeny to We formally tested the presence of phylogenetic signal in anchor shape by applying Adams's multivariate Kmult test (Adams, 2014a), implemented using the physignal function (10 000 iterations) in the geomorph package (Version 2.1.1; Adams and Otarola-Castillo, 2013).
The Kmult statistic is a multivariate generalization of Blomberg's K statistic (Blomberg et al., 2003). The phenotype (a continuous trait) interest in the lineages of a given phylogeny is assumed to evolve in phylomorphospace according to Brownian motion. In the absence of phylogenetic signal, Kmult = 0 (i.e. phenotypic variation is independent of lineages). At Kmult =1, phenotypic variation between taxa in the same lineage conforms to expectation under Brownian motion evolution. Values of Kmult < 1 correspond to phenotypic variation that is larger than Manuscript to be reviewed expected between taxa of the same lineage, and values of Kmult >1 correspond to phenotypic variation that is smaller than expected between taxa of the same lineage (Adams, 2014a).
Statistical significance is determined via a permutation procedure under assumption of phylogenetic signal absence.

Analysis of Anchor Shape and Size Evolution
In studying anchor shape and size evolution, we were primarily concerned with trends occurring in different clades of the phylogeny of the 13 Ligophorus species. To control for the effect of body size in subsequent phylogenetic regression analysis of anchor shape against anchor size, it was necessary to first test for collinearity of body size and anchor size (Mundry, 2014).
Since body size was prone to distortion during fixation, we used the median of body length and body width of each species to reduce the impact of outliers. For analysis, the logarithm (base 10) of the product of median body length and median body width was used.
The GPA landmark coordinates of the ancestral anchor were estimated using the maximum likelihood method as implemented in fastAnc function from the phytools package (Version 0.4-21; Revell, 2012) . Anchor shape change associated with a clade is statistically supported if mean directional change deviates significantly from uniformity in a set of landmarks. We visualized directional deviation in the 11 landmarks of both ventral and dorsal anchors using circular plots (Agostinelli and Lund, 2013; implemented in the circular package, Version 0.4-7). We then performed Rayleigh's test (Batschelet, 1981) to test for evidence against directional uniformity in each landmark. The strength of statistical evidence against mean directional uniformity in each landmark was assessed using p-value. Wireframelollipop plots (Klingenberg, 2013) were used to graphically summarize the mean change in direction and mean magnitude of landmark displacement from root ancestor landmark configuration. For both regression analyses, p-values were computed via a resampling procedure with 10 000 iterations.

Covariation of Anchor Shape and Size with Male Copulatory Organ Morphology
Rohde and Hobbs (1986) hypothesized that the reproductive barrier among congeneric species that share the same host can be maintained in monogenean parasites by their having different copulatory organ morphology when attachment organs are similar (thus occupying Manuscript to be reviewed similar microhabitats); conversely, when parasites have dissimilar attachment organs (thus occupying different microhabitats), the morphology of their copulatory organs would not show important differences, since the lack of proximity puts less evolutionary pressure on the parasites to evolve morphologically different copulatory organ. Qualitative evidence with limited number of congeneric species (Lambert and Maillard, 1975;Roubal, 1981;Rohde et al., 1994) supported the hypothesis's feasibility, as well as later studies with more species Morand, 2008, 2015). Quantitative evaluations using larger species assemblage that relied on traditional morphometric data are available, but the interpretation of their results in support of the hypothesis was obscured by either the problem of using inflated degrees of freedom in regression analysis (e.g. Šimková et al., 2002) or failure to control for the effect of phylogeny (e.g. Jarkovský et al., 2004). With the development of new tools for geometric morphometric and phylogenetic comparative methods, we are in a position to retest the Rohde-Hobbs hypothesis. To this end, we compared the size of the male copulatory organ (mean tube length data from Lim, 2012, 2015; and three of its selected morphological characters (Table 2: position of copulatory organ entrance at the main lobe of accesory piece; accesory piece of male copulatory complex; shape of accesory piece of male copulatory complex) against anchor shape and size variation. Ancestral node positions were estimated as before using the fastAnc function in the phytools package. (2) outer root longer than point X -√

Morphological Integration Analysis
The roots of the anchor are bases for muscle attachment. Biomechanically, force exerted through the muscles and transmitted to the point compartment controls the anchor's grip strength on the gills. Because of this, we may expect the anchor to function as a single, fully integrated module (Klingenberg, 2008) on a priori grounds. Anchor shape is strongly constrained by either phylogeny or convergent evolution. If the latter's effect was weak, then suitable morphological characters based on variation in anchor shape can be expected to be systematically useful.
To date, only few morphological integration analyses in monogeneans have been done. analysis (Rohlf and Corti, 2000). Here, we extended their morphological integration analysis to the interspecific level in Ligophorus. We applied the phylogeny-aware partial least squares method based on the evolutionary covariance matrix (Adams and Felice, 2014) to estimate the extent of morphological integration between the ventral and dorsal anchors, as well as that of the root compartment (L1 to L4 and L11) and the point compartment (L5 to L10) within and between the ventral and dorsal anchors.

New Morphological Characters from Morphometric Variables
A continuous morphometric variable can be discretized and treated as a morphological character with two or more states for use in a cladistic analysis (Thiele, 1993;Rae, 1998;Wiens, 2001). In doing so, the taxonomist relies on experience and intuition to select promising   Table S4 for character state matrix) was chosen. We replaced the morphological characters derived from traditional morphometric measurements of anchors with new candidates derived from geometric morphometric analysis to assess their phylogenetic informativeness. To this end, we compared how well-resolved the resulting maximum parsimony trees (using PAUP; Swofford, 2002) were. Tree search (initial tree obtained via stepwise addition) was performed using the heuristic search option. Branch-swapping was done using the tree bisection and reconnection algorithm. Tree reliability was assessed using 1000 bootstrap replicates and branches were collapsed if bootstrap support was below 50%.

Molecular Phylogeny
The GTR + G model was the best DNA substitution model that did not incorporate the proportion of invariable sites. After mid-point rooting, the estimated ML tree ( Fig. 3; 10000 bootstrap replicates) contained two major clades. One of them consisted of species infecting M.
buchanani (Clade I), and the other consisted of species infecting Liza subviridis (Clade II).
Bootstrap support was high for most internal nodes, except the most recent common ancestor node of L. parvicopulatrix-L. bantingensis, and L. grandis-L. fenestrum (between 50-60%).  Table S5 gives summary statistics of anchor size, anchor shape, body size and male copulatory organ size for the 13 Ligophorus species.

Manuscript to be reviewed
Scatter plots of GPA landmark configuration of all specimens for the ventral and dorsal anchors are given in Fig. 4 (see Fig. S4 to Fig. S16 for species-specific alignments). Figure 5 shows the PCA plots of PC2 against PC1, and PC3 against PC1 for shape variables of ventral and dorsal anchors (See Fig. S17 for a three-dimensional PCA plot). The first three PC accounted for 84% and 79% of total shape variation in the ventral and dorsal anchors, respectively. To interpret these three PCs, we simultaneously compared the scatter plots of the GPA landmark configurations with the heat map of shape variable loadings ( Fig. S18 and S19

Cluster Analysis of Geometric Morphometric Data
The cluster heat map (Fig. 6) shows that variation in anchor shape alone allows the samples to be clustered unambiguously into 12 clusters corresponding to 12 of the Ligophorus species, confirming that between species variation is much larger than within species variation.
On the other hand, with only eight specimens used, the clustering outcome was ambiguous for L.

FIGURE 6.
Consistent with the detection of significant phylogenetic signal in anchor shape, hierarchical clustering revealed two major clades whose members were exactly the same as those of Clade I and Clade II. The quality of clustering using specimens that passed quality control was improved especially for species with anchors that have very similar shapes such as L. Manuscript to be reviewed navjotsodhii and L. chelatus (Misclassification error of 2.5% with filtering (Table S6); 3.8% with no filtering (Table S7)).
For each shape variable, we labelled the samples according to their membership in Clade I or Clade II, and then ranked the shape variables in descending order using the two-sample t-

Anchor Shape and Size Evolution
The wireframe-lollipop graphs (Fig. 7)  Manuscript to be reviewed distance is expected to be of value for discriminating the two clades. To be easy to measure, the landmarks should be of Type I. Thus, the distance from LM1 to LM3 and from LM1 to LM5 were found to be good candidates. We found the common practice of using the inner root length (IL) and the outer root length (OL) (distance from LM1 to LM7 and LM3 to LM7, respectively) to be suboptimal since both LM3 and LM7 had almost parallel directional changes (Figs. 7a, 7c, 7d), whereas mean directional change in LM1 and LM7 did not show any clear patterns of divergence in one clade and convergence in another (Fig. 7) to be able to show large variation between Clade I and Clade II. Figure 8 (see also Fig. S22 and Fig. S23) shows that it is possible to define cut-offs for the LM1-LM3 (15μm) and LM1-LM5 (25μm) distances that result in discrimination of Clade I from Clade II, but no reasonable cut-offs for IL and OL lead to similar results. Results from the Adams-Collyer phylogenetic regression indicated that the interaction of body size and anchor size was not statistically significant in the dorsal anchor (p-value = 0.2), but anchor size was a significant predictor of anchor shape (p-value = 0.01). For the ventral anchor, the interaction of body size and anchor size was a significant predictor of anchor shape (p-value = 0.02). Since body size and anchor shape were not significantly correlated, we may expect similar anchor shape to be found across a range of body sizes (Fig. S24).

Patterns of Morphometric and Morphological Variation in the Male Copulatory Organ and Anchor
Where the male copulatory organ was similar in size among closely-related species with similar anchor shape and size, its morphology varied (I and II in Fig.10). Ligophorus belanaki, L.
careyensis, L. navjotsodhii and L. chelatus shared a most recent common ancestor, whose ancestral character states for the three male copulatory organ characters could be inferred as 113 on a parsimony criterion. The divergence of L. careyensis from L. belanaki did not involve major Manuscript to be reviewed changes in anchor shape, size and size of male copulatory organ, but on the latter's morphology, which acquired three changes to become 002. Similarly, The most recent common ancestor of L.
navjotsodhii and L. chelatus probably evolved character states 000 or 001 from 113, and divergence of these two species was associated with a change in the third character state, with only relatively minor change in either anchor shape and size or copulatory organ size.
In contrast, where the male copulatory organ was similar in morphology among closelyrelated species with similar anchor shape and size, its size varied (III in Fig. 10). Ligophorus kederai, L. grandis, L. kedahensis and L. johorensis have similar anchor shape and the same character states 114 for the morphology of their male copulatory organ. Consistent with prediction from the Rohde-Hobbs hypothesis, substantial variation in the size of their male copulatory organ was observed. It is possible for size and morphological variation to co-occur in the male copulatory organ, as shown in the divergence of L. grandis and L. fenestrum from their common ancestor. Finally, species with similar anchor shape (IV in Fig. 10) showed large variation in both their male copulatory organ size and morphology.

Morphological Integration
The shapes of both ventral and dorsal anchors were strongly and significantly correlated Manuscript to be reviewed 0.03). Figure S25 provides a graphical summary of the results of the morphological integration analysis.

Phylogenetically Informative Morphometric Variables
For the current 13 Ligophorus species, the maximum parsimony tree estimated using the set of morphological characters containing discretized LM1-LM3 and LM1-LM5 distances and anchor shape was better resolved (Fig. 11) . Clade I and Clade II were clearly identified, and the partially resolved relationships within each clade were also congruent with those of the molecular phylogeny's. In contrast, using morphological characters of the anchors derived from traditional morphometrics as in Sarabeev and Desdevises (2014) produced a maximum parsimony tree that was mostly reticulate and unable to distinguish Clade I and Clade II. Manuscript to be reviewed

Data Quality Control
We are not aware of any geometric morphometric analyses of anchors in monogeneans that currently implement specimen quality control. Specimen quality introduces an important source of non-biological variation into observed anchor shape variation, the impact of which depends on whether the data would be analyzed at the intra or interspecific level. Thus, while the inclusion of specimens that failed quality control into hierarchical clustering did not fundamentally change species discrimination outcome in this study, it is important to control for this confounder where intraspecific variation can be expected to impact conclusions of an analysis, for example, when investigating mean directional change in landmarks of anchors ( Fig.   7), or testing for association between intraspecific anchor shape variation and evolutionary potential of a species (Rodríguez-González et al., 2015a) . In the current study, we observed up to 50% loss of specimens (L. fenestrum) due to low quality score. Assuming this value optimistically as an upper bound, then at least 40 specimens per species would have to be obtained in order to anticipate 20 or more specimens passing quality control. An ideal case like this may not always be possible, since sampling trips do not always yield sufficient study material.
The assumption of symmetry in the left and right side of the ventral and dorsal anchors used in the quality control procedure has a caveat that should be noted. Poor quality specimens that arise as a consequence of incongruent image and object planes are confounded with specimens that show fluctuating asymmetry. If fluctuating asymmetry (Graham et al., 2010) is common, then a large proportion of specimens would have been discarded using the quality control procedure, but this is not the case in the present study, where about 82% ( 437/530) of all specimens were retained. Careful observations over multiple well-prepared mounts should help Manuscript to be reviewed researchers decide the appropriateness of implementing the quality control step for specimens from a particular monogenean genus.
Variation in quality scores can be attributed to several sources, such as the method of slide preparation, the quality of camera lens and software used for capturing images, and the skill and experience of the data gatherer. In this study, a single data gatherer (O.Y.M. Soo) prepared and acquired the landmark data, using the same compound microscope and computer. Because of this, we expect other factors to explain the poor quality scores. Interestingly, we note that species that had larger proportion of specimens failing quality control tend to have body and/or anchor size that were relatively large or small. In the case of L. grandis (log10 median body size 2.1 SD larger from total mean; dorsal anchor size 1.0 SD larger than total mean), L. fenestrum (log10 median body size 1.6 SD larger than total mean; dorsal anchor size 0.9 SD larger than total mean), and L.
bantingensis (log10 median body size 0.9 SD smaller than total mean; ventral anchor size 2.7 SD smaller than total mean), we observed about 30%, 50% and 45% of the specimens failing quality control (Q-score < 10; Fig. S2), respectively. A possible explanation may be that the robust anchors of L. fenestrum and L. grandis have uneven thickness at the root and point regions, which makes them difficult to evenly flatten on slides. The large body bulk may also further hinder effective flattening. Monogenean anchor thickness is not usually measured but may be indirectly inferred through 3D-modelling (Teo et al., 2010;Teo et al., 2013). Since size and physical inertia are positvely correlated, the small body and anchor size of L. bantingensis make specimen orientation on the slide sensitive to variation in force applied during slide flattening. Manuscript to be reviewed Different species were found to have varying levels of intraspecific phenotypic plasticity in this study. While within species shape variation in both ventral and dorsal anchors was large in some species (L. kedahensis, L. parvicopulatrix), it was limited in others (L. grandis, L. liewi).

Integrative Geometric Morphometric Analysis Supports the Rohde-Hobbs Hypothesis
Evidence for supporting the Rohde-Hobbs hypothesis has traditionally come from integrating spatial distribution data of monogeneans on gill microhabitats (e.g. Rohde, 1977;Ramasamy et al., 1985;Koskivaara et al., 1992) with morphological data of the monogenean species (e.g. Fig. 6.3 in Šimková and Rohde, 2013). These efforts were very laborious, but crucially established anchor shape-microhabitat association. Benefitting from such insights, our current integrative geometric morphometric analysis was able to reveal patterns consistent with the hypothesis's predictions on how male copulatory organ size and morphology vary with respect to anchor shape (Fig. 10), despite the absence of spatial distribution data for the 13 Ligophorus species across gill microhabitats.

Manuscript to be reviewed
In their intraspecific study of morphological integration between the root and point compartments in L. cephali, Rodríguez-González et al. (2015a) reported only modest degree of integration in the same anchor, but stronger compartmental integration between the ventral and dorsal anchors. On the other hand, using interspecific data, we demonstrated a much stronger degree of integration between the root and the point compartments within anchors, and showed relatively weaker integration of the root compartments between ventral and dorsal anchors.
Intuitively, intraspecific compartmental integration within the same anchor is expected to be high, so a possible explanation for the discrepancy may be the lack of a quality control procedure for filtering poor quality slides. Without the latter, it seems difficult to rule out the possibility that the observed intraspecific anchor shape variation in L. cephali may contain non-trivial amount of artifactual noise.
Generally, a certain degree of homoplasy may be expected in the morphology of attachment organs in parasites, on grounds that functional requirements for attaching to the host and adapting to within-host microhabitats would override shape constraints imposed by phylogeny . Such form of adaptive evolution can cause Kmult to become less than 1 (Blomberg et al. 2003). In the present Ligophorus phylogeny, the Kmult value of 0.948 is significantly larger than 0 (p-value 0.0003), but slightly biased downwards from the expected value under Brownian motion evolution for anchor shape. The reason seems to be that, although anchor shape is more or less lineage-dependent in the present set of Ligophorus species, three species: L. liewi (Clade I), L. funnelus (Clade II) and L. parvicopulatrix (Clade II), have similar shaft shape (Fig. 9), despite being from different lineages, thus lowering Kmult. To put the magnitude of Kmult in perspective, its univariate version -Blomberg's K statistic, has been found to be generally less than 1 for numerous primate traits, exceeding 1 only for brain size (Kamilar and Cooper, 2013). The phylogeny (Bakke et al., 2002;Vanhove et al., 2015), immunophysiology (Buchmann and Lindenstrøm, 2002), ecology (Šimková et al., 2006), and behaviour of the fish hosts Mendlová and Šimková, 2014 ) all contribute to determine patterns of infection and transmission, host specificity, and subsequent mode of speciation in their monogenean parasites (Littlewood et al., 1997;Cribb et al., 2002;Whittington and Kearn, 2011;Vanhove and Huyse, 2015). Therefore, it is not surprising that different monogenean families may show different degrees of morphology-phylogeny covariation. In Lamellodiscus (Family: Diplectanidae), the lack of cospeciation between host-parasite  suggests that adaptation to host after host-switching would lead to convergent evolution in some of the morphologies (but see Machkewskyi et al., 2014). Indeed, Poisot et al. (2011) found morphological features of sclerotized haptor and male copulatory organs to be weakly associated with phylogeny in Lamellodiscus. On the other hand, in Cichlidogyrus (Family: Ancyrocephalidae), sclerite (including anchors) shape was found to contain significant phylogenetic signal . However, this conclusion needs to be qualified in light of the recent demonstration of a host-switching event in Cichlidogyrus (Messu Mandeng et al., 2015), which was missed in Vignon et al. (2011) because the host range (at the familial level of the fish hosts) of Cichlidoygrus was not taken into consideration. For Ligophorus, the significant phylogenetic signal found in the anchors does not appear to need similar qualification, since its host range is currently known to be restricted to mullets. In an examination of over 1000 fish from 10 fish families in the Black Sea, Öztürk and Özer (2014) observed that Ligophorus was restricted to mullets (Mugil cephalus and Liza aurata), with almost 100% prevalence in the hosts.
Nevertheless, since it has been reported that a combination of host-switching (within fish host family) and intra-host speciation events probably generated the present diversity of Ligophorus Manuscript to be reviewed regarding phylogenetic signal in the anchors requires additional tests using samples of Ligophorus species from that region.
The anchor morphology-phylogeny correlation shown for Ligophorus in the present study suggests new morphological variables from the anchors that are potentially useful in morphometric analysis and also morphological phylogenetics. In particular, we suggest more active exploration of the newly proposed morphometric variables: the distance between the inner root point and the outer root point (LM1 -LM3), and the distance between the inner root point and the dent point (LM1 -LM5), as alternatives to existing usage of the inner root length (IL) and the outer root length (OL) in traditional morphometric analysis. The discretization of the two proposed morphometric variables: anchor shaft shape and anchor shaft shape, can also be further considered in future morphological phylogenetic analysis.

Anchor Shape and Size Correlation with Host and Ecological Factors
In the present study, the modest levels of anchor shape-size covariation revealed through the Adams-Collyer phylogenetic regression analysis for the 13 Ligophorus species suggest that, apart from the effect of shared ancestry, anchor shape-size covariation is likely non-trivially constrained by additional factors, one of which could be their biomechanical compatibility.
Another factor is host size and ecology. On average, the larger Moolgarda buchanani (body length range 35-48cm) harbored larger Ligophorus species, whereas the smaller Liza subviridis (body length range 25-30cm) harbored smaller Ligophorus species.
The current analysis showed that L. bantingensis evolved anchors of a smaller size but retained the sickle-shape shaft common in Clade I, which is associated with larger anchor size.
The observed small anchor size is consistent with the hypothesis that small or medium-sized attachment organs are generally associated with a generalist lifestyle, as these sizes expand the Manuscript to be reviewed host range of monogeneans to small or medium-sized hosts, which are generally more diversified than larger hosts .  (Koehn and Crook, 2013). Indeed, for skin-parasitic fish monogeneans, their attachment appendages must be able to withstand dislodgement by strong water currents pushing against the fish surface as the fish host performs burst swimming (Kearn, 2014).
Concurrently, it seems plausible that burst swimming may also, for a short time, cause strong water currents to flow through the gill chamber, thus imposing a selective pressure for more robust anchors in monogeneans that live on the gills.
The U-shaped root groove in some species infecting M. buchanani (e.g. L. liewi, L. kederai and L. fenestrum), which could result from the accretion of sclerotic material in the space between the inner and the outer root (Klaus Rohde, pers. comm., 2014), may have substantially expanded the root base, providing space for connection with more muscle tissues. This would likely have resulted in anchors with stronger contraction strength, necessitating the evolution of the shorter but more robust sickle-shaped shaft, hence the finding of tight morphological integration between the root and point compartments.
Interestingly, U-shaped roots and sickle-shaped shaft are also common in Cichlidogyrus species , which infect the cichlids in Africa. Ecologically, it seems unlikely that such robust shapes would have evolved in freshwater environments, where cichlids are mostly found. If the ancestors of cichlids and their monogenean parasites were marine in origin, or host switching occurred with contact between salt-tolerant cichlids and the marine ancestor of the monogenean parasite, the observed robust shapes would then be inherited, with their shapes constrained by phylogeny. While some studies (Murray, 2001) suggested that the ancestors of cichlids were likely marine, others were sceptical (Chakrabarty, 2004;Sparks and Smith, 2005).
Joint consideration of the morphology and phylogeny of monogenean fauna of cichlids and other marine fishes is probably required  to assess the competing claims. For example, a recent molecular phylogenetic analysis using 28S rRNA sequence (Tan, 2013;Fig. S26) indicated that Ligophorus (marine) and Cichlidogyrus (mostly freshwater) species shared the most recent common ancestor.
In three species (L. grandis, L. funnelus and L. liewi), fenestration in the anchor base seems to be an invariant character state, as all examined specimens (n=22,50,32, respectively) showed consistent presence of fenestration. Presently, the ecological significance of fenestrations in anchors is unclear. Some progress may be possible with biomechanical studies (e.g. Wong and Gorb, 2013) that compare whether fenestrated and non-fenestrated anchors differ significantly in their resistance against turbulence and strong water currents. The present phylogenetic analysis suggests that fenestration is not a synapomorphic character state (Fig. 3), but a clearer picture requires more extensive taxa sampling.

Outlook for Geometric Morphometric Analysis of Sclerotized Haptoral Structures in Monogenean Parasites
The use of less sophisticated methods could have also answered the biological problems that were studied in the present work. For example, to infer the presence of phylogenetic signal in attachment organs, or to investigate anchor size correlation with parasite and host phylogeny, one could map qualitative shape and size categories of the organ of interest onto the molecular phylogeny of the parasites, as done in Šimková et al. (2006). Validation of the Rohde-Hobbs Manuscript to be reviewed hypothesis could be done by correlating parasite abundance data in different gill microhabitats with their attachment organ and male copulatory organ morphology. The lack of a unified quantitative framework, however, forces the researcher to use different data types to address each question. We have shown in the present work that, if geometric morphometric data is used as a starting point, all the preceding questions can be simultaneously answered. In fact, the geometric morphometric approach is the key that unlocks the potential of shape information in the organs of interest for broader research questions, such as species and biogeographical population discrimination (Vignon and Sasal, 2010), morphological integration (e.g. Vignon et al. 2011;Rodríguez-González et al., 2015a), and more recently, canalization and developmental stability (Llopis-Belenguer et al. 2015), all of which cannot be answered satisfactorily using traditional morphometric or qualitative methods.
Our present study used large numbers of samples, averaging about 35 per species.
Compared to laborious measuring of selected lengths as done in traditional morphometrics, data acquisition is far more efficient with a landmark digitization software such as TPSDIG2, which simultaneously captures shape and size variation information. This improved efficiency is important -by greatly reducing the tedium associated with measuring many lengths per specimen, there is more incentive to sample more extensively.
Although the geometric morphometric approach has been strongly advocated by Vignon and Sasal (2010) as an effective means to pursue systematics research in monogeneans, the scientific community still lack integrative tools that would make it easy to share data and adopt a common analysis pipeline to ease comparison of old and new results. Here, the monogeneaGM R package that we have developed enables substantial number of shape and size variables from large number of samples to be analyzed efficiently to answer multiple questions, ranging from systematic value of the sclerotized haptoral structures to understanding patterns of phenotypic and phylogenetic correlation. We hope the development of monogeneaGM will contribute to reducing the bottleneck for large scale data analysis in the study of monogeneans. Indeed, as there has been a surge in alpha taxonomy and systematic biology studies of Ligophorus in recent years (Abdallah et al., 2009;Siquier and Ostrowski de Núñez, 2009;Blasco-Costa et al., 2012;Dmitrieva et al., 2012Dmitrieva et al., , 2013Lim, 2012, 2015;El Hafidi et al., 2013;Kritsky et al., 2013;Sarabeev and Desdevises, 2014;Marchiori et al., 2015;Rodríguez-González et al., 2015b;, the analysis can only get more interesting as data from other species from other hosts in other geographical regions are added. We expect the data analysis tools in monogeneaGM to continue to evolve to handle complexities of data analysis when this happens. The use of two-dimensional landmark data implies that the analysis of anchor size and shape evolution is necessarily approximate, since some of the potential biological variation in anchor morphometry may only be adequately captured in three dimensions (Galli et al., 2007).
Nevertheless, given the wealth of corroborative inference regarding anchor shape and size evolution that have been obtained in the current study, it appears that no general loss of interpretability arises from usage of two-dimensional data for geometric morphometric analysis in Ligophorus.

Future Prospects
Adequate taxa coverage remains an important factor for accurate phylogenetic inference (Sanderson et al., 2010). The three markers used in this study are the most common ones reported for other known Ligophorus species in the GenBank database, hence their continued use supports efforts at expanding taxa sampling of molecular sequences. Indeed, a research program in Ligophorus systematics that expands taxa coverage of anchor geometric morphometric and sequence data opens up the possibility of using parasite phylogeny and anchor morphometry to test hypotheses of host genealogy and ecology (Nieberding and Olivieri, 2007) in grey mullets Manuscript to be reviewed (Teleostei: Mugilidae), a speciose fish family that is economically important (Durand et al., 2012). Moreover, analysis of patterns of congruence between the phylogenies of the fish host species and their Ligophorus parasites can provide insights into prevalence of host switching (Zietara and Lumme, 2002;Vanhove and Huyse, 2015) and thence the relative importance of allopatric and sympatric speciation (Huyse et al., 2003) in shaping the diversity of this genus. It is also possible to expand this analysis in a biogeographic context by sampling different geographical populations of a host species, since some Ligophorus species have been reported to be useful biological markers of geographical fish host populations (El Hafidi et al., 2013). Lastly, we suggest that the present methods of generating and treating geometric morphometric data from the anchors, and possibly other sclerotized hard parts, can be extended to other monogenean families with little difficulty.

CONCLUSIONS
The presence of significant phylogenetic signal in the anchors makes the quantitative analysis of their shape and size variables useful in answering species discrimination and evolutionary problems in Ligophorus specifically, and in other monogenean genera more generally. In this study, we inferred two major host-specific clades from DNA sequence data, which corroborated well with clades inferred from geometric morphometric data of anchors. We further extracted size information through the first principal component of size variables based on all pairwise Euclidean distances between landmarks, and showed that the Ligophorus species infecting Moolgarda buchanani generally evolved larger anchors compared to those infecting Liza subviridis. Anchor shape was correlated with anchor size after controlling for the effect of phylogeny. Subsequently, through the analysis of directional change, we discovered two new morphological characters based on the length between the inner root length and the outer root length, and the length between the inner root point and the dent point, which proved more Manuscript to be reviewed phylogenetically informative than existing characters based on the inner length and the outer length. The Rohde-Hobbs hypothesis was validated in the 13 species of Ligophorus considered, suggesting the exploitation of different microhabitats and subsequent evolution of reproductive barriers in the form of highly variable male copulatory organ size and morphology after intrahost speciation. Finally, we demonstrated evidence for significant interspecific morphological integration of the root and point compartments within the anchors, as well as integration of the same compartment between the ventral and dorsal anchors.

SOFTWARE AVAILABILITY
The monogeneaGM package is available for download at https://cran.rproject.org/web/packages/monogeneaGM/. Analyses in this study can be replicated using the R scripts deposited in Supplemental Material C1.

ACKNOWLEDGMENTS
This work is dedicated to the memory of Lee Hong Susan Lim (1952Lim ( -2014see Gibson and Ng, 2014), who conceived the present project. Susan Lim contributed actively to global monogenean systematics for decades, and was instrumental at developing and transmitting this art in Malaysia. Klaus Rohde read the early version of the manuscript and provided helpful feedback. We thank Thian Liang Cheow for contributing R codes for data processing. Maarten Vanhove helpfully pointed out some very recent studies that were useful for the Discussion section. Finally, we are grateful to Jean-Lou Justine (Academic Editor), Timothée Poisot, and an anonymous reviewer for their constructive comments which stimulated important improvements to the present work.