Assessing population-level morphometric variation of the Mountain Mullet Agonostomus monticola ( Teleostei : Mugilidae ) across its Middle American distribution

Population-level morphometric variation of the Mountain Mullet (Agonostomus monticola) was assessed in 419 adult specimens from 25 sample sites (river basins) across its Middle American distribution (Pacific and Atlantic-Caribbean drainages). This analysis was based on 36 standardized linear measurements and 19 landmarks on geometric morphometrics approach. Discriminant function analysis (DFA) revealed 19 linear morphological characters with significant variation among groups. Geometrically, the most notable changes were associated to the curvature of the frontal region of the head, the anterior and posterior insertion of the first dorsal and anal fins. The resulting grouping based on the DFA and geometric morphometrics techniques (Pacific-A, Pacific-B and NE México-Caribbean) were similar to those previously recovered by genetic techniques, where the Pacific-B (Ayuquila river basin) was the most different group. Our results provide morphological evidence for considering Agonostomus monticola as a complex of evolutionary entities, represented by two forms in the Pacific Ocean and another in the Atlantic Ocean.


Introduction
The family Mugilidae, commonly known as "mullets," includes 75 currently recognized species (Nelson, 2006;Froese, Pauly, 2016).Most are largely marine and inhabit mainly estuarine-lagoon environments of tropical and subtropical regions of the world (Thomson, 1997;Castro-Aguirre et al., 1999); however, a few species in the family are diadromous and move between freshwater and marine waters (Harrison, Howes, 1991).Taxonomic nomenclature among members of this family has been a common theme of discussion (cf.Ghasemzadeh, González-Castro, 2016), where at least 280 nominal species have been described of which only 75 are currently recognized as valid (Harrison et al., 2007;Durand et al., 2012).The few revisionary studies of the family Mugilidae (e.g., Harrison, Howes, 1991; Morphological variation of the Agonostomus monticola Neotropical Ichthyology, 15(4): e170036, 2017 2 e170036 [2] Thomson, 1997) have fallen into frequent contradictions on the taxonomy of its members, in part because of the conservative morphology (Aurelle et al., 2008) and limited number of specimens examined along the distributional range of each species.However, the recent application of new techniques and procedures (molecular taxonomy and geometric morphometrics) in the study of the family Mugilidae have provided a different perspective on phylogeny of the taxa in this group, including the presence of cryptic species (Durand, Borsa, 2015;Ghasemzadeh, González-Castro, 2016).
Recently, molecular studies of A. monticola have indicated that the species may be more diverse than currently recognized.In this sense, Durand et al. (2012) conducted a family-wide phylogenetic study of Mugilidae using multiple mitochondrial markers.They recovered a non-monophyletic Agonostomus, and demonstrated the existence of multiple lineages within the widespread A. monticola.Later, McMahan et al. (2013) conducted a phylogeographic study of A. monticola based also on molecular data (mitochondrial and nuclear) and recovered four lineages, with origins from the late to middle Miocene (14.7 to 7.0 million years).
Despite the fact that A. monticola has been traditionally considered as one morphotype in America (Nelson, 2006), the widespread distribution coupled with the results from molecular phylogenetic studies suggests that this species is in need of morphological study, as multiple lineages may be present.No previous study for evaluating morphological variation throughout the distributional range of this species is known (Castro-Aguirre et al., 1999); therefore, this is the first assessment to characterize and quantify the morphological variation of A. monticola in its Middle American distribution.Our main question is to determine whether the morphometric variation of A. monticola through its Middle American distribution is corresponding with the lineages that have been determined by previous genetic studies.

Material and Methods
Voucher specimens of A. monticola examined in this study were came from 25 localities throughout Middle America including coastal states from the Mexican Pacific (Sinaloa to Chiapas, including Baja California Sur; referred to as Pacific-A), from the Atlantic of northeastern México (Tamaulipas and San Luis Potosí; referred to as NE México) and from the Caribbean region of Honduras to Panama including Jamaica (Fig. 1).Additionally, we referred specimens from Ayuquila drainages in Jalisco, México as Pacific-B following McMahan et al. (2013)  Thirty-six morphometric characters (linear measurements) were analyzed for comparative analysis: Lca (caudal longitude), Do (eye diameter), Am (minimum depth), Amx (maximum deptht), Lv (pelvic fin length), Ld1  All measurements were made on the left side of specimens in millimeters (mm), taken by digital calipers (precision 0.01 mm) connected to a computer.Nonlinear regression was used to remove the size component from shape measurements (allometry) and to homogenize variance (Elliott et al., 1995) Quantification of body shape was assessed using landmark-based geometric morphometrics (Rohlf, 2005).A Nikon digital camera (16 megapixels) was used to record images from the left side of each specimen (N= 419), with a ruler with 1-cm gradations.The x and y coordinates of 19 landmarks, of which three were semi-landmarks (L3, L4 and L17; Fig. 2) were used to capture variation in body shape.All landmarks selected follow the conventional rules of homology, reliability and the ability to describe the geometry of the form of interest (Zeldicth et al., 2012).In order to ensure accurate placement of the landmarks on the specimen, pins were placed on these points, which were digitized using tpsDIG v2.10 (Rohlf, 2005).Using MORPHOJ v1.05f (Klingenberg, 2011), the general Procrustes superimposition procedure was conducted to Morphological variation of the Agonostomus monticola Neotropical Ichthyology, 15(4): e170036, 2017 4 e170036 [4] remove all non-shape variation (position, orientation, and size biases) and projects the data to the tangent space by orthogonal projection.In order to assess allometric shape variation, a regression between shape data and log centroid size was performed.The total amount of shape variation was expressed as a percentage of the total variation.Centroid size was computed as the square root of the sum of squares of the distances from all landmarks to their centroid (Zeldicth et al., 2012).Covariance matrices are created from datasets of shape data after Procrustes superimposition and residuals from the allometric regression.Finally, a canonical variate analysis (CVA) was performed in MORPHOJ v1.05f to discover shape features that best distinguish between groups.All individuals were analyzed together given the absence of secondary sexual dimorphism.

Results
The mean (M), standard deviation (SD) and coefficient of variation (CV) for the 36 linear body measurements (n = 419) are shown in Tab. 1. Results showed that body measurements with the highest CV were pelvic fin base (L16-19, 22.8%), posterior insertion of the pelvic fin to anal fin origin (L19-15, 19.3%), eye diameter (Do, 17.6%) and distance from the superior origin of the operculum to insertion of the pectoral fin (L5-6, 17.5%).The linear measurements with the lowest CV were posterior insertion of the anal fin to origin of the second dorsal fin (L14-10, 5.3%), origin of the anal fin to origin of the second dorsal fin (L15-10, 5.5%) and distance from the snout tip to pelvic fin insertion (L1-16, 6.0%).
Tab. 1. Mean (M), standard deviation (S.D.) and coefficient of variation (C.V.) of standardized morphometric characters for Agonostomus monticola, based on Hubbs, Lagler (1958) and Bookstein et al. (1985)   The DFA returned a global Wilk´s lambda value of λ= 0.11 (p < 0.0), with an overall correlation value (r 2 ) of 0.867, which indicated robustness of the model.Nineteen of 35 variables (linear measurements) were statistically significant, of which the highest values in the first canonical root (CR1) were associated with maximum depth (Amx, -0.885), minimum depth (Am, -0.578), and anterior insertion of the first dorsal fin and pelvic fin origin (L8-16, -0.574).The total percentage of variation explained in CR1 was 73%, while the second canonical root (CR2) accounted for 21% of the total variation (Tab.2).Based on squared Mahalanobis distances (D 2 ), the general accuracy of classification of individuals was 91%, where the highest classifications were obtained for individuals from the Caribbean (95.36%;Tab. 3, Fig. 3).Fig. 4 depicts the phenogram derived from a cluster analysis using the criterion of Chebyshev distance expressed as a percentage.This showed that individuals of A. monticola from Caribbean and NE México drainages share a considerable extent of homogeneity.Nevertheless, specimens from the Pacific-A group were more similar to those from the Atlantic (Caribbean and NE México), while the Pacific-B population had a larger distance to other groups.

Tab. 3. Mahalanobis squared distances between
Agonostomus monticola groups using: (a) linear measurements (Hubbs, Lagler, 1958;Bookstein, 1982) and (b) landmark-based method (Rohlf (2005)).Geometric morphometric analysis indicated 7% of global shape change was due to allometric variation, which was assessed by means of the regression between shape data and log centroid size.A deformation grid (Fig. 5) exhibits in the mean shape of A. monticola, where the most notable changes were shifts of semi-landmark 3 (curvature of the frontal region of the head), landmark 8 (anterior insertion of first dorsal fin), landmark 9 (posterior insertion of first dorsal fin), landmark 14 (posterior insertion of the anal fin) and landmark 15 (anterior insertion of the anal fin).In the canonical variate analysis CVA, the first canonical variate (cv1) explained 73.5% of the variation, while the second canonical (cv2) explained 21% (Fig. 6).A p-value obtained from permutation tests (100000 permutation rounds) for Mahalanobis distances among groups was significant (< 0.001; Tab. 3).Fig. 6 shows multivariate relationships among groups, similar to that obtained by the DFA using linear body measurements.

Discussion
The present study focused on morphological variation of Agonostomus monticola throughout its Middle American distribution shows significant differences for at least 19 linear body measurements, four landmarks and one semilandmark.Characters with the highest discriminate values are associated with distances in the body depth at level of the region comprised between the first and second dorsal fin and the anal fin.Likewise, the ranges of variation for some linear morphological characters examined in this study are similar to those previously reported (cf.Thomson, 1997;Camacho-Rodríguez, Lozano-Vilano, 2004), who used non-standardized measurements.Some of the linear measurements of this analysis are consistent with those used in morphological studies for some taxa of the family Mugilidae and have found to be effective for taxonomical discrimination (e.g., Ibáñez-Aguirre, Lleonart, 1996;Cousseau et al., 2005;Ibáñez-Aguirre et al., 2006;González-Castro et al., 2008, 2012;Konan et al., 2014).These last authors determined for the genus Mugil that the landmarks associated to the head region had a higher discriminative power.In our study, the most discriminate variables were associated with the body depth region.
In the case of the geometric morphometry method, authors as Corti, Crosetti (1996) and González-Castro et al. (2012) have successfully used this procedure for the morphological discrimination of species of the genus Mugil.The first authors determined a high variation explained by the first canonical variable, with the highest contribution of the landmarks associated to the first dorsal fin; while the second authors found a higher discrimination in the landmarks of head length, caudal peduncle depth, body depth at the origin of the second dorsal fin, pectoral fin length and anal fin length.
A similar result was obtained in our analysis, except that these authors did not include the use of semi landmarks as we did.Therefore, we provide a set of morphological features (e.g., curvature of the frontal region of the head, the anterior and posterior insertion of the first dorsal and anal fins) that might be useful for the discrimination of forms of this species using multivariate statistical techniques.
Despite the apparent morphological similarity of A. monticola throughout its distributional range (cf.Schultz, 1946;Harrison, Howes, 1991;Thomson, 1997), we are able to distinguish the presence of different shapes in this taxon due to the exploration of new characters related to body morphology.In this sense, Miller et al. (2005) pointed out the potential presence of more than one group within this taxon (referring to a clade) in Mexican inland waters.Likewise, McMahan et al. (2013) using nDNA and mtDNA sequences recovered at least four clades for this species through its Middle American distribution (two in the Atlantic and two in the Pacific), which were considered for the classification of groups (lineages) in the present analyses.These authors mentioned the possible existence of a complex of cryptic taxa within this species but no morphological differences were reported.However, we found trends of morphological differentiation between some of these groups, especially the Pacific-B (Ayuquila drainage) and those of the Caribbean and NE México.Thus, 7 e170036 [7] we found consistency in the morphometric relationships that are corresponding to geographic areas, supporting the hypothesis of lineages successfully separated.
Corroborating prior work (e.g., Harrison, Howes, 1991;Thomson, 1997), we find a high similarity between the populations of the Caribbean and NE México; however, molecular evidence supports the existence of two different groups within this region (Durand et al., 2012;McMahan et al., 2013).Otherwise, our analysis shows that individuals from the Atlantic slope (NE México and Caribbean) are morphologically similar.
NE México and Caribbean populations are clustered in analyses and both are morphologically most similar to Pacific-A while Pacific-B (Ayuquila drainage) branches as more divergent from these three groups.The first branch of the phenogram contains specimens of A. monticola from Pacific-A, with this group encompassing the wider ranges of distribution that could be explained by dispersion events of individuals coming from the same reproductive stock in a pelagic marine environment.It has been suggested that adult A. monticola migrate from mountain streams to the open ocean for spawning (Anderson, 1957;Cruz, 1987;Phillip, 1993) because early stages have been found closer to the mouth of the rivers.In this pelagic environment, early life stages are dispersed by surface currents promoting population mixing (Franco-Gordo et al., 2002) as suggested for other mugilids in the Mexican Pacific (Ibáñez et al., 2012).
In summary, we argue a possible scenario for the geographical disjunction of the two forms of A. monticola in the Pacific (Pacific-A and Pacific-B).The morphological and phylogenetic differences detected between them could be explained by genetic isolation in this region, as reported by McMahan et al. (2013).It possible that genetic isolation of Pacific-B lineage (Ayuquila drainage) is related to the fact of reproduction of adults and development of juveniles in freshwater fully conditions without migration to coastal waters.Finally, we highlight the need of studies of comparative osteology and population genetics for these lineages reported here, as well as to focus assessments on the life history of the Ayuquila population (Pacific-B), in order to robustly clarify the taxonomy and evolutionary relationships of lineages within this species.
, computed for each linear character as Ms = Mo (Ls/Lo) b ; where Ms= standardized measurement (mm), Mo = measured character length (mm), Ls = mean (arithmetic) standard length (mm) for all individuals from any group examined, Lo = standard length (mm) of the specimen; parameter "b" was calculated by the logarithmic (base 10) regression equation, log Mo= log a + b log Lo, using all specimens in any group.Descriptive statistics (mean and standard deviation) of all standardized linear measurements were recorded for each group.The coefficient of variation (CV) was computed for each standardized linear character according to the equation: CV = 100 × SD/X, where SD is Standard Deviation and X is the mean of the standardized linear measurements of morphometric characters in each group.All standardized linear measurements were compared among groups by means of discriminant function analysis (DFA).Compared groups included Pacific-A (all Mexican populations in the Pacific slope from Sinaloa and Baja California Sur to southern Chiapas, with the exception of the Ayuquila drainage in Jalisco, n = 202), Pacific-B (Ayuquila drainage, n = 27), Northeastern México (Tamaulipas and San Luis Potosí, n = 17), and Caribbean (including Central America countries and Greater Antilles, n = 180).Statistical significance of discrimination among groups was determined using Wilk's lambda (λ), which ranges from 0.0 (perfect discrimination power) to 1.0 (absence of discrimination power).Standardized coefficients of canonical variables were determined by estimating the contribution of each variable to each canonical function.Finally, a phenogram was built based on squared Mahalanobis distances.All statistical analyses were performed using Statistica 6.0 (StatSoft Inc Tulsa OK, 2002).

Fig. 3 .Fig. 4 .
Fig. 3. Plot showing the classification of individualsAgonostomus monticola resulting from the discriminant function analysis using linear body measurements.

Fig. 5 .
Fig. 5.A deformation grid shows mean shapes of all specimens using the residual data for Agonostomus monticola.

Fig. 6 .
Fig. 6.Two-dimensional ordination based on Canonical variate analysis (CVA) of Agonostomus monticola specimens from the four groups.
. All specimens examined in this study are from fish collections of the Centro Universitario de Ciencias Biológicas y Agropecuarias (CUCBA-UDG) in Guadalajara, Jalisco, México; The Field Museum of Natural History (FMNH), in Chicago, Illinois, USA; Louisiana State University Museum of Natural Science (LSUMZ) in Baton Rouge, Louisiana, USA; Southeastern Louisiana University (SLU) in Hammond, Louisiana, USA; Universidad Autónoma de Baja California (UABC) in Ensenada, Baja California, México; Universidad Autónoma de Nuevo León (UANL) in Monterrey, Nuevo León, México; University of Southern Mississippi (USM) in Hattiesburg, Mississippi, USA.A list of fish material examined appears in Appendix (S1 -Available only as online supplementary file accessed with the online version of the article at http://www.scielo.br/ni).
Tab. 2. Standardized coefficients of canonical variables (cv) and values of lambda of Wilks, significance (p-level), cumulative probability (cum.pro.) and tolerance (toler.)resulting from the discriminate function analysis of the morphological analysis of Agonostomus monticola (see Tab. 1 for abbreviation of variables).Highest values in bold.