Introduction

The wide range of morphological variation exhibited by Middle American cichlids is perhaps one of their most conspicuous characteristics, which has often been correlated to the great environmental heterogeneity that characterizes this region (^{López-Fernández et al., 2010}; ^{Říčan et al., 2016}). The Grijalva-Usumacinta River basin is shared between Guatemala and Mexico and likely one of the best examples of a watershed with such heterogeneity in habitats and environments. It has a complex and dynamic geological, volcanic, orographic, and physiographic history with wide variations in climate, vegetation, and aquatic ecosystems (^{Kolb & Galicia, 2012}; ^{Vaca et al., 2019}). These conditions have been attributed as major contributors to the high fish diversity found in this watershed (^{Matamoros et al., 2015}; ^{McMahan et al., 2017}). The Grijalva-Usumacinta River basin contains over 30 species of cichlids, the highest concentration of species within this family in Middle America, which exhibit great morphological and functional diversity. Additionally, this basin also has the highest number of endemic (^{Matamoros et al., 2015}) and non-endemic species of cichlids with restricted geographic distributions (^{Gómez-González et al., 2018}; ^{Říčan et al., 2016}; ^{Soria-Barreto et al., 2019}). The exceptional species richness and morphological variation in this region have historically led to a series of taxonomic and systematic problems reflected in several phylogenetic hypotheses and taxonomic changes in this group (^{López-Fernández et al., 2010}; ^{McMahan et al., 2015}).

Among Central American cichlids, the genus Vieja is one of the most taxonomically complex. Its widespread distribution and high levels of morphological variability have historically prevented a clear understanding of its evolutionary history (^{Gómez-González et al., 2018}; ^{McMahan et al., 2010}, ²⁰¹⁵, ²⁰¹⁹; ^{Říčan et al., 2016}). It was only recently (see ^{McMahan et al., 2015}) determined that this genus comprises 8 valid species whose distribution extends from the Papaloapan River basin in Mexico to the Chagres River in Panama on the Atlantic slope, and from the Isthmus of Tehuantepec in Mexico to Lake Coatepeque in El Salvador on the Pacific slope (^{Fricke et al., 2020}; ^{Gómez-González et al., 2018}; ^{McMahan et al., 2015}, ²⁰¹⁷; ^{Říčan et al., 2016}).

Although higher-level relationships of Vieja and its close relatives have largely been resolved through molecular and morphological studies, species relationships, limits and distributions remain problematic, especially among closely related species (^{McMahan et al., 2010}, ²⁰¹⁵; ^{Říčan et al., 2016}). While morphological characters frequently show high levels of overlap between species, there are little or no reliable diagnostic characters in the original species descriptions and the available taxonomic keys present ambiguous characters (^{Gómez-González et al., 2018}; ^{McMahan et al., 2019}). The color patterns of live specimens have often been used to distinguish species; however, this practice also yields ambiguous results due to specimen preservation, geographic and ontogenetic variation, and ecophenotypic plasticity (^{Miller et al., 2005}). Therefore, the incorporation and description of additional morphological characters in comparative and phylogenetic analyses are highly recommended to achieve improved taxonomic clarification and distinction among species of Vieja (^{Gómez-González et al., 2018}; ^{McMahan et al., 2019}; ^{Soria-Barreto et al., 2011}).

Morphological characters related to feeding are useful discriminant characters in cichlids (^{Trewavas, 1983}; ^{Witte & van Oijen, 1990}). Similar to other groups of fishes, the ecological diversification of cichlids has involved resource partitioning. In this process, the adaptive capacity of trophic morphology is crucial and includes form, function, and behavior (^{Burress, 2015}). The trophic morphology of cichlids is closely related to the exploitation and partitioning of food resources, which are potential mechanisms for ecological divergence (^{Albertson et al., 2003a}, ^b; ^{Hulsey, 2009}). Several functional morphological studies have shown that the shape of bone structures of the trophic apparatus can determine biomechanical performance in relation to the capture and processing of food (^{Albertson & Kocher, 2001}; ^{Albertson et al., 2003a}; ^{Meyer, 1989}; ^{Waltzek & Wainwright, 2003}). The traits of these feeding structures have been used to taxonomically discriminate several genera and species of fish, including cichlids (^{Kassam et al., 2004}; ^{Trewavas, 1983}; ^{Witte & van Oijen, 1990}; ^{Witte & Witte-Maas, 1987}). Moreover, features such as the shape and number of teeth in the lower pharyngeal jaw (LPJ) and oral jaw have traditionally been included in descriptions of Neotropical and African cichlids (^{Alonso et al., 2019}; ^{Oliver, 2018}).

The closely related species V. bifasciata (^{Steindachner, 1864}), V. breidohri (^{Werner & Stawikowski, 1987}), and V. hartwegi (^{Taylor & Miller, 1980}) possess interesting phenotypes for the study of trophic morphology. Molecular phylogenetic analyses have shown few genetic differences among these 3 species (^{Gómez-González et al., 2018}; ^{McMahan et al., 2010}; ^{Říčan et al., 2016}). Although these species are very similar, certain differences in color pattern, body shape, and pharyngeal teeth have been observed (^{Gómez-González et al., 2018}). Vieja breidohri is sympatric with V. hartwegi in the middle and upper reaches of the Grijalva River basin, while V. bifasciata is widespread in the entire Usumacinta River basin and does not coexist with its close congeners (^{Miller et al., 2005}).

Since sympatric cichlids frequently exhibit trophic characteristics related to their diversification (^{Burress, 2015}, ²⁰¹⁶), differences in trophic anatomy are expected among these closely related species, especially in structures that are functionally important to food capture and processing. The main goal of this study was to analyze patterns of variation in the bones that form the oral jaws and LPJs of V. bifasciata, V. breidohri, and V. hartwegi to evaluate their utility as discriminating characters for species delimitation. Additionally, we explored possible functional implications of observed morphological variation, particularly among sympatric species.

Materials and methods

For the selection of specimens for analysis, we considered their number, size, and conservation status in the fish collection. We also attempted to compare the same number of specimens and those of similar size.

A total of 55 adult specimens of the genus Vieja were obtained from the fish collection of El Colegio de la Frontera Sur (ECOSC), including V. bifasciata (n = 20, standard length, SL 62.4 - 148 mm, ECOSC 55, 731, 1621, 2381, 2387, 2675, 2746, 2760, 3353, 3362, 3878, 7509), V. breidohri (n = 21, SL 61.07 - 115.35 mm, ECOSC 4584, 6835, 6840, 6868), and V. hartwegi (n = 14, SL 59.38 - 108.29 mm, ECOSC 6838, 7455). All V. breidohri and V. hartwegi specimens were collected from the Belisario Domínguez reservoir and Lagos de Colón, located in the upper reaches of the Grijalva River in Chiapas, Mexico (Fig. 1). Moreover, all V. bifasciata specimens originated from localities in the middle reaches of the Usumacinta River. Their taxonomic determination was corroborated following ^{Miller et al. (2005)} and their original descriptions (^{Gómez-González et al., 2018}; ^{Steindachner, 1864}; ^{Taylor & Miller, 1980}; ^{Werner & Stawikowski, 1987}). Due to high levels of overlap in morphometric and meristic characters, the key characters used to separate these species are basically body color patterns (^{Gómez-González et al., 2018}). Preserved specimens of V. bifasciata have a curved midlateral stripe extending from the caudal fin base to the opercle, running at the mid-body level. In V. hartwegi, the longitudinal stripe is straight and has a dark blotch at the caudal fin base. Vieja breidohri has an interrupted midlateral stripe with dark blotches extending from the caudal fin base to the opercle.

Figure 1 Study area and sampling localities of examined material. 

Specimens were cleared and stained following ^{Taylor & Van Dyke (1985)}. The premaxilla, dentary, and articular from the oral jaw, as well as the LPJ (Fig. 2), were dissected following standardized protocols. We selected these bones because they are important in trophic morphology and have taxonomic value. The maxilla was discarded because it was difficult to recognize homologous points to define landmarks and describe their shape. Each bone was photographed with a reference scale using a Canon 70D camera mounted on an EsteREO Discovery V12 Zeiss microscope. Only LPJs were photographed dorsally. For V. breidohri, 19 premaxillas and 12 LPJs were available, while 18 premaxillas and 17 LPJs were available for V. bifasciata. Broken bones were removed from the analysis.

Figure 2 Lateral view of the oral jaws (a) and dorsal view of the lower pharyngeal jaw (LPJ) (b) of generalized shape for the genus Vieja. 1, Ascendent process of premaxilla; 2, dentigerous process of premaxilla; 3, ascendent process of dentary; 4, dentigerous area of dentary; 5, ventral process of dentary; 6, coronoid process; 7, articular facet; 8, horizontal process; 9, lateral process of LPJ; 10, dentigerous area of LPJ; 11, anterior process of LPJ. 

We defined 2D landmarks and semi-landmarks to describe the shape of each bone using MakeFan8 (^{Sheets, 2014}) and TpsDig2 ver. 2.26 software (Rolf, 2016). The landmark and semi-landmark configurations for each bone were as follows: premaxilla (3 landmarks, 21 semilandmarks), dentary (6, 10), articular (5, 16), and LPJ (7, 10) (Figs. 3-6). We used open curves in each bone to describe the shape of the contours between landmarks: premaxilla (3 curves), dentary (2 curves), articular (2 curves) and LPJ (3 curves). Semi-landmarks were aligned by sliding points along the curves using the minimum Procrustes distance criterion in SemiLand8 software (^{Sheets, 2014}). To remove the effects of scale, orientation, and position, the landmark and semi-landmark configurations were superimposed using the generalized Procrustes analysis in CoordGen8 software (^{Sheets, 2014}).

Figure 3 CVA ordination of Vieja species based on premaxilla shape. Deformation grids are associated to the most negative and positive values of the first 2 factors (CV1 and CV2). Premaxilla with landmarks (black dots) and semi-landmarks (white dots) used in the analysis. 

Figure 4 CVA ordination of Vieja species based on dentary shape. Deformation grids are associated to the most negative and positive values of the first 2 factors (CV1 and CV2). Dentary with landmarks (black dots) and semi-landmarks (white dots) used in the analysis 

Figure 5 CVA ordination of Vieja species based on articular shape. Deformation grids are associated to the most negative and positive values of the first 2 factors (CV1 and CV2). Articular with landmarks (black dots) and semi-landmarks (white dots) used in the analysis 

Figure 6 CVA ordination of Vieja species based on LPJ shape. Deformation grids are associated to the most negative and positive values of the first 2 factors (CV1 and CV2). LPJ with landmarks (black dots) and semi-landmarks (white dots) used in the analysis 

To evaluate the effect of allometry on intra- and interspecific bone shape variation, we performed linear univariate regression with the Procrustes distances of each specimen as the dependent variable and the logtransformed values of bone shape centroid size as the independent variable (^{Zelditch et al., 2012}). This analysis tested the null hypothesis that there is no relationship between size and shape. The analysis was executed using Regress8 software with 100 permutations (^{Sheets, 2014}). Superimposed landmarks and semi-landmarks were aligned and used to obtain partial warps with the 3 smallest specimens as references. The significance of the regression parameters was estimated through bootstrap permutations (1,000 replicates). We used the slope value (m) and explained variance (%) to determine the influence of size on bone shape. Lower values of regression slope (m) and explained variance (%) indicate no linearity between both variables and a weak effect of size on shape (^{Zelditch et al., 2012}).

Partial warps were used in a canonical variate analysis (CVA) and multivariate analysis of variance (MANOVA) to view and quantify the intra- and interspecific bone shape differences in morphospace using CvaGen8 software (^{Sheets, 2014}). Bartlett’s test was performed to determine the significance of CVA scores based on Wilk’s lambda (λ) values. Bone shape variation was visualized using deformation grids associated with the canonical variates. Partial Procrustes distances between the mean shapes of groups were calculated to perform paired comparisons. The interspecific differences of each bone were assessed using Goodall’s F test with Bonferroni correction and 100 permutations using Twogroup8 software (^{Sheets, 2014}).

Finally, we tested for significant differences in the number of teeth in the LPJ to determine the effect of LPJ size on this number. The number of teeth in LPJs were counted and an analysis of covariance (ANCOVA) was performed with species as the independent variables and the number of teeth in the LPJ as the dependent variable, with centroid size as a covariate. This procedure was implemented in SPSS 15.0 (2006).

Results

The allometric analysis did not show a significant relationship between size and bone shape. At the intraspecific level, the values of the regression slope (m) were low and the percentage of shape variation explaining size (E) varied from 6.03 to 36.98%. Premaxilla: V. bifasciata (m = 0.044, E = 33.11%), V. breidohri (m = -0.019, E = 36.98%), and V. hartwegi (m = -0.003, E = 24.76%). Dentary: V. bifasciata (m = 0.03, E = 12.04%), V. breidohri (m = 0.004, E = 7.6%), and V. hartwegi (m = 0.004, E = 14.71%). Articular: V. bifasciata (m = -0.008, E = 13.48%), V. breidohri (m = -0.003, E = 17.3%), and V. hartwegi (m = -0.007, E = 8.34%). LPJ: V. bifasciata (m = 0.004, E = 6.03%), V. breidohri (m = 0.005, E = 23.91%), and V. hartwegi (m = 0.022, E = 9.44%). At the interspecific level, we detected relatively low linearity values in the regression analyses and in the variation explained for the articular (m = -0.01, E = 7.85%), dentary (m = 0.013, E = 6.09%), premaxilla (m = 0.0004, E = 21.02%), and LPJ (m = -0.009, E = 8.38%).

Although the CVA revealed significant differences in shape among the 4 bones in the 3 studied species (Table 1, p < 0.05), the obtained Wilks lambda (λ) values were relatively low. The premaxilla and LPJ showed the greatest differentiation among species (Figs. 3, 6) and explained variance (59.8 and 20.48%, respectively) (Table 1). The dentary exhibited significant differences in the first 2 axes (Fig. 4) with less explained variance (5.52%; Table 1). The articular only recovered a significant CV1 explaining 10.19% of the variance (Fig. 5). Paired comparisons between species also indicated significant differences in the shape of the 4 bones (Table 2).

Deformation grids showed changes in the size and orientation of the ascending process of the premaxilla as well as the height of the dentigerous process. Vieja breidohri was recovered in the negative axis of CV1 and their ascending process was elongated and thin at its base, while the dentigerous process was wider at the anterior portion. The positive axis contained V. bifasciata, in which the ascending process was short and wide at its base and the dentigerous process was narrow toward its posterior part. Vieja hartwegi possessed intermediate characteristics (Fig. 3).

In CV2, changes were associated with the shape of the ascending process and the height of the dentigerous process. Vieja hartwegi was located along the negative axis, where the ascending process was slightly longer and oriented posteriorly while the anterior section of the dentigerous process was higher with its posterior part inclined dorsally (Fig. 3).

In the dentary, changes were associated with the length of the ascending process and the heights of the dentigerous area and ventral process. Vieja breidohri and V. hartwegi were recovered along the negative axis of the CV1 and showed the shortest ascending process and the highest dentigerous area and ventral process. Vieja bifasciata was located along the positive region of CV1 with the longest ascending process and relatively short dentigerous area and ventral process. At the negative axis of CV2, V. breidohri was differentiated from V. hartwegi by a lower height of the dentigerous area, with V. bifasciata located between them (Fig. 4).

For the articular, the most notable changes were observed in the length of the coronoid process, height of the horizontal process, and size of the articular. Vieja hartwegi was located at the negative axis of CV1 with the highest coronoid process, widest horizontal process, and narrowest articular facet. Vieja bifasciata was located at the positive end of CV1 with the lowest coronoid process, narrowest horizontal process, and widest articular facet. Vieja breidohri possessed intermediate articular characteristics (Fig. 5).

The LPJ showed changes in the length of the anterior process, the position of the lateral processes, and the size of the dentigerous area. Vieja breidohri was in the negative axis of CV1 with the lateral processes displaced, anterior process and anterior part of the dentigerous area lengthened, and posterior lateral edges thinned. Vieja hartwegi was at the positive end of CV1 and showed lateral processes displaced forward, anterior process shortened, and posterior lateral edges of the dentigerous area widened. Vieja bifasciata had an intermediate form. In CV2, the most important change was in the shape of the dentigerous area. Vieja breidohri and V. hartwegi were in the negative region, showing the dentigerous area with the widest posterior lateral edges and the shortest anterior border. Vieja bifasciata was located in the positive region, with the narrowest dentigerous area at the back and an elongated front edge (Fig. 6).

Based on the analysis of covariance (ANCOVA), no significant effect of centroid size was found (p > 0.05, F = 2.387) upon comparing the number of teeth in the LPJ. Significant differences in V. bifasciata were observed when compared to V. breidohri and V. hartwegi (p < 0.05, F = 12.96). Vieja bifasciata had 130 to 210 teeth (x-= 170), V. breidohri had 103 to 178 (x-=129), and V. hartwegi had 105 to 172 (x-= 143) (Fig. 7).

Figure. 7 Box plots showing maximum, median, minimum values of number of pharyngeal teeth in the LPJ. Quartiles correspond to 25%. 

Discussion

Our morphometric analysis showed that the shape of the bones that make up the oral jaw and LPJ differ significantly among the 3 analyzed species of the genus Vieja. These findings support the taxonomic distinction of these 3 species and can function as discriminating characters in addition to observed differentiation in color pattern, body shape, and the type of pharyngeal teeth (^{Gómez-González et al., 2018}).

The results of the intra- and interspecific allometry tests suggest that there is no significant influence of allometry in this study. Frequently, a linear relationship between both variables can affect the morphometric comparison between species or groups. Therefore, it is very important to evaluate such relationships through this type of study and remove or account for the effect of allometry in analyses (^{Zelditch et al., 2012}).

The shape and variation of the premaxilla explained the greatest variation among species, with the most notable differences observed in the length, width, and orientation of the ascending process. Although the length of the ascending process is functionally associated with an increase in the protrusion distance from the mouth, it is coupled with a decrease in bite force (^{Barel, 1983}; ^{Burress et al., 2020}). These functional considerations, as well as the variation observed in the premaxilla, allow us to predict that V. bifasciata has the least amount of protrusion of the mouth and V. breidohri has the greatest. Moreover, V. hartwegi should possess a functionally intermediate form.

The dentary and articular are fused to form the lower oral jaw. The dentary is located in the anterior position, supporting the oral teeth and participating directly in the capture of food (^{Albertson & Kocher, 2001}). In this bone, the greatest variation was observed in the height of the dentigerous area and the length of the ascending and ventral processes. All of these characters are functionally associated with bite force during food capture, with a greater bite force present when the dentigerous area is high and the length of the dorsal process is lower (^{Albertson et al., 2003b}). Therefore, it can be inferred that bite force is likely higher in V. breidohri and V. hartwegi than in V. bifasciata.

In the articular, the majority of variation was found in the length of the coronoid process and the size of the articular facet. The coronoid process works as a lever during the closing of the mouth, and there is a direct relationship between its length and bite force (^{Albertson & Kocher, 2001}; ^{Otten, 1981}). Under these conditions, V. hartwegi is expected to have the highest bite force. Regarding the shape of the articular facet, jaw support increases during the opening of the mouth when the ventral edge is extended (^{Albertson & Kocher, 2001}; ^{Otten, 1981}). Based on this relationship, we expect that V. breidohri and V. bifasciata, which exhibit wider articular facets, have greater jaw support than V. hartwegi.

LPJ shape exhibited a high level of variation, with differences in the orientation of the lateral processes and the length of the anterior process. Functionally, the pharyngeal jaw processes food before it passes to the digestive tract (^{Barel, 1983}; ^{Burress et al., 2020}; ^{Hulsey et al., 2019}; ^{Liem, 1973}). The shape of this structure is related to the bite force of the pharyngeal jaw, which increases when the distance between the lateral processes is greater and the anterior process is short (^{Burress et al., 2019}; ^{Muschick et al., 2011}). According to our results, and from a functional perspective, V. hartwegi should have the strongest pharyngeal bite, followed by V. bifasciata, while V. breidohri should have the weakest.

Finally, the number and shape of pharyngeal teeth play a fundamental role in cichlid feeding (^{Burress, 2016}; ^{Liem, 1973}). Vieja bifasciata was observed to have more teeth, while the differences between V. breidohri and V. hartwegi were not significant in this regard. All species of Vieja have been classified as omnivores (^{Říčan et al., 2016}). Omnivorous species typically possess more pharyngeal teeth, which are usually used to manipulate and select food using the pharyngeal jaw (^{Burress, 2016}). Therefore, we expect that V. bifasciata can manipulate food more efficiently than V. hartwegi and V. breidohri. Besides the number of teeth, another important aspect is the shape of the teeth and their arrangement in the LPJ. Notably, teeth can have various shapes and sizes. Those located on the back of the LPJ are often stronger and more specialized, while those on the front of the LPJ are used to grasp and transport food (^{Burress, 2016}; ^{Hellig et al., 2010}; ^{Liem, 1973}). The shape and arrangement of these teeth were previously confirmed in a qualitative review of the LPJ (^{Gómez-González et al., 2018}; ^{Taylor & Miller, 1980}). On the LPJ, teeth are conical in shape with the cusp facing the front and a second small cusp at the base. The anteriorly located teeth are small and increase in size toward the back, while the medial teeth are arranged in series with the posterior teeth being noticeably larger and more robust (^{Gómez-González et al., 2018}).

Inferences regarding bite force can be tested via biomechanical analysis (^{Hulsey et al., 2005}), but also, they can be explained by information on the diet and feeding behavior of species. The diet of V. bifasciata includes soft items such as detritus, plants, and algae (^{Pease et al., 2018}; ^{Soria-Barreto et al., 2019}); in contrast, the diet of V. hartwegi includes plants and invertebrates (^{Conkel, 1993}). Although the diet of V. breidohri remains unknown, the type of dentition on its LPJ suggests a diet based on harder prey such as mollusks and other invertebrates (^{Konings, 1989}). To support this, we reviewed the stomachs contents of the fishes used in this analysis. The diet of V. breidohri in Lagos de Colón is composed mainly of mollusks, with a smaller proportion of seeds and fish scales. This information confirms that this species is omnivorous with a tendency to consume mollusks. The shape and type teeth of the LPJ are the most conspicuous morphofunctional traits related with dietary tendencies.

The differences in osteological characters observed in the present study contradict the scarce genetic differentiation reported for these species, especially between V. breidohri and V. hartwegi (^{Gómez-González et al., 2018}; ^{McMahan et al., 2010}; ^{Říčan et al., 2016}). Under this scenario, these morphological differences may represent a case of trophic polymorphism similar to that of the widely documented case of the endemic species Herichthys minckleyi (^{Cohen et al., 2005}; ^{Hulsey et al., 2005}, ²⁰⁰⁶; ^{Sage & Selander, 1975}) and the Amphilophus spp. complex (^{Barluenga & Meyer, 2010}; ^{Muschick et al., 2011}), where trophic morphology is related to various ecological factors that include the distribution of food resources.

The LPJ is one of the most commonly analyzed structures in different groups of cichlids because it has diversified under sympatric conditions and its variation has been reflected in trophic adaptations that have favored the differential use of food resources and segregation of ecological niches (^{Burress, 2016}; ^{Takahashi & Koblmuller, 2011}). In this sense, our results for V. breidohri and V. hartwegi are remarkably congruent with those reported for polymorphic species of the genus Amphilophus and Herichthys, which have no substantial differences in the structures that make up the oral jaw but exhibit differences in the form and dentition of the LPJ (^{Burress, 2016}; ^{Hulsey, 2006}; ^{Hulsey et al., 2005}, ²⁰¹⁶; ^{Muschick et al., 2011}; ^{Pérez-Miranda et al., 2019}).

In this study, the degree of differentiation and directions of change in the shape of the 4 bones do not have a defined pattern. A lack of integration between the components of both jaws may result from the modular organization of trophic structures (^{Burress et al., 2020}; ^{Cooper et al., 2011}). While the oral jaw has the functional purpose of capturing food, the pharyngeal jaw processes food before it moves on to the digestive tract (^{Barel, 1983}; ^{Burress et al., 2020}; ^{Hulsey et al., 2019}; ^{Liem, 1973}). This functional organization has also been observed in other parts of the skull for some cichlids (^{Parsons et al., 2018}), even between the structures that control the opening and closing of the jaw (^{Albertson et al., 2005}; ^{Parsons et al., 2011}, ²⁰¹²).

Vieja bifasciata generally exhibits greater morphological differences, which are likely associated with geographic isolation and ecological interactions with other cichlids that share their area of distribution. In the case of the sympatric species V. breidohri and V. hartwegi, gradual differentiation was observed that sometimes overlaps. This morphological overlap has also been observed in the skulls of Herichthys species, where it is considered the result of selective and genetic restrictions. However, head shape differences have been observed between sympatric sister species (^{Pérez-Miranda et al., 2020}).

In cichlids there is evidence of trophic structures from oral and mandibular jaws that have an underlying genetic basis (^{Fruciano et al., 2016}). Making interpretations about the taxonomic value of these characters could be risky since their expression could result from phenotypic plasticity linked to ecological factors such as the availability of food resources or habitat preferences (^{Gómez-González et al., 2018}; ^{Říčan et al., 2016}). Therefore, we support the need for further morphofunctional and ecomorphological studies to assess relationships between diet and the phenotypic expression of the trophic apparatus in cichlid fishes (^{Feilich & López-Fernández, 2019}). This is particularly important for V. hartwegi, which shows enormous morphological variation throughout its geographic distribution and during growth (^{Gómez-González et al., 2018}). Likewise, it is necessary to perform further genetic studies to assess the degree of differentiation between species and its effect on the phenotypic expression of oral and pharyngeal jaws.

From a biomechanical perspective, minor morphological differences can affect the capture and processing of food (^{Barel, 1983}; ^{Cooper et al., 2010}; ^{Kassam et al., 2004}). The trophic anatomy of cichlids is one of the most complex and interesting functional systems since it can be highly variable and responds to multiple environmental and historical factors (^{Burress, 2016}).