The highly diverse ichthyofauna of northern Middle America has a complex biogeographic history. Frequent geological, volcanic, and climatic events from the Late Cretaceous to the Miocene and Pleistocene determined the diversification and distribution of fishes in the region (Elías et al. 2021). Biogeographic patterns and a large amount of endemism among freshwater species represent the clearest evidence of the hydrological history of northern Middle America (Albert et al. 2020; Elías et al. 2021).

The Grijalva–Usumacinta hydrological system provides an excellent model for understanding the effects of geological and climatic events on the evolution of fish communities in northern Middle America. This hydrological system is characterized by its diversity of fish species and large amount of endemism, which is predominantly observed in the Cichlidae and Poeciliidae families (Gómez-González et al. 2015; Velázquez-Velázquez et al. 2016; Soria-Barreto et al. 2018; Sánchez et al. 2019; Álvarez-Pliego et al. 2021). Constant events involving vicariance resulting from the capture of rivers and changes in sea level favored the diversification of several lineages through geographic isolation (Albert et al. 2020; Elías et al. 2021). The construction of dams on the Grijalva River can be considered an example to test the recent effects of hydrological and environmental changes on the dynamics of fish populations, particularly the consequences of a reduction or interruption in genetic flow for morphological and functional differentiation (Sánchez et al. 2019).

The lower regions of the Grijalva and Usumacinta basins share many species of the same ichthyofauna (Macossay-Cortez et al. 2011; Soria-Barreto et al. 2018). Nevertheless, the structure of the community in the upper regions of these rivers notably differs (Gómez-González et al. 2015; Velázquez-Velázquez et al. 2016; Soria-Barreto et al. 2018). Studies of systematic, biogeography and population genetics have revealed that the upper regions of the two basins have independent evolutionary histories (Elías et al. 2021; Beltrán-López et al. 2021; Terán-Martínez et al. 2021). The presence of relic endemic species such as the catfish of Chiapas (Lacantunia enigmatica Rodiles-Hernández, Hendrickson et Lundberg, 2005) and the livebearer (Xenodexia ctenolepis Hubbs, 1950) of the upper Usumacinta are evidence of the region’s unique biogeographic history (Rodiles-Hernández et al. 2005; Hrbek et al. 2007; Elías et al. 2021).

Although cichlids and poeciliids are the most diverse and abundant fish families within the Grijalva–Usumacinta system, other families also reflect the effects of the region’s historic events. Such is the case of the genus Dorosoma, for which genetic evidence shows that cryptic diversity exists throughout the distribution of the species Dorosoma petenense (Günther, 1867) in Middle America, which consists of several lineages (Elías et al. 2021). Two of these lineages converge in the Grijalva–Usumacinta basins. The lineage with the broadest distribution inhabits the region from the Río Grande to the lower Grijalva and Usumacinta rivers in Mexico and Guatemala, while the other is restricted to the upper Usumacinta River and Petén Itzá Lake in Guatemala (Elías et al. 2021). Curiously, in other Dorosoma species with similar distribution patterns, no significant genetic differences have been found (Elías et al. 2021). In this sense, the migratory behavior and salinity tolerance of each species could help to explain differences in genetic segregation and morphological differentiation (Bloom and Egan 2018). Dorosoma petenense is the most widespread and tolerant of different salinity ranges (Elías et al. 2021). In contrast, Dorosoma anale Meek, 1904 is considered a freshwater fish with little tolerance for salinity changes (Castro-Aguirre et al. 1999; Elías et al. 2021).

Based on the biogeographic and molecular precedents of the ichthyofauna in northern Middle America, we proposed an analysis of the intraspecific morphological variation of the shads Dorosoma anale and D. petenense throughout the Grijalva and Usumacinta rivers in Mexico. This study used linear biometric and geometric morphometric methods. Notably, both analyses are complementary and have been widely used in ichthyology to identify intra- and interspecific morphological differences and describe patterns of variation (Kerschbaumer and Sturmbauer 2011; Tripathy 2020; Lishchenko and Jones 2021). Morphological differentiation was expected due to historical vicariance events that occurred during river separation and sea level changes in northern Middle America. Furthermore, recent biological, environmental, and hydrological changes caused by human activities were expected to affect the phenotypic expression of morphological differences.

A total of 262 adult specimens were analyzed, corresponding to the species Dorosoma anale (n = 136, 71 males and 65 females) and D. petenense (n = 126, 44 males, 82 females). The specimens were deposited in the Fish Collection of El Colegio de la Frontera Sur (ECOSC; Table 1). The specimens originated from 18 sites distributed throughout the Grijalva and Usumacinta basins in Mexico (Fig. 1 and Table 2). For the Grijalva basin, specimens from sites located in the middle region were used. For the Usumacinta basin, specimens were selected from representative sites of the upper, middle, and lower regions, following the criteria proposed by Soria-Barreto et al. (2018; see Table 2). Sexual maturity was determined by reviewing the gonads.

Specimens were photographed from the left side of the body with a Sony DSC-HX300 digital camera (10 megapixels) using a 10 mm reference scale. To characterize the body shape, we used a configuration of 18 fixed landmarks (Fig. 2 and Table 2) digitalized using tpsDig2 software, version 2.16 (Rohlf 2015).

Morphometric and statistical analyses. To analyze geographic intraspecific variation, specimens were classified into four groups according to the collection site. One group consisted of specimens from the Grijalva basin, while the other three were from the upper, middle, and lower regions of the Usumacinta basin (Fig. 1 and Table 2).

An analysis of geometric morphometrics was conducted using MorphoJ software, version 1.07a (Klingenberg 2011). Based on the coordinates obtained from digitization, a generalized procrustes analysis was performed (Rohlf and Slice 1990; Dryden and Mardia 1998) to eliminate variation among the configurations of the reference points due to the effects of position, orientation, and specimen size. Additionally, to eliminate the effect of allometry due to variation in specimen size, multivariate regression was conducted based on the procrustes coordinates (shape variables) with the values of the logarithm of the centroid (size variable).

With the residual values of the multivariate regression, a principal components analysis (PCA) was performed to evaluate intraspecific variation. The first two principal components were used to explore the distribution of the specimens in the morphospace and describe variation in body shape based on the deformation grids. Later, we conducted a canonical variate analysis (CVA) to determine whether significant differences in body shape exist among the four groups. Additionally, we carried out paired comparisons based on the procrustes distances. Finally, we applied a discriminant function analysis (DFA) to perform cross-validation to determine the percentage of classification of the specimens in each group based on the Mahalanobis distances. All tests subjected the data to 10,000 permutations, when appropriate.

Additionally, based on the deformation grids, we identified the body sections for which greater variation existed. We then took linear measurements to evaluate whether they are discriminant among the four groups. Measurements were obtained by using the CoordGen8 program (IMP; Sheets 2014) with the photographs. The configuration of landmarks used in the geometric morphometric analysis followed the protocols of Hubbs and Lagler (1947), Whitehead (1985), and Farré et al. (2016). Measurements were standardized according to the method of Elliott et al. (1995) to eliminate the effect of allometry caused by variation in specimen size. Later, we conducted statistical analyses (ANOVA and Kruskal–Wallis) and a posteriori tests (Tukey, Mann–Whitney) to determine whether significant differences existed among the groups. For the measurements showing significant differences, box plots were elaborated to visualize their variation, which was expressed in proportions with respect to the standard length. For statistical analyses and the elaboration of box plots, we used the PAST program, version 3.14 (Hammer et al. 2001).

Museum catalog information. Catalogue number of the specimens used in the morphometric analysis. Dorosoma anale: ECOSC 612, 658, 1286, 1737, 3492, 4426, 6708, 6714, 10713, 10714, 11748, 11752, 12549, 12555, 12665, 12790, 13521, 13533 al 13535, 13546, 13549, 13561, 13564, 13565, 13976, 14290, 14300; D. petenense: ECOSC 7339, 8698–8707, 9882–9891, 12616, 12619, 12657, 12669, 12716, 12722, 13702, 13723, 13738, 14547, 14679, 14680.

Intraspecific variation in Dorosoma anale. In the PCA, the first two components explained 42.5% of the total variance. In the morphospace, no formation of groups was observed given the extensive overlap among specimens (Fig. 3A). Throughout PC1 (28.35%), greater variation was obtained for specimens from the upper region were the most broadly distributed throughout this component. Although much less variation occurred for PC2 (14.18%), it was not possible to distinguish groups (Fig. 3A). Based on the deformation grids for PC1, specimens in the negative axis were observed to have larger heads, shorter caudal peduncles, and anteriorly displaced dorsal and pelvic fins. For the positive axis, specimens had shorter heads, longer caudal peduncles, and dorsal and pelvic fins that were slightly more forward (Fig. 3A). For PC2, the most evident variations were in body height and the fact that specimens reached maximum body depth on the positive axis, while they tended to be less deep toward the negative axis (Fig. 3A).

Meanwhile, the CVA and paired tests revealed significant differences among groups (P < 0.05). Nevertheless, the only differences among the four groups were between the Grijalva group and the three groups of the Usumacinta basin. For the Usumacinta groups, only the middle and upper regions were significantly different (P < 0.05; Table 2). Cross-validation based on the DFA indicated that the assignment percentages were highest in the Grijalva group, while the highest percentages for the Usumacinta basin were found in the group of the upper region (Table 3).

Based on the variation in body shape observed in the deformation grids, we selected five linear measurements to evaluate their capacity to discriminate among the groups. The selected measurements were as follows: I) anterior margin of the upper mandible to the posterior margin of the operculum (landmarks 1–17); II) posterior margin of the supraoccipital crest to the anterior insertion of the anal fin (2–9); III) anterior insertion of the dorsal fin to the anterior insertion of the pelvic fin (3–10); IV) posterior insertion of the dorsal fin to the anterior insertion of the anal fin (4–9); V) anterior insertion of the anal fin to the posterior insertion of the anal fin (8–9).

In D. anale, the statistical analyses (ANOVA, Kruskal–Wallis) and respective a posteriori tests (Tukey, Mann–Whitney) revealed that only the following three measurements could discriminate at least one of the groups (P < 0.05): I) anterior margin of the upper mandible to the posterior margin of the operculum (1–17); III) anterior insertion of the dorsal fin to the anterior insertion of the pelvic fin (3–10); IV) posterior insertion of the dorsal fin to the anterior insertion of the anal fin (4–9). The box plot of these three measurements (expressed in proportions) allowed us to determine that the groups with the greatest variation were from the upper regions of the Usumacinta and Grijalva basin (Fig. 3B).

Intraspecific variation in Dorosoma petenense. In the PCA, the first two components explained 42.73% of the total variance. In PC1 (25.19%), we observed a substantial overlap of the four groups on the negative axis. However, on the positive axis, specimens from the upper Usumacinta appeared to diverge, especially from Site 4 (Miramar Lagoon, Fig. 4A). Although no clear separation was found among the groups in PC2 (17.53%), the majority of specimens from the upper Usumacinta were located toward the positive extreme (Fig. 4A). The negative end of PC1 represented specimens with relatively deep bodies and fins that were anteriorly displaced relative to specimens that possessed positive scores (Fig. 4). In the deformation grids of PC2, shallow-bodied specimens were in the negative axis, while deeper-bodied specimens were located toward the positive part (Fig. 4B).

The CVA and paired tests showed significant differences among the four groups (P < 0.05; Table 3. The greatest procrustes distances were obtained for the upper Usumacinta group (Table 2). The DFA corroborated that the specimens of the upper Usumacinta were the most different, while the cross-validation revealed that they had the highest percentage of allocation (Table 3).

Multivariate analyses (ANOVA, Kruskal–Wallis) of the five linear measurements and the corresponding a posteriori tests showed that four measurements are discriminant (P < 0.05): II) posterior margin of the supraoccipital crest to the anterior insertion of the anal fin (2–9); III) anterior insertion of the dorsal fin to the anterior insertion of the pelvic fin (3–10); IV) posterior insertion of the dorsal fin to the anterior insertion of the anal fin (4–9); V) anterior insertion of the anal fin to the posterior insertion of the anal fin (8–9). Statistical differences were observed in the groups of the upper Usumacinta and the Grijalva. The diagrams of the four measurements (expressed in proportions) suggest that the groups that differed most were from the upper Usumacinta and the Grijalva (Fig. 4B).

Little intraspecific morphological differentiation was observed within Dorosoma anale and D. petenense throughout their distribution in the Grijalva and Usumacinta basins. Nevertheless, patterns of variation and morphological differences were identified in some of the geographic groups, which allowed us to assume that some regional historic and/or ecological processes were involved in creating and maintaining the phenotypic differentiation in both species. In the Clupeiformes and other fish taxa, migratory behavior and tolerance to salinity could have important implications for morphological differentiation (Aguirre and Bell 2012; Bloom and Egan 2018; Jiménez-Prado and Aguirre 2021). While it was difficult to test this hypothesis based on our results, this possibility should not be discarded.

In the morphospace of D. anale, no separation among groups by geographic location was observed. However, among specimens of the upper Usumacinta basin, variation existed in terms of head size, body depth, and fin position. This could be related to the type of habitat since the separation was observed among specimens from the river (sites 1, 2, and 18) and lake habitats (sites 3 and 4). Additionally, comparisons among geographic groups based on the statistical tests showed that the specimens of the Grijalva basin were differentiated by having shallower bodies.

A similar pattern of variation was observed in the morphospace of D. petenense, while the overlap between the four geographic groups was also found. Nevertheless, specimens of the upper Usumacinta tended to be differentiated from the rest of the groups by having a deeper body. Within the upper Usumacinta group, the separation from specimens of the Miramar Lagoon was notable (Fig. 1 and Table 1). The latter had shorter anal fins and their pelvic fins were further back along the body. Statistical tests confirmed that the upper Usumacinta group was the most different; however, specimens from the Grijalva tended to exhibit slender bodies.

In both species, the greatest intraspecific morphological variation was principally observed in body shape, head size, and fin position. For many diverse species of fish, it has been demonstrated that variation in these anatomic attributes has functional importance and has been correlated with environmental factors such as water current speed, habitat structure, and the presence of predators (Langerhans 2008; Langerhans et al. 2003; Bravi et al. 2013; Araújo et al. 2014; Peris-Tamayo et al. 2020). This has also been documented in certain fishes of Middle America, such as characins (Santos and Araújo 2015; Garita-Alvarado et al. 2018, 2021), cichlids (Feilich 2016; Barrientos-Villalobos et al. 2018; Gómez-González et al. 2018; Aguilar-Contreras et al. 2021), and poeciliids (Araújo et al. 2014; Jourdan et al. 2016). In many fish taxa, phenotypic plasticity is an important precursor to morphological differentiation through adaptation to changing environmental conditions (Oufiero and Whitlow 2016).

Nevertheless, despite morphological evidence indicating that ecological-environmental factors may be promoting phenotypic differentiation in both species between the Grijalva and upper Usumacinta groups, the effect of geographic isolation and distance should also be considered—particularly for specimens of the Grijalva, which are the most geographically isolated (Sánchez et al. 2019). Furthermore, the specimens analyzed come from dams on the principal course of the Grijalva River. Notably, there is evidence of environmental changes and reductions in gene flow being caused by dam construction, which may induce rapid morphological modifications in fish (Michel et al. 2008; Franssen 2011; Perazzo et al. 2019; Gilbert et al. 2020). Additionally, while the Grijalva and Usumacinta basins share some biological and ecological characteristics, each basin has an independent biogeographical history (Elías et al. 2021). Likewise, they have been subjected to different modifications and anthropic activities (Lázaro-Vázquez et al. 2018; Sánchez et al. 2019). Also, it should be noted that hydrogeomorphological and habitat changes can affect the migratory behavior of these species—especially in the lower regions of both basins (Bloom and Egan 2018).

The morphological differentiation among D. petenense specimens in the upper Usumacinta appears to corroborate part of the hydrological history of the basin (Rosen 1967, 1970, 1979; Elías et al. 2021). Although this geographic group was the most highly differentiated within the study area, a notable separation among specimens from the Miramar Lagoon was evident in this group. This morphological evidence supports the hypothesis that diverse lineages exist within D. petenense in the Usumacinta basin (Elías et al. 2021). As presented in other biogeographic studies of the fish of northern Middle America, the territory that includes the upper Usumacinta has a particular independent evolutionary history characterized by geological events that promoted the geographic isolation of fish populations (Elías et al. 2021), such as the reversal of river current, collapse of stream beds, and underground connectivity among rivers (Rosen 1967, 1970, 1979).

Although morphological variation in the same anatomical attributes was found in both species, D. anale is less variable than D. petenense, mainly observed in the upper Usumacinta populations. Contrary to our expectations, despite being closely related and having similar ecological requirements, the magnitude and direction of morphological changes were distinct. This has also been found for other groups of fish for which the level of morphological variation among species is related to the level of dietary specialization (Ornelas-García et al. 2018; Kentao and Jearranaiprepame 2021). Regarding cichlids, generalist species were found to have greater variation in body shape than specialists (Kassam et al. 2004; Gilbert et al. 2020). It is also known that changes in head size and eye position may be related to dietary conduct, especially in terms of size and diet type (Haas et al. 2010; Kentao and Jearranaiprepame 2021).

Differences in the patterns of variation found in D. anale and D. petenense once again demonstrated that the phenotypic expression of morphological characteristics is a product of the interaction among diverse biological, environmental, and historic processes (Michel et al. 2008; Franssen et al. 2013; Bracciali et al. 2016). Thus, there is a need to continue conducting studies of the fish communities within the Grijalva and Usumacinta basins. As previously demonstrated, northern Middle America is an extraordinary natural laboratory for evaluating and understanding past and recent processes that have promoted morphological diversification among diverse species of Neotropical fish.

There is little intraspecific variation in body shape in D. anale and D. petenense. However, we uncovered statistically significant differences between specimens of both species from the Grijalva basins and the upper Usumacinta. Morphological differentiation was based on body height, head length, pelvic fin position, and anal fin length. Nevertheless, even though variation was observed for the same attributes in both species, the direction and magnitude differed. Since the variation in these morphological attributes seems to be related to biological, environmental, and geographic factors, it could serve to define ecotypes. For both species, morphological differences among specimens from the Grijalva basin could be due to geographic isolation. Meanwhile, differentiation among D. petenense specimens from the upper Usumacinta appears to support the hypothesis regarding the existence of two lineages in the Usumacinta basin. Additionally, in D. petenense, differentiation was detected among specimens from the Miramar Lagoon. Notably, there is a need for further taxonomic and biogeographic studies of the ichthyofauna in northern Middle America to better comprehend their diversity and the processes related to their evolution—particularly in the Grijalva and Usumacinta basins, which possess some of the most interesting and complex fish communities of the Neotropics.