The Mayan Cichlid (Mayaheros urophthalmus) is usually considered to be a complex of 18 subspecies, most of which are endemic to the Yucatan Peninsula and were diagnosed by phenetic analyses based on traditional morphometrics and color pattern. However, morphological differences can be due to environmental conditions rather than taxonomic distinctiveness. We evaluated, by means of a geometric morphometric analysis, two hypotheses for shape differences in 20 natural populations of M. urophthalmus, including five subspecies recently raised to species status: M. alborus, M. cienagae, M. conchitae, M. mayorum, and M. zebra. The geographical distribution and three types of aquatic environment (River, Lagoon or Pond, and Cenote) were used as classificatory variables. In addition, a molecular analysis of two concatenated fragments of mitochondrial DNA (mtDNA), cytochrome b (cytb) and cytochrome oxidase subunit I (COI), showed genetic differentiation among some populations (FST = 0.36). Whereas the geometric morphometric analysis found significant differences among all aquatic environments, patterns of body shape of M. urophthalmus are more consistent with ecophenotypic variation than with genetic differentiation due to geographic isolation by distance. We think that there is currently no evidence to raise the traditionally recognized subspecies of M. urophthalmus to the species level.
THE Mayan Cichlid (Mayaheros urophthalmus) is a freshwater fish that ranges from southern Veracruz, Mexico to Nicaragua in Central America, inhabiting rivers, lakes, ponds, and cenotes (sinkholes) in the Yucatan Peninsula (YP). It occurs as an exotic species also in Oaxaca (Mexico), Florida (USA), and the Philippines (Loftus and Kushlan, 1987; Espinosa-Pérez et al., 1993; Ordóñez et al., 2015).
Based on morphological characters and color pattern, Hubbs (1935, 1936, 1938) described 12 subspecies of Cichlasoma (=Mayaheros) urophthalmus, seven of them endemic to YP: C. u. aguadae, C. u. amarum, C. u. cienagae, C. u. conchitae, C. u. mayorum, C. u. trispilum, and C. u. zebra. Barrientos-Medina (1999, 2005) performed phenetic analyses of the “C. urophthalmus complex” in the YP based on traditional morphometrics, meristics, and coloration. He found statistically significant differences between putative species (Hubbsian subspecies and undescribed forms); however, all the proposed new species, or subspecies raised to species level, overlap in these quantatitive traits.
For this study, some of the former subspecies of Hubbs (1936) were reanalyzed: M. alborus was diagnosed chiefly by Barrientos-Medina (2005) as having a body depth 44% SL (standard length), a suborbital bone wider than 22% HL (head length), 31–32 scales on the lateral line, and 7 vertical bars on sides (wider than the interspaces); M. cienagae is 43–50% SL in body depth, suborbital width lower than 20% HL, anal-fin base 24–27% SL, 28–31 scales on the lateral line, and 7 vertical bars on sides (narrower than the interspaces); M. conchitae is also 44% SL in body depth, has 25 scales on the lateral line, and 15 or more dorsal spines; M. mayorum has a moderately deep body (40–45% SL), suborbital width 14–19% HL, 9–12 predorsal scales, 29–31 scales on the lateral line, bars as M. cienagae, and the caudal ocellus joined to the last bar; and M. zebra also has a body depth 40–46% SL, suborbital width 13–24% HL, 12 or fewer predorsal scales, 28–32 on the lateral line, bars as M. cienagae, and the caudal ocellus not joined to the last bar. All of these traits apply also to other putative species (e.g., the unnamed forms from Chetumal Bay and Laguna de Términos). A possible exception is the low dorsal-fin spine count (14 spines) of M. aguadae, which is known only from the holotype (no topotypes known) and is not analyzed here.
In spite of the lack of diagnoses beyond significantly different but overlapping continuous characters, Barrientos-Medina's (2005) conclusion was that up to 18 species should be recognized, including raising all subspecies to species status. The proposal was followed by Kullander (2003) and more recently by Říčan et al. (2016), who described the genus Mayaheros for the “urophthalmus complex” and claimed that “[t]he range of M. urophthalmus sensu stricto is limited to the Lake Petén Itzá and contiguous zones, being substituted in YP, in additional parts of México, Guatemala, and Honduras by morphologically similar species, endemic and restricted in their distribution.”
However, Alfaro-Bates (1989) had previously compared two subspecies, C. u. mayorum and C. u. zebra, and she concluded that morphological differences between them were due to the ecological conditions of water bodies (i.e., the different morphs were ecophenotypes). Indeed, several of the local variants that Hubbs (1936) would later formally describe as subspecies were examined previously by Evermann and Goldsborough (1902:158), who remarked: “…Although it is easy to pick out the individual fish from any one of these localities, we do not find any structural difference of value. The color differences are due simply to the character of the water in which they were found. The difference in form is simply a question of food supply…”.
The taxonomic value of morphological and coloration characters must be considered with caution, especially in the absence of a cladistic framework, because these traits can be environmentally induced: “…all too often […] non-genetic factors have led to the descriptions of numerous species and subspecies of no taxonomic value” (Quicke, 2013:33). In addition, Razo-Mendivil et al. (2013) found low levels of genetic divergence 0.4% (Dxy = 0.004%) among 26 populations of Mayan Cichlid, based on cytochrome b sequences (599 bp). There was a weak association with geographical isolation by distance.
In this paper, we revisit the morphological and genetic diversity in 15 natural populations of M. urophthalmus and five nominal species (former Hubbsian subspecies: M. alborus, M. cienagae, M. conchitae, M. mayorum, and M. zebra, although only the first two were included in our genetic analysis, see Appendix 1 and Fig. 1), with emphasis on YP, where most of the putative species are distributed.
Fig. 1.
Study area, Yucatan Peninsula. Localities of M. urophthalmus taken into account for the geometric morphometric analysis, in gray numbers. Localities with genetic data, marked with an asterisk. Localities from nominal species included in the morphometric analysis: M. alborus (15), M. cienagae (10), M. conchitae (8), M. mayorum (6), and M. zebra (9).
The shape was analyzed from the perspective of the landmark-based geometric morphometrics method (Adams et al., 2004). Geometric morphometrics considers a constellation of discrete anatomical loci, each one described by 2- or 3-dimensional Cartesian coordinates. The spatial relationship between landmarks supplies a graphical visualization of differences and allows quantitative descriptions and comparisons.
We explored two possible causal agents of morphometric differences in the studied populations of M. urophthalmus. First, we treated the sampled populations as independent units capable of reflecting morphological differences due to a process of genetic differentiation that resulted from isolation by distance. We wanted to explore whether isolation by distance was correlated with the observed geometric morphometric variation. Consequently, the geographic demarcation of populations (i.e., locality) was treated as a classificatory variable for the analysis. Second, we classified the sampled populations into three types of aquatic environment (1 River, 2 Lagoon or pond, 3 Cenote) of M. urophthalmus. This way we intended to infer whether shape differences could be interpreted as ecophenotypes.
In addition to the shape analysis, we implemented a molecular approach using two concatenated fragments of mtDNA, cytb, and COI, both widely used in phylogeography and in delimitation of species, because they contain both conserved and variable regions (Ward et al., 2009). Our expectation is that a pattern consistent with accepting previous species hypotheses would imply geographic distance (isolation) rather than environment as a better explanation of morphological differences. We acknowledge that adaptation to different habitats can also be a motor of speciation, but then we would expect that morphological and genetic disjunctions should strongly correlate as well.
MATERIALS AND METHODS
Collection sites.—We surveyed 15 sites (Fig. 1) in the YP (Mexico) and Belize, with the aim of incorporating the spatial genetic variation of M. urophthalmus on the YP, because most of the species recognized by Barrientos-Medina (1999, 2005), i.e., the subspecies of Hubbs (1935, 1936, 1938), are distributed on YP. Tissue samples (tips of dorsal fins) were obtained from fish caught using cast nets and traps and stored in 90% ethanol for transport to the laboratory. Voucher specimens were deposited in the fish collection at El Colegio de la Frontera Sur, Chetumal, Mexico (acronym ECO-CH). Although only seven localities were used for both the morphometric and the molecular analyses, both analyses include material from the Caribbean and Gulf of Mexico versants, as well as north and south, of the YP (Fig. 1). However, for the genetic analysis only two out of the five former Hubbsian subspecies were included, M. cienagae from Progreso (Yucatán) and M. alborus from Río Zapata (Tabasco), because some type localities have now disappeared, for example, the Conchita and Holpuch cenotes, that were filled with debris and paved for the construction of a road and a parking lot, respectively, and others are in protected archaeological sites as Cenote Xtolok at Chichén Itzá and Cenote Xlakah at Dzibilchaltún, from where fresh samples were not available.
Morphometric methods.—To evaluate the body shape variation of populations of M. urophthalmus, 152 adult specimens between 90 and 160 mm SL were examined from three ichthyological collections (see Appendix 1). The specimens were photographed and digitized on the left side of the body, under the same light intensity and adjusting the lens in the smaller fish to recognize landmarks with precision. We used a digital camera (SONY Cyber-shot, 10.1 megapixels) mounted on a tripod; all pictures included a ruler as a reference.
We set 11 landmarks and 13 semilandmarks (non-homologous points used for recover morphometric data along curves) for the supracephalic profile, for a total of 24 reference points (Fig. 2), using TpsDig 2.17 (Rohlf, 2013). To superimpose landmarks and semilandmarks, and to remove size as a variable via Procrustes superimposition in tpsRelw, first we built a slider file in tpsUtil, and then we selected a chord-min d2 as a slide method that uses the distance-minimization from the consensus to specimen approach or Procrustes distance. After that we saved the aligned data and centroid sizes (CS) and checked for digitizing errors; we conducted a visual inspection of relative warps as recommended by Zelditch et al. (2004) via a principal components analysis.
Fig. 2.
Constellation of landmarks and semilandmarks for the geometric morphometric analysis showing Procrustes superimposition of landmarks and semilandmarks. Anatomical positions of landmarks: Landmark 1—intersection between the posterior base of dorsal fin and caudal peduncle; Landmark 2—the most distal pore of the end of lateral line; Landmark 3—intersection between the posterior base of the pelvic fin and caudal peduncle; Landmark 4—intersection between anterior base of pelvic fin and belly; Landmark 5—apical intersection between subopercle and interopercle; Landmark 6—commissure of dentary at the more posterior end at its junction with the anguloarticular, quadrates, and ectometapterygoid bone; Landmarks 7 and 8—delimitation of the width of the commissure of dentary at its junction with the maxilla; Landmarks 9 and 10—delimitation of the width of the dentary at the more anterior end; Landmarks 10 and 12—delimitation of the width of the premaxilla at the more anterior end; Landmark 11—dorsal end of edge between opercle and cheek; Landmarks 12 to 24—the curve of 13 semilandmarks beginning at the intersection between the premaxilla and ascending process at semilandmark 12, the curve traced upwards until the base of dorsal fin at semilandmark 24.
To test for differences in mean CS among all populations of M. urophthalmus, we performed a one-way NPMANOVA based on pairwise Euclidean distances between individuals, with 9999 permutations, using PAST 2.17c (Hammer et al., 2001). Posteriorly, as a post-hoc test, pairwise comparisons between groups were based on Hotelling's tests with a sequential Bonferroni correction.
Afterward, a canonical variates analysis (CVA) of shape was performed, based on two-dimensional vectors (the partial-warps scores), with the software IMP8 (Sheets, 2003). A visual presentation of shape differences described by the CVA was produced by regressing the shapes on the first two canonical vectors (Zelditch et al., 2004); this permitted the splines of the shape change to be associated with positive and negative values of canonical vectors. The plots of the deformation implied by the CVA axes were visualized as thin-plate spline with IMP8 (Sheets, 2003). The regression shows all the changes in shape correlated with the CVA axis score to a factor of deformation of 10% on each analysis, but the CVA axis itself does not show all the correlated change in shape.
Finally, we performed two Mantel's tests using PAST 2.17c (Hammer et al., 2001) to infer the possible correlation of body shape variation with geographic distances and with genetic distances between populations. For this, we employed the Euclidean distances between CS, and Slatkin's linearized FST, with 5000 permutations, at P < 0.05.
Molecular methods.—The present concatenated dataset is composed of two mtDNA fragments, the protein-coding cytb gene (1061 bp) and a portion (616 bp) of the protein-coding gene COI. The concatenated molecular matrix includes 1677 aligned positions from 81 individuals from 15 natural populations of M. urophthalmus (Fig. 1). Approximately 2 mm3 of tissue was ground in DNeasy Blood & Tissue Kit (QIAGEN) for total genomic DNA extraction following the manufacturer's protocol. Farias et al. (2001) and Chakrabarty (2006) describe oligonucleotides used for cytb and COI DNA amplification, respectively. The process was conducted in Peltier-effect thermocyclers (MultiGene OptiMax Thermal Cycler) using the following parameters: one initial cycle at 95°C for 120 s, followed by 35 cycles of 95°C for 20 s, 50°C for 20 s, 72°C for 60 s, with one final cycle at 72°C for 240 s. All PCR reactions were conducted along with positive and negative controls to detect potential false positives due to contamination. Products of PCR were visualized on 2% agarose gels stained with ethidium bromide. Successful amplifications were purified using Wizard® SV Gel and PCR Clean-Up System (Promega). Purified products of PCR were sent for sequencing in both directions to Macrogen (South Korea). Sequence files were analyzed with the aid of the program BioEdit Sequence Alignment v7.0.9 (Hall, 1999). All sequences were deposited in GenBank under the following inclusive accession numbers: MF741939 to MF741974, and MF776666 to MF776701.
Genetic variation within and among populations was evaluated using empirical descriptive values such as haplotype diversity (h), nucleotide diversity (π), and segregated sites (Nei and Kumar, 2000). Population stasis was evaluated using Fu's F (Fu, 1997) and Tajima's D (Tajima, 1989) tests. In addition, a non-hierarchical analysis of molecular variance (AMOVA) was performed. The fraction of the genetic variation distributed among populations was estimated using the FST statistic (Lynch and Crease, 1990). For its computation, a hierarchical analysis of variance among populations was used (Weir and Cockerham, 1984). To obtain confidence intervals for the FST estimates, one million non-parametric permutations of haplotypes were performed. Alternatively, a ΦST (Excoffier et al., 1992) test under a Tamura-Nei correction was calculated.
Finally, for inferring the possible pattern of long distance isolation, a Mantel's test under the Slatkin's linearized FST matrix was performed and Euclidean distances were used (5000 permutations, P < 0.05). The calculations of all the above descriptors were performed with the programs DNASP 6 (Rozas et al., 2017), Arlequin 3.5.2.2 (Excoffier and Lischer, 2010), and PAST 2.17c (Hammer et al., 2001).
To choose the model of molecular evolution that best fitted our concatenated sequences data under the Akaike information criterion (AIC, the model selected = TIM3+I p-inv = 0.9620), we used jModelTest 2.1.10 (Darriba et al., 2012). To infer evolutionary relationships among haplotypes, a Bayesian analysis was performed using MrBayes 3.2 (Ronquist et al., 2012), with the following parameters: Lset Nst = 6, revmatpr = (0.3049, 7.4737, 1.0000, 0.3049, 2.5344, 1.0000), statefreqpr = (0.2395, 0.3215, 0.1525, 0.2866), and pinvarpr = 0.9620, as suggested by jModeltest. The search ran for 15 million generations, with four parallel chains and sampling every 1000th generation. The burn-in value was set to 25%. A majority-rule consensus tree after burn-in was used to estimate the posterior probabilities for each node. To polarize the tree, we employed the species Petenia splendida and Darienheros calobrensis as outgroups (we took existing sequences from GenBank with the following accession numbers: AF370679, EU751899, AY843381, GU817255). The ingroup was not constrained as monophyletic for the analysis.
As an alternative to the Bayesian trees, a phylogenetic network for the haplotypes was inferred by means of NETWORK 5.0.0 (Bandelt et al., 1999), a median-joining (MJ) algorithm based on genetic distances was calculated with an epsilon value = 10, and as a postprocessing to purge superfluous links and median vectors from the network, a Maximum Parsimony calculation (MP) was executed (Polzin and Daneshmand, 2003).
RESULTS
Geometric morphometrics, geography, and environment.—The CS of body shape variation of all populations of M. urophthalmus showed a normal distribution (Monte Carlo test, P = 0.0091; P-P plot, PP = 0.978). Using sampled sites as a classificatory variable of shape, there were significant differences of CS means among some, but not all, populations (F19,3.4e07 = 33.56, P = 0.0001; Table 1). The nominal species M. alborus, M. cienagae, and M. zebra overlapped with the others in the CVA scatterplot (Fig. 3). Mayaheros conchitae from the destroyed cenotes Conchita and Huolpoch, and M. mayorum from the cenote Xtolok at archaeological zone Chichén Itzá, were peripheral in the CVA scatterplot, as well as individuals from River Palizada, from Payo Obispo in Chetumal, and from Pucteito; M. conchitae did not show any difference.
Table 1.
Non-parametric PermANOVA on Euclidean distances of centroid size for independent analysis of geographic location and aquatic environment as classificatory variables (9999 permutations).
Fig. 3.
CVA 1 vs. CVA 2 by population at geometric morphometric space delimited by minimum convex polygons, and the grid of morphometric space of CVA with a factor of deformation = 0.1; (A) axis 1 of CVA and (B) axis 2 of CVA.
In contrast, when sampled populations were classified according to aquatic environment (River, Lagoon-pond, Cenote), there were significant differences in CS means among all habitats (F2,3.4e07 = 12.52, P = 0.0001; Tables 1, 2). All populations, including the nominal species (M. alborus, M. cienagae, M. conchitae, M. mayorum, and M. zebra), overlapped within their corresponding type of environment in the CVA scatterplot (Fig. 4). The widest overlap occurred between lentic environments (i.e., Lagoon-pond and Cenote).
Table 2.
Pairwise a-posteriori tests among aquatic environments, P values of sequential Bonferroni correction significances above the diagonal and F values below.
Fig. 4.
Axis of CVA under aquatic environment approach with ellipses of the probability of 95%, and the grid of morphometric space of CVA with a factor of deformation = 0.1; (A) axis 1 of CVA and (B) axis 2 of CVA.
The thin-plate spline under the assumption of sampled sites as independent populations (Fig. 3) displayed a slight dorsoventral contraction of the body, plus an expansion of the snout in landmarks 7, 8, 9, and 10 on axis 1 of the CVA (which explained 36.74% of the variability). In addition, landmarks 18 to 24 showed a dorsoventral contraction and a horizontal expansion of the grid at the head. In axis 2 of the CVA (15.49% of the variability), the major grid deformation is on landmark 24, relative to the insertion of the dorsal fin, and on landmarks 5, 6, 7, and 8, relative to the snout, dentary, and maxilla. However, to reach 95% of total variation, nine axes of the CVA were needed. Mantel's tests on Euclidean distances between CS versus geographic distances and versus genetic distances among populations were not significant (r2 = –0.025, P = 0.597; r2 = 0.48, P = 0.093, respectively), confirming that there was no significant geographic or genetic effect on body shape.
The thin-plate spline by the environment (Fig. 4) displayed a contraction of the grid at the level of dorsal fin in respect to the elongation of the caudal peduncle (landmarks 1, 2, and 3), a displacement upwards of the maxillary position of landmark 10, and an expansion downwards of landmarks 5, 6, 7, and 8. An anteroposterior-ventral axis expansion also occurred at the head. The most important deformation of body shape occurred in axis 1 of the CVA, which explained the 70.79% of the variability of the body shape. In addition, the vectors of landmarks 2 and 10 are related to placement of the end of the lower lateral line and length of the premaxilla, respectively. Axis 2 of the CVA (29.21% of the variability) showed a major deformation over the first half of the body. The vectors of landmarks 2, 3, and 4 were ventrally projected, with the greatest body depth at landmarks 17 to 24. These two axes explain 100% of the variability (Fig. 4).
Genetics and geography.—The concatenated mitochondrial sequences recovered 36 haplotypes (H) for populations taken into account in the genetic analysis, of which 28 were singletons; H12 was the most widespread haplotype, present in five of the 15 populations, whereas H6 was the most common haplotype, present in 14 individuals of three populations (Table 3). The relative nucleotide composition observed in the cytb and COI fragments was similar, characterized by sequences rich in cytosine (33.6% cytb, 30.1% COI), with intermediate percentages for thymine (28.3%, 29.3%) and adenine (24.3%, 22.4%), and poor in guanine (13.7%, 18.0%).
Table 3.
Genetic diversity within populations of M. urophthalmus; n = number of sequences, S = segregating sites, π = nucleotide diversity, h = number of haplotypes, Hd = haplotype diversity, F = Fu's F, Fp = Fu's test p value, D = Tajima's D, Dp = Tajima's Dp value, θπ = average pairwise nucleotide diversity, Ne = female effective population size. * = significant differences.
The highest haplotypic diversity was recorded at Río Palizada (1), Río Emiliano Zapata (1), and Lake Cobá (0.9), whereas the lowest diversity occurred at Laguna Manatí (0.2). Similarly, the highest nucleotide diversity (π) was recorded at Río Palizada (0.0035) and Sabancuy (0.0033), where the average number of nucleotide differences per site were highest, and the lowest π was found at Río Hondo (0.0004; Table 3).
The neutrality tests suggested a population expansion, based on large negative Fs values (Fig. 5), indicating an excess of rare mutations. At the intrapopulation level, Cobá (Fs = –1.937, P = 0.011), Río Palizada (Fs = –2.517, P = 0.017), and Sabancuy (Fs = –1.812, P = 0.02) exhibited demographic expansion. Also, Tajima's test of selective neutrality proved to be significant in the Progreso Río Canotaje population (D = –1.390, P = 0.05; Table 3).
Fig. 5.
Topology of haplotypic relationships of 15 populations of the M. urophthalmus complex inferred at a network based on two concatenated mitochondrial protein gene fragments (cytb and COI, 1677 bp) and showing the north and south components.
Molecular variance was greater within populations (68.11%, Table 4) than among populations (FST = 0.318), showing a genetic structuration. However, based on the exact differentiation test under the hypothesis of random distribution of the individuals between all pairs of populations and based on the ΦST test, significant differences occurred among some populations, namely between Cenote Azul, Laguna Manatí, Cenote Popolvuh, and the other populations (Table 5).
Table 4.
Analysis of molecular variance based on two concatenated mitochondrial gene fragments (cytb and COI, 1677 bp) of M. urophthalmus.
Table 5.
Conventional population pairwise FST values from haplotype frequencies below the diagonal and Tamura-Nei ΦST values above, both under a sequential Bonferroni correction. * Significant values.
The Bayesian tree corroborates the monophyletic relationships among the populations considered in the analysis, including two nominal species, M. cienagae and M. alborus, with high Bayesian posterior probability (100; Fig. 6). However, the Bayesian probability support values at most internodes were weak, the topology shows poor resolution, probably due to gene flow or incomplete lineage sorting among populations (haplotypes 17, 22, 23, 30 from Progreso Río Canotaje, Celestún, and Cenote San Juan del Río). Nevertheless, the Bayesian topology was highly consistent with the haplotype network inferred by statistical parsimony at median-joining results. The haplotype network displayed a high proportion of singletons (Fig. 5). However, the Mantel's test revealed no relationship between linearized (FST) values and geographic distance (r2 = –0.0028, P = 0.523).
DISCUSSION
The phenotypic differences between populations of a species may be due to genetic variation, phenotypic plasticity, or the interaction of these. The high morphological variability of M. urophthalmus has led to taxonomic decisions concerning its splitting into several subspecies or even species. However, we evaluated both morphology and molecular variation of populations of M. urophthalmus in YP, where several putative species were recognized in the absence of autapomorphies or at least consistent diagnoses by Barrientos-Medina (1999, 2005) and endorsed by Říčan et al. (2016), who did not include any of them in their phylogeny, and our results are more consistent with an ecophenotypic variation interpretation than with a genetic effect derived from the isolation by distance with posterior morphological differentiation. Elevated gene flow limits speciation and facilitates the maintenance or evolution of phenotypic plasticity (Crispo, 2008), and gene flow was moderate to high among most populations of M. urophthalmus. The exceptions were Cenote Azul, Laguna Manatí, and Cenote Popolvuh, between which gene flow was low; none of them, however, were supported as different species (data not shown). The Bayesian posterior probability values of these internodes were weak, and none of these populations were recovered as monophyletic. It is true that plasticity in small, isolated populations can be lost due to genetic erosion (e.g., Luquet et al., 2011), but M. urophthalmus presents ecophenotypic variation consistent with population expansion, and the current pattern reflects intraspecific rather than interspecific structure, with no deeper phylogeographic structure.
Aquatic environments, on the other hand, can show great spatial and temporal variation, and many fish species show extreme morphological differences between contrasting habitats because water is a dense medium that affects performance characteristics (Langerhans et al., 2003). Environmental conditions also can cause body shape changes in freshwater fish species, such as as body depth, e.g., fish from lotic habitats have a more hydrodynamic body shape and fish from lentic habitats have greater maneuverability (Santos and Araújo, 2015). Also, lotic habitats have a high concentration of dissolved oxygen, whereas lakes and swamps are usually lower, even hypoxic at times, and this environmental condition may influence phenotypic plasticity (Crispo and Chapman, 2010). In addition, the displacement of landmarks in the body shape of M. urophthalmus suggests a change in the position of the mouth with respect to habitat. Morphs correlated with the use of resources may differ in jaw length and bluntness of snout; this has been associated with the type of prey and the ability for recognizing available prey at their habitats (Smith and Skulason, 1996). Some cichlids are an example of spectacular adaptive radiation and phenotypic plasticity that concurrently has generated a lot of body shapes, coloration patterns, and trophic morphs due to resource use (Takahashi and Koblmüller, 2011).
Morphological and molecular analyses do not always agree in New World cichlids (Farias et al., 2001; Chakrabarty, 2007; Schmitter-Soto, 2007; Říčan et al., 2008). These disagreements are the consequence of an enormous ecomorphological versatility and probable convergence of traits associated with trophic ecology and habitat use. In this sense, the morphometric differences we found among populations of M. urophthalmus from different aquatic environments are coupled with low genetic differentiation (FST = 0.31), as found also by Razo-Mendivil et al. (2013), 0.4% with cytb. In addition to the poor resolution and support at most internodes of the Bayesian analysis and networks of statistical parsimony, it is safe to state that the populations included in this study as M. urophthalmus comprise a single widely distributed species, with only an incipient divergence among north and south populations, that was not observed in our morphometric results.
The subspecies concept that Hubbs (1936) was following when he formally named so many “local variants” becomes apparent when he admits: “…I feel confident that other collections from the hundreds of unsampled waters within the general range of this species [M. urophthalmus] will produce a full series of intermediate races” (Hubbs, 1936:268); hence, a subspecies, in his view, is a variety that intergrades with other such variants. The subspecies category has since fallen out of favor with ichthyologists: if a taxon is diagnosable, then it is usually recognized as a species, under some versions of the phylogenetic species concept (see Mayden, 1997).
According to Lowe and Allendorf (2010), the m (number of migrants per generation) and FST values are the approximate values associated with the connectivity between populations under the island migration model. In this context, the three populations that displayed significant differences in FST and ΦST values (Cenote Azul, Cenote Popolvuh, and Laguna Manatí) may be considered, in the sense of Lowe and Allendorf (2010), in the category of “adaptive connectivity,” which implies that there is sufficient gene flow to spread advantageous alleles among these three populations and the others; however, this expectation needs to be corroborated, especially given that no adaptive gene is identified. On the other hand, the pattern of the haplotype network may be indicative of a population that grew from a few founders or experienced a bottleneck previous to expansion (Slatkin and Hudson, 1991).
Contrary to Barrientos-Medina (2005), Razo-Mendivil et al. (2013), who included the type locality of M. urophthalmus (Lake Petén Itzá, Guatemala), did not find evidence to support recognition of any other species of Mayaheros in YP. We concur with the latter authors: at least the putative species (former Hubbsian subspecies) that we examined morphometrically and/or genetically, i.e., M. alborus, M. cienagae, M. conchitae, M. mayorum, and M. zebra, should be considered junior synonyms of M. urophthalmus, pending additional data (e.g., an autapomorphy or at least a consistent diagnosis of characters in combination). As a taxonomic corollary, we think that raising subspecies to species status is a decision that requires the same amount and quality of evidence as the description of any new taxon.
ACKNOWLEDGMENTS
The Mexican Consejo Nacional de Ciencia y Tecnología financed the research via a postdoctoral grant to the first author (290847). L. Chumba, H. Espinosa, and M. Valdez kindly lent material under their curatorial care. R. Méndez, G. I. Cruz, and E. Guillén helped in fieldwork.
LITERATURE CITED
Appendices
Appendix 1.
Material examined for geometric morphometric analysis. Non-named forms appear as “M. urophthalmus.” The locality name includes the type of environment: River, Lagoon or Pond, and Cenote. The localities of named forms are the type localities; full description in Hubbs (1936). n = sample size. Acronyms: ECO-CH, Colección Ictiológica de El Colegio de la Frontera Sur, Unidad Chetumal; CZOO-CCBA-UADY, Colección Zoológica del Campus de Ciencias Biológicas y Agropecuarias de la Universidad Autónoma de Yucatán; CNP IBUNAM, Colección Nacional de Peces, Instituto de Biología de la Universidad Nacional Autónoma de México.
