This study assesses the patterns of variation in body shape, and relations of morphological similarity among species of the Theraps–Paraneetroplus clade in order to determine whether body shape may be a trait in phylogenetic relations. A total of 208 specimens belonging to 10 species of the Theraps–Paraneetroplus clade were examined. The left side of each specimen was photographed; in each photograph, 27 fixed landmarks were placed to identify patterns in body shape variation. Images were processed by using geometric morphometrics, followed by a phylogenetic principal component analysis. The phylogenetic signal for body shape was then calculated. To determine the relations in morphological similarity, a dendrogram was created using the unweighted pair group method and arithmetic mean values, while a Procrustes ANOVA and post-hoc test were used to evaluate significant differences between species and habitats. We found three morphological groups that differed in body length and depth, head size, and the position of the mouth and eyes. The body shape analysis recovered the morphotypes of seven species, and statistical differences were demonstrated in eight species. Based on traits associated with cranial morphology, Wajpamheros nourissati (Allgayer, 1989) differed the most among the species examined. No phylogenetic signal was found for body shape; this trait shows independence from ancestral relatedness, indicating that there is little congruence between morphological and genetic interspecific patterns. As evidenced by the consistently convergent morphology of the species in the Theraps–Paraneetroplus clade, the diversification of the group is related to an ecological opportunity for habitat use and the exploitation of food resources. Although no phylogenetic signal was detected for body shape, there appears to be an order associated with cranial morphology-based phylogeny. However, it is important to evaluate the intraspecific morphologic plasticity produced by ecological segregation or partitioning of resources. Therefore, future morphological evolutionary studies should consider cranial structures related to the capture and processing of food.

diversification, geometric morphometrics, morphological convergence, phylogenetic signal, Usumacinta province

Among Neotropical freshwater fishes, evidence of diversity suggests that allopatric speciation models frequently apply to several clades and that there are few cases of sympatric speciation stemming from adaptive processes (Albert et al. 2020). The main historical processes include river capture and sea-level oscillations, which fragment and merge fluvial networks. In this scenario involving geographic changes and ecological heterogeneity, phenotypic variation has been an important attribute in morphological diversification and environmental adaptation (Albert et al. 2020). Freshwater fish orders such as Siluriformes, Characiformes, Cyprinodontiformes, Gymnotiformes, and Cichliformes are the best examples of morphological diversity in the Neotropical region due to their high species richness and abundance (Albert et al. 2020; Elías et al. 2020).

In Neotropical cichlids, the ability to use new or newly available resources (i.e., ecological opportunity) has been an important mechanism in diversification (Arbour and López-Fernández 2016; Říčan et al. 2016). Studies related to morphological diversity in South American cichlids have demonstrated that body shape and size variation have been the main axes of diversification, most notably the constant presence of morphological convergence between lineages (López-Fernández et al. 2010, 2013). Morphological convergence is interpretable as evidence that natural selection has selected similar traits, thus providing strong evidence for the adaptive quality of said traits (Elmer and Meyer 2011; Losos 2011; Burress 2015).

Despite advances in knowledge regarding the evolution of diverse groups of cichlids worldwide, there are still lineages with incipient research, such as Middle American heroine cichlids. Middle America harbors approximately 124 cichlid species (Říčan et al. 2011; Matamoros et al. 2015) and includes areas considered to be centers of endemism and high diversity, such as the San Juan and Usumacinta ichthyological provinces (sensu Říčan et al. 2016). From evolutionary evidence, it has been assumed that their diversification was promoted by ecological opportunity and resource partitioning (López-Fernández et al. 2012; Burress 2016; Říčan et al. 2016). This is supported by the diversity of body shapes and the specialization of trophic anatomy, particularly of the oral and pharyngeal jaws (Liem 1973; Meyer 1993; Salzburger 2009; Burress 2016). Notably, this has promoted frequent cases of morphological convergence and the low phylogenetic signal of diagnostic characters. This is the main reason for the unclear and complex taxonomy of cichlids in Middle America (Stiassny 1991; Říčan et al. 2008, 2016; McMahan et al. 2013).

Among Middle American cichlids, the Theraps–Paraneetroplus clade (sensu Říčan et al. 2016) is notable due to the presence of species therein with highly variable and frequently convergent morphology, which is most evident in their body shapes and characteristics associated with food capture (Soria-Barreto and Rodiles-Hernández 2008; Soria-Barreto et al. 2011, 2019). This group of fish is estimated to have originated under sympatric conditions approximately 7.3 mya (Miller et al. 2005; Říčan et al. 2016) and is thought to result from the event of ancient adaptive radiation (Arbour and López-Fernández 2016; Albert et al. 2020). The Theraps–Paraneetroplus clade includes 25 species belonging to 10 genera. Moreover, the distribution of this clade is located in the Usumacinta ichthyological province, comprising the hydrological basins of Papaloapan, Coatzacoalcos, Grijalva, Usumacinta, and northern Belize (Říčan et al. 2016).

In addition to ecomorphological evidence indicating that the phenotypic expression of morphological attributes in some clade members is associated with habitat type and feeding (Soria-Barreto et al. 2019), studies on the systematics and evolution of Middle American cichlids have demonstrated the existence of convergent morphological characters between several species of the Theraps–Paraneetroplus clade (López-Fernández et al. 2014; McMahan et al. 2015; Říčan et al. 2016). As the same body shape patterns are recurrent among species that exploit similar habitats (McMahan et al. 2015; Říčan et al. 2016), the existence of lentic and lotic ecomorphological patterns has been previously proposed; species of the genera Cincelichthys, Kihnichthys, Oscura, and Vieja represent the lentic ecomorphotype (i.e., with short and deep bodies), and the species Theraps, Wajpamheros, Chuco, Rheoheros, and Paraneetroplus represent the lotic ecomorphotype (i.e., with elongated and slender bodies). Notably, both ecomorphotypes are present in the genus Maskaheros (see Říčan et al. 2016).

In the Theraps–Paraneetroplus clade, convergent morphological characteristics seem to support the hypothesis of diversification via ecological opportunity and resource partitioning, which contrasts with the hypothesis proposed by phylogenetic systematics and the theory of evolutionary non-independence (Felsenstein 1985). In the non-independence hypothesis, it would be expected that body shapes within the species and genera of this clade would have a diversification pattern similar to that of phylogeny and taxa sharing an ancestor that is most morphologically similar. In this way, species with a lotic body shape should share an ancestor, which should be similarly true for the lentic body shape. In contrast, if the non-independence hypothesis is rejected, then the morphological patterns should not be statistical dependents of the common ancestry (Revell et al. 2008). Therefore, the presently reported study aims to describe and compare the body shape variation patterns in 10 species of the Theraps–Paraneetroplus clade. Furthermore, the phylogenetic signal is obtained to measure the statistical non-independence of the morphologic trait values of the species due to their phylogenetic relatedness (Revell et al. 2008).

For this purpose, geometric morphometrics and comparative phylogenetic methods are used as analytical tools because they are commonly used to study the evolution of biological morphology. Geometric morphometrics can be used to identify variation in the pure shape of organisms, and separate the variation and size of individuals by analyzing shapes in multivariate space (Adams et al. 2004; Zelditch et al. 2004; Aguirre and Jiménez-Prado 2018). Notably, comparative phylogenetic methods can be used to analyze morphological characteristics and their significance in species diversification from a phylogenetic perspective (Pagel and Harvey 1988; Adams and Collyer 2018; Borges et al. 2019; Villalobos-Leiva and Benítez 2020). Recently, phylogenetic approaches have made it possible to understand the significance of morphological variation and changes in the dynamics of biological communities, particularly in speciation, adaptation, and extinction. This is essential to predict the effect of natural and anthropogenic changes on ecosystem processes. It also represents a fundamental step towards the management and conservation of biodiversity on the planet (Cavender-Bares et al. 2009).

To analyze the morphological variation among members of the Theraps–Paraneetroplus clade, the presently reported study included a total of 208 specimens (females and males of similar size) that correspond to 10 species, representative of each genus in the clade. All the specimens were deposited at the Fish Collection of El Colegio de la Frontera Sur, San Cristóbal (ECOSC). Theraps clade: Chuco intermedium (Günther, 1862) (abbreviation and number of specimens: Chin, n = 21); Cincelichthys pearsei (Hubbs, 1936) (Cipe, n = 24); Kihnichthys ufermanni (Allgayer, 2002) (Kiuf, n = 23); Theraps irregularis Günther, 1862 (Thir, n = 20); Wajpamheros nourissati (Allgayer, 1989) (Wano, n = 24); Paraneetroplus clade: Maskaheros argenteus (Allgayer, 1991) (Maar, n = 25); Oscura heterospila (Hubbs, 1936) (Oshe, n = 20); Paraneetroplus bulleri Regan, 1905 (Pabu, n = 3); Rheoheros lentiginosus (Steindachner, 1864) (Rhle, n = 25); Vieja hartwegi (Taylor et Miller, 1980) (Viha, n = 23).

Museum catalogue information. Chuco intermedium (Chin) ECOSC 103, 314(4), 334, 395(5), 440, 473(3), 815(2), 12747(2), 4892(2); Cincelichthys pearsei (Cipe) ECOSC 204, 229(2), 299, 300, 337, 444, 719(4), 849, 1049, 1512(3), 1055, 2352, 2546, 2575, 4422(2), 4436 (2); Kihnichthys ufermanni (Kiuf) ECOSC 90, 186, 233, 406, 409, 613, 675, 769, 1729, 1230, 1536(3), 1548, 1557(2), 1867, 1873, 2118(2), 2298, 4687, 7618; Maskaheros argenteus (Maar) ECOSC 386, 698, 741, 1280, 1448,1472, 1481, 1502, 1606, 1747, 1771, 1998, 2020, 2163, 2174, 2395, 2555, 2577, 4716, 4747(2), 4806(2), 4821, 7774; Oscura heterospila (Oshe) ECOSC 2338, 2720, 3053, 3054, 3491, 3505, 3777, 4563, 6709, 7826, 8465, 9070, 9080, 9267, 9318, 9816, 9849, 10164, 10165, 13757; Paraneetroplus bulleri (Pabu) ECOSC 12018(3); Rheoheros lentiginosus (Rhle) ECOSC 646, 853(3), 869(2), 1471(2), 1503(2), 1874, 1900, 2296(4), 2389, 2515(3), 2549, 2559, 7789, 4695, 12748; Theraps irregularis (Thir) ECOSC 245, 254, 817(2), 1255, 1780, 1967, 2133, 2626, 4725, 4729, 4809(9); Vieja hartwegi (Viha) ECOSC 4445(3), 4546, 6838, 6857(4), 7468, 7542(4), 7543(4), 7548(2), 7549(2), 12340; Wajpamheros nourissati (Wano) ECOSC 532(2), 684(2), 820(2), 893, 1237, 1288(2), 1546(2), 1847, 1289(2), 2082(2), 2105, 2280, 2651, 4744, 4888, 7336, 7453.

Morphometric analysis. Specimens were photographed on their left side using a Canon (EOS 70D) digital camera. The camera was mounted on a tripod to standardize the distance from the specimen. A 1-cm scale was placed on each photograph. To describe and compare body shapes, a geometric morphometric analysis was performed. In each photo, 27 fixed landmarks were placed using the configuration provided by Mejía et al. (2015) with two additional landmarks (Fig. 1). Image digitization and processing were performed using the software tpsUtil ver. 1.70 (Rohlf 2018) and tpsDig ver. 2.26 (Rohlf 2017).

Then, in order to eliminate variation caused by the size, rotation, and displacement of the specimens, a generalized Procrustes analysis (Goodall 1991) was performed (Aguirre and Jiménez-Prado 2018). The mean body shape configuration of each species was obtained in the same manner. In both cases, the “gpagen” function of the Geomorph ver. 4.0.0 library (Adams and Otárola-Castillo 2013; Adams et al. 2016) was used in R software (R Core Development Team 2017).

Figure 1. 

Location of fixed landmarks in species of the Theraps–Paraneetroplus clade (image modified from Mejía et al. 2015). 1. Anterior end of the lower maxilla, 2. Anterior end of the upper maxilla, 3. Length of the ascending premaxillary process, 4. End of the supraoccipital bone, 5. Start of the dorsal fin, 6. Last spine of the dorsal fin, 7. End of the dorsal fin, 8. Upper boundary of the caudal fin, 9. Center of the caudal fin, 10. Base of the caudal fin, 11. End of the anal fin, 12. Last spine of the anal fin, 13. Origin of the anal fin, 14. Origin of the pelvic fin, 15. Cleitral fusion, 16. Posterior end of the lower maxilla, 17. Posterior end of the upper lip, 18. Maximum point of curvature at the preoperculum, 19. Upper end of the preoperculum, 20. Upper end of the operculum, 21. Most posterior end at the operculum, 22. Dorsal insertion of the pectoral fin, 23. Ventral insertion of the pectoral fin, 24. Upper extreme of the sphenotic orbit, 25. Base of the sphenotic orbit, 26. Left extreme of the sphenotic orbit, 27. Right extreme of the sphenotic orbit.

Shape analysis. To reduce morphological variation related to phylogeny and differences in specimen size (allometry), regression of the Procrustes coordinates and centroid size was conducted using the “phyl.resid” function of the Phytools package ver. 0.7–80 in R software (Revell 2009, 2012). In order to then identify patterns in body shape variation, a phylogenetic principal component analysis (pPCA) was performed based on the Procrustes coordinates of the mean configurations of the 10 species. The nDNA molecular phylogeny based on ddRAD sequences proposed by Říčan et al. (2016) for Middle American cichlids was used in both instances. The body shape variation of the species was displayed on the first three pPCA axes. Additionally, deformation grids were obtained to visualize and describe the morphological variation among species in morphospace. All analyses were conducted in R software using the Phytools package ver. 0.7–80 (Revell 2009, 2012).

Additionally, the phylogenetic signal for body shape was computed by using the Kmult statistic (K) across 1000 permutations via the Geomorph package 4.0.0 in R software (Adams et al. 2021), K values <1 indicate a low phylogenetic signal, while K values >1 indicate a strong phylogenetic signal (Adams 2014). To determine the relation involved with morphological similarity among the 10 species, a dendrogram was constructed using the unweighted pair group method using mean values in Past 4.05 software (Hammer et al. 2001) based on the Mahalanobis distances obtained in MorphoJ 1.07a software (Klingenberg 2011).

Finally, to determine whether significant differences in body shape exist among species and between lentic (LE) and lotic (LO) habitats, a Procrustes ANOVA was used as implemented in Geomorph 4.0.0 (Adams et al. 2021) using the procD.lm function over 1000 permutations with the Procrustes coordinates and centroid size. Furthermore, to establish whether there are significant differences in body shape among species, a pairwise comparison test was performed using the pairwise function as implemented in library RRPP (Collyer and Adams 2018) over 1000 permutations.

The pPCA indicated that the first three components explained 80% of the total variance (PC1: 52.6%; PC2: 17.6%; PC3: 11.6%). Species located on the positive axis of PC1 include T. irregularis, P. bulleri, and R. lentiginosus; O. heterospila, M. argenteus, C. pearsei, and K. ufermanni were located on the negative axis; and C. intermedium, V. hartwegi, and W. nourissati were located on the middle axis. The deformation grids showed variation among species on the positive axis related to decreased body height, elongation of the caudal peduncle, a convex base of the anal fin, and a narrow distal section. Species on the negative axis showed a deep body height, shortened caudal peduncle, and concave base of the anal fin. According to the morphotypes, the lotic species T. irregularis, P. bulleri, and R. lentiginosus were on the positive axis, while lentic species on the negative axis included O. heterospila, M. argenteus, C. pearsei, and K. ufermanni. However, the lotic species C. intermedium and W. nourissati, as well as the lentic species V. hartwegi, were in the middle of the axis (Fig. 2A, B)

Figure 2. 

A) Phylomorphospace formed by PC1 and PC2. B) Phylomorphospace formed by PC1 and PC3. Black circles represent the mean body shape configuration for each species, and white circles the ancestral state. Black lines represent Theraps clade, and grey lines Paraneetroplus clade. LE represent lentic ecomorphotype, LO represent lotic ecomorphotype, LE/LO represent both ecomorphotypes. Deformation grids are associated to the most negative and positive values of the PC1 and PC2. Abbreviations: Chin-LO = Chuco intermedium; Cipe-LE = Cincelichthys pearsei; Kiuf-LE = Kihnichthys ufermanni; Thir-LO = Theraps irregularis; Wano-LO = Wajpamheros nourissati; Maar-LE/LO = Maskaheros argenteus; Oshe-LE = Oscura heterospila; Pabu-LE = Paraneetroplus bulleri; Rhle-LO = Rheoheros lentiginosus; Viha-LE = Vieja hartwegi.

In PC2, species found in the positive axis included C. pearsei, K. ufermanni, C. intermedium, T. irregularis, and V. hartwegi, while O. heterospila, W. nourissati, R. lentiginosus, M. argenteus, and P. bulleri were found on the negative axis (Fig. 2A). The deformation grids showed that the most remarkable deformation occurred in the cephalic region. Species on the positive axis exhibited heads with straight profiles and mouths in a terminal position, while the eyes were displaced posteriorly and slightly enlarged. Species on the negative axis showed rounded heads and a ventral mouth position, while the eyes were both smaller in size and were displaced anteriorly (Fig. 2A, 2B).

In PC3, W. nourissati was the most differentiated on the positive axis, exhibiting accentuated variation in the cephalic region with increased head size and a notable anteroventral displacement of the mouth. Additionally, the eyes and pectoral fins of W. nourissati were displaced posteriorly. On the negative axis, the remaining species were equally distributed with short heads, rounded profiles, and small mouths (Fig. 2B).

The dendrogram based on Mahalanobis distances showed that C. intermedium differed the most in body shape, followed by the C. pearsei and K. ufermanni groups. The remaining species formed two groups—one composed of P. bulleri, R. lentiginosus, and T. irregularis and another composed of M. argenteus and O. heterospila. The second most similar species were W. nourissati and V. hartwegi. The phylogenetic signal value of the Kmult statistic was 0.765, with a significance value of P = 0.308 and displaying no significant effect of phylogeny on body shape under the Brownian motion evolutionary model (Fig. 3).

Figure 3. 

Dendrogram based on Mahalanobis distances (Cophenetic correlation = 0.79). LE represent lentic ecomorphotype, LO represent lotic ecomorphotype, combined LE/LO represent both ecomorphotypes.

The Procrustes ANOVA identified significant differences in body shape among species (F = 34.62, R² = 0.59, P < 0.01) and habitats for the comparison of Procrustes coordinates (F = 56.24, R² = 0.20, P < 0.001). However, the comparison using the centroid size of each species failed to recover significant differences among species (F = 1.06, R² = 0.043, P= 0.181) and habitats (F = 1.20, R² = 0.005, P = 0.264). However, the pairwise comparison test showed statistical differences between all the species (P < 0.05) excluding V. hartwegi and P. bulleri (P = 0.055) (Table 1).

In the pPCA, species were placed in morphospace based on body length and depth, followed by variation in the cephalic region based on changes in head size and profile, mouth position, and eye position and size. Additionally, some morphological changes were evident based on the position and size of the fins. The variation observed in these morphological characteristics has been closely associated with the environments, locomotion, and feeding of fish (Wootton 1990; Muschick et al. 2012; Feilich 2016; Fugi unpublished¹). In African and South American cichlids, variation in these morphological characteristics has been fundamental in adaptive radiation events because it facilitates phenotypic and lineage diversification through the exploitation of ecological opportunities (Kocher et al. 1993; Cooper et al. 2010; Muschick et al. 2014; Arbour and López-Fernández 2016). The morphological body variation patterns observed in South American cichlids are similar to those observed in species of the Theraps–Paraneetroplus clade, which supports the hypothesis of Middle American cichlid diversification via ecological opportunity as proposed in previous works (Arbour and López-Fernández 2016; Feilich 2016; Říčan et al. 2016).

In the morphospace and similarity analysis, the most notable groupings were T. irregularis, R. lentiginosus, and P. bulleri, all of which presented elongated bodies and caudal peduncles as well as pelvic fins that were positioned ventrally. Ecomorphological studies have reported that these characteristics are functionally associated with high-velocity current environments (Lowe-McConnell 1991; Lauder and Tytell 2005; Pease et al. 2012; Feilich 2016). Elongated bodies are more hydrodynamic, while fins in the ventral position permit fish to maintain stability (Drucker et al. 2005; Lauder and Tytell 2005; Pease et al. 2012; Feilich 2016; Han et al. 2020). Soria-Barreto and Rodiles-Hernández (2008) reported the same morphological pattern for T. irregularis and R. lentiginosus in the Lacandon rainforest (Montes Azules Biosphere Reserve, Mexico), and Gómez-González et al. (2018) in V. hartwegi. Additionally, the mouth positions of the three species tend to be sub-terminal or ventral, which is associated with foraging for algae and invertebrates on the bottom, particularly on the surfaces of rocks with silt and sand substrates (Keast and Webb 1966; Miller et al. 2005; Artigas-Azas 2005b).

The group comprising O. heterospila, K. ufermanni, M. argenteus, and C. pearsei occupies another part of the morphospace. These species displayed deep bodies, short heads, shortened caudal peduncles, and mouths that were generally in a terminal position. These morphological characteristics are associated with environments where current velocity ranges from medium to slow, and have different types of substrates (i.e., rock, gravel, sand, and mud) (Berbel-Filho et al. 2016; Feilich 2016; Říčan et al. 2016). Although these species are generally herbivorous and omnivorous, some may be detritivorous (Miller et al. 2005; Soria-Barreto et al. 2019; Říčan et al. 2016). Furthermore, V. hartwegi, C. intermedium, and W. nourissati have an intermediate morphology with shallow elongated bodies and shortened caudal peduncles. These three species likely share ecological niches and inhabit areas with moderate-to-high velocities (Miller et al. 2005; Pease et al. 2012; McMahan et al. 2015; Říčan et al. 2016; Gómez-González et al. 2018). Additionally, these species present more ventral snouts, facilitating feeding on aquatic invertebrates, detritus, algae, and vegetation (Soria-Barreto et al. 2019).

Wajpamheros nourissati , the outlier of the cichlids analyzed, is the only species that has long and thick lips. This is a characteristic associated with feeding between cracks, on rock surfaces, and on the substrate (Winemiller et al. 1995; Artigas-Azas 2005a; López-Fernández et al. 2014; Říčan et al. 2016). This characteristic has been reported in Middle and South American cichlids, which feed on benthic components and live in habitats with a variety of substrates including sand, silt, and organic matter in fine and coarse particles (Barlow and Munsey 1976; Moreira and Zuanon 2002; Hahn and Cunha 2005; Elmer et al. 2010; López-Fernández et al. 2012). Additionally, W. nourissati has well-developed pectoral fins that provide stability in habitats with currents of moderate velocity (Allgayer 1989). These morphological traits can be functionally important for exploiting specific habitats and alimentary resources while reducing competition with other sympatric cichlids.

The variation associated with the cephalic region revealed that the majority of species occupy a position in the morphospace that corresponds to their respective clades. For example, members of the Theraps clade (i.e., C. pearsei, C. intermedium, K. ufermanni, and T. irregularis) are positioned toward the superior part of the morphospace, while members of the Paraneetroplus clade (i.e., M. argenteus, O. heterospila, P. bulleri, and R. lentiginosus) are located toward the lower region. However, convergence between the two clades was also observed, largely in the variation of the cephalic characteristics of V. hartwegi and W. nourissati. This notably coincides with the variation of the cephalic characteristic in members of the opposite clade in both cases. Comparative studies have revealed patterns of rapid diversification between lineages and phenotypes through habitat- and diet-related morphological diversity (Arbour and López-Fernández 2016; Feilich 2016). Furthermore, it has been shown that the early radiation of certain Neotropical cichlid tribes in South America occurred rapidly, which resulted in a pattern of conflicting divergence (López-Fernández et al. 2013; Astudillo-Clavijo et al. 2015).

Evidence of the low phylogenetic signal of the body shape thus disproves the hypothesis of evolutionary non-independence. However, as the value of the phylogenetic signal can be determined by several processes and evolutionary rates (Revell et al. 2008), it raises the question of what specific causes determine body shape. In the evolution of quantitative traits, there is evidence that a low phylogenetic signal is related to speciation events caused by divergent selection and adaptive processes (Revell et al. 2008). In several Neotropical cichlid clades, this is the main cause for diversification events related to ecological opportunity (López-Fernández et al. 2012, 2013; McMahan et al. 2015; Feilich 2016; Říčan et al. 2016).

In this context, some species of both clades (Theraps–Paraneetroplus) exhibited high convergence in the variation associated with body length and depth, fins position, and eye size and position. This is evident in C. intermedium, V. hartwegi, and W. nourissati, which showed similar body shapes despite being considered as different ecomorphotypes. The remaining species exhibited a morphological pattern consistent with what was expected from the corresponding habitat. The lentic species (C. pearsei, K. ufermanni, M. argenteus, and O. heterospila) showed a clear morphologically divergent pattern with respect to the lotic species (R. lentiginosus, T. irregularis, and P. bulleri). In the case of C. intermedium, V. hartwegi, and W. nourissati, it is probable that habitat preference is determined by other morphological traits that are likely linked to the cephalic region and associated with feeding behavior. This adaptive process has been documented in other fish groups; for example, in the Pomacentridae family, variations in cephalic traits explain feeding behavior and diet and are an important predictor of trophic habit (Aguilar-Medrano et al. 2011).

Based on the analysis of body shape and head characteristics, the lack of congruence in the order of species in morphospace does not fully support the ecomorphological classification described by Říčan et al. (2016). However, the results support the hypothesis of functional independence between the cephalic and postcranial regions identified for Middle American cichlids (Říčan et al. 2016). In regard to the cephalic region, the occurrence of five ecomorphotypes associated with feeding behavior was identified. Meanwhile, the postcranial region showed the lentic and lotic ecomorphotypes associated with their respective environments. In both instances, molecular phylogeny supported the recurrent evolution of diverse cephalic and postcranial ecomorphotypes among Middle American cichlids (Říčan et al. 2016). Additionally, modularity has been considered an evolutionary factor leading to patterns of variation among cichlids. Modularity studies in African cichlids have centered on modules associated with feeding: cichlid species that feed by suction have functional modules and those that feed by foraging have developmental modules (Parsons et al. 2012).

Species of the Theraps–Paraneetroplus clade show highly varied and convergent morphologies that are largely promoted by ecological opportunities associated with habitat and feeding preferences (Albertson and Kocher 2001; McKaye et al. 2002; Kassam et al. 2003; López-Fernández et al. 2013; Feilich 2016; Říčan et al. 2016). Body size and depth are features that appear unrelated to phylogenetic relations, whereas the head features of the majority of species exhibit patterns of variation associated with their phylogeny. However, further analyses of the shape, role, and evolution of these morphological attributes are necessary to understand their significance in the diversification of Middle American cichlids.

In this paper, the morphological variation of the Theraps–Paraneetroplus clade was divided into three groups within the morphospace. These variations are related to body length, body height, head shape, mouth position, and eye size and position. In particular, W. nourissati was the most divergent species due to its cranial characteristics, which are largely related to the shape of the head and the position of the mouth. Body size and height were the variables that best described the position of the species in the morphospace. This facilitates the recovery of lotic and lentic ecomorphotypes for seven species, with only C. intermedium, V. hartwegi, and W. nourissati showing incongruity. The presence of the phylogenetic signal disproves the non-independence hypothesis, but evidence suggests that body shape results from adaptive processes related to ecological opportunity. Although the results of the ANOVA with Procrustes distances and size of the centroid were contradictory, groupings in the morphospace and dendrogram were consistent with the ecomorphotypes and phylogeny. The pairwise comparison test showed statistical differences between all species, with the exception of V. hartwegi–P. bulleri. The morphological patterns found support the taxonomical validity of each species and can be used to describe body shape at the genus level. Future morphological evolution studies should consider cranial structures related to the capture and processing of food.

We thank the anonymous reviewers who provided helpful comments that improved the manuscript. Financial support for this study was received from the Project: “Conectividad y diversidad funcional de la cuenca del río Usumacinta” (Fondo de Investigación Científica y Desarrollo Tecnológico de El Colegio de la Frontera Sur, FID-784), coordinated by RRH. This manuscript is the result of the master thesis of YEAC.