Introduction

The Sonoyta pupfish (Cyprinodon eremus) (^{Miller and Fuiman, 1987}) is distributed in the Sonoyta River basin and the Quitobaquito Springs in northwestern México and the southwestern United States of America (USA) (^{Echelle et al., 2000}; ^{Miller et al., 2009}; ^{Minckley and Marsh, 2009}). The Quitobaquito Springs population is stable, however, the Sonoyta River population is decreasing. These declines are the result of exotic fish introductions, drought, and groundwater pumping that dramatically affects native aquatic habitats (^{Miller and Fuiman, 1987}; ^{Miller et al., 2009}; ^{Minckley and Marsh, 2009}). Currently, this trend of degradation and desiccation continues in the Sonoyta River promoted by the need for water for human growth (^{Miller et al., 2009}; ^{Minckley and Marsh, 2009}; ^{Minckley et al., 2013}). The only perennial water flow in the river that persists, is about 1 km in length at the Agua Dulce or Papalote locality (USON 0222, Table 1), and is maintained during the dry season by underground flow of shallow waters (^{Minckley et al., 2013}).

C. eremus is considered as endangered by the International Union for Conservation of Nature and the United States government (^{NatureServe et al., 2019}). Strategies to manage Cyprinodon spp. include protecting their habitat and developing refuge populations to increase or re-establish wild populations in cases of extirpation (^{Minckley et al., 1991}; ^{Minckley, 1995}; ^{Koike et al., 2008}). Refuges for C. eremus have been established in México to hamper the species’ gradual extinction (^{Minckley et al., 2013}). These conservation efforts also echo an initiative that sought to recover native fish in the arid southwest during the 1960s (^{Minckley, 1995}). In 1988, the first refuge population was created from the Sonoyta River populations in an artificial pond with only lentic habitat located at Centro Ecológico de Sonora (CES) in Hermosillo, Sonora (^{Marsh and Sada, 1993}). Between 2007 and 2011, another five refuges were established: one in the Biological Station and another in Schuk Toak Visitor Center (translocated to the Biological Station refuge), both of them at El Pinacate y Gran Desierto de Altar Biosphere Reserve (RBEPGDA), one at Centro Intercultural de Estudios de Desiertos y Océanos (CEDO), one at Colegio de Bachilleres in Sonoyta (COBACH), and another at Quitovac that include fish stocked in the springs (^{Minckley et al., 2013}). Currently, only the refuges at the RBEPGDA, COBACH, CEDO, and CES remain. The refuge in CES is the most populated (>1000 fish) for C. eremus in México. All the refuges were founded in México without an evaluation of their genetic variability.

Previously published papers on the pupfish family (Cyprinodontidae) have reported that populations with 15-30 years of isolation in distinct habitats develop morphological and genetic variations in the translocated populations (^{Collyer et al., 2005}; ^{Wilcox and Martin, 2006}; ^{Collyer et al., 2007}; ^{Koike et al., 2008}; ^{Lema, 2008}; ^{Collyer et al., 2011}; ^{Finger et al., 2013}; ^{Collyer et al., 2015}; ^{Black et al., 2017}). These morphological changes may be attributed to phenotypic plasticity in response to alterations in the environmental conditions of the new habitats (^{Collyer et al., 2005}; ^{Wilcox and Martin, 2006}; ^{Collyer et al., 2007}; ^{Lema 2008}; ^{Collyer et al., 2015}; ^{Black et al., 2017}), including evolutionary processes that occur on an ecological time scale (^{Collyer et al., 2007}; ^{Collyer et al., 2011}). Translocating individuals to be used as reproductive stocks poses a risk for species management due to phenotypic differentiation (^{Collyer et al., 2005}; ^{Wilcox and Martin, 2006}; ^{Collyer et al., 2007}; ^{Lema, 2008}; ^{Collyer et al., 2011}) and genetic adaptation to captivity (^{Frankham, 2008}), potentially reducing the survivability of the captive population during reintroduction in a wild environment.

Considering that the Sonoyta pupfish population of the Sonoyta River has decreased drastically and the habitat historical water flow does not exist, the CES refuge population represents an opportunity for the recovery of the species. However, it has been isolated since 1988 in an artificial pond under distinct environmental conditions compared with the wild populations. In this regard, we performed multivariate morphometric analysis (^{Blackith and Reyment, 1971}; ^{Reyment, 1982}) to characterize the morphological discrepancies between Sonoyta River wild C. eremus populations and CES refuge collected and founded from the same sample, respectively, with the goal of evidence the morphological variation after 29 years of isolation. The results obtained will help in redesigning, increasing the shelter area, and creating future management plans for species conservation.

Material and methods

Morphometric analysis

Based on ^{Hubbs and Lagler (2004)} and the box-truss protocol (^{Strauss and Bookstein, 1982}; ^{Bookstein et al., 1985}; Table 2; Figure 1), 35 morphological distances were measured considering that these characters underwent variation in both sexes of the genus Cyprinodon, as a result of the distinct habitats (^{Humphries et al., 1981}; ^{Collyer et al., 2005}; ²⁰¹⁵; ^{Black et al., 2017}). Seven meristic characters were also counted (Table 2) based on the description of C. eremus by ^{Miller and Fuiman (1987)}. Females and males were separately analyzed due to the sexual dimorphism of cyprinodontids. Each specimen was examined with a digital caliper (precision 0.01 mm) connected to a personal computer.

Table 2 Morphological distances modified from ^{Humphries et al. (1981)}, ^{Miller and Fuiman (1987)}, ^{Hubbs and Lagler (2004)}, ^{Collyer et al. (2005)}, ^{Collyer et al. (2015)}, ^{Black et al. (2017)}, and additional measures based on the box truss protocol (^{Strauss and Bookstein, 1982}; ^{Bookstein et al., 1985}) and quantified meristic characters for Cyprinodon eremus based on ^{Miller and Fuiman (1987)}.  

Tabla 2 Distancias morfológicas modificadas de ^{Humphries et al. (1981}), ^{Miller and Fuiman (1987)}, ^{Hubbs and Lagler (2004)}, ^{Collyer et al. (2005)}, ^{Collyer et al. (2015)}, ^{Black et al. (2017)} y medidas adicionales basadas en el protocolo box truss (^{Strauss and Bookstein, 1982}; ^{Bookstein et al., 1985}) y caracteres merísticos cuantificados para Cyprinodon eremus basados en ^{Miller y Fuiman (1987)}. 

Code	Morphometric character
M1-3	Dorsal length of head
M1-23	Upper lip - Center of the eye
M1-17	Head length
M1-12	Preanal length
M1-14	Prepelvic length
M1-16	Ventral length of head
M1-4	Length of upper jaw
M2-23	Lower lip - Center of the eye
M2-10	Standard length
M3-5	Occiput - Dorsal fin origin
M3-14	Occiput - Pelvic fin origin
M3-16	Occiput - Isthmus
M5-7	Length of depressed dorsal fin
M5-6	Length of dorsal fin base
M5-11	Dorsal fin origin - Base of the last anal fin ray
M5-14	Body depth
M5-16	Dorsal fin origin - Isthmus
M6-8	Dorsal length of caudal peduncle
M6-9	Base of the last dorsal fin ray - Ventral base of the caudal fin
M6-11	Anterior depth of caudal peduncle
M6-14	Base of the last dorsal fin ray - Pelvic fin origin
M8-9	Depth of caudal peduncle
M8-11	Dorsal base of the caudal fin - Base of the last anal fin ray
M9-11	Ventral length of caudal peduncle
M11-14	Base of the last anal fin ray - Pelvic fin origin
M12-13	Length of depressed anal fin
M14-15	Length of pelvic fin
M14-16	Pelvic fin origin - Isthmus
M18-19	Length of pectoral fin
M18-20	Length of pectoral fin base
M21-22	Eye diameter
A1	Interorbital width
A2	Head width
A3	Width of gape
A4	Body width
No.	Meristic character
1	Dorsal fin rays
2	Anal fin rays
3	Pectoral fin rays
4	Pelvic fin rays
5	Caudal fin rays
6	Scales from the dorsal fin origin to anal fin origin
7	Caudal peduncle scale count

Fig. 1 Landmarks for Box Truss protocol used in Cyprinodon eremus analysis. Black dots represent landmarks for distance measurements and the white dots represent the reference marks for the width measurements (Table 2 contains explanation of measurement codes).  

Fig. 1 Marcas para el protocolo Box Truss utilizado en el análisis de Cyprinodon eremus. Los puntos negros representan puntos de referencia para las medidas de distancia y los puntos blancos representan las marcas de referencia para las medidas de ancho (Tabla 2 contiene explicación de los códigos de medición). 

The regression described by ^{Elliot et al. (1995)} was performed to standardize the body measurements of each specimen. This regression model removes the component of size from the shape measurements (allometry), thereby homogenizing their variances (^{Jolicoeur, 1963}). Each character was standardized using the equation Ms = Mo (Ls/Lt)^b, where Ms = standardized measurement of the character; Mo = original measurement of the character (mm); Ls = average standard length (mm) of all the specimens from all the examined taxa; Lt = standard length (mm) of the specimen; and, “b” was estimated for each character from the observed data via the nonlinear regression equation M = aL^b. The parameter “b” was estimated as the slope of the regression log Mo on log Lt, using data from all specimens. The parameter “a” is the non-standardized measurement of the character (mm) and L = Ls/Lt.

The standardized morphometric data and meristic values of all the specimens were used to perform discriminant function analysis (DFA) and principal component analysis (PCA), which are the most used analyzes in multivariate morphometrics (^{Humphries et al., 1981}; ^{Reyment, 1982}; ^{Turan, 1999}). In the case of the PCA it does not require an a priori assignment of individuals into groups, but rather summarizes in linear combinations, called Principal Components, the variables that describe the shape variation in the combined sample (^{Humphries et al., 1981}; ^{Turan, 1999}). On the DFA, the individuals are assigned a priori into groups to calculate the function that better discriminates between the groups (^{Humphries et al., 1981}; ^{Turan, 1999}). The DFA was performed independently for females and males via a forward stepwise form using Statistica 5.0 software (StatSoft, Inc., Tulsa). It was performed to determine the combination of variables that optimally discriminated between wild and refuge populations. Statistically significant differences between populations were determined using Wilks’ lambda (λ), which oscillates from 0.0 (perfect discrimination power) to 1.0 (absence of discrimination). Values with p < 0.05 obtained in the DFA were considered statistically significant. The PCA was performed using the “factoextra” (^{Kassambara and Mundt, 2020}) and “FactoMineR” (^{Lê et al., 2008}) R packages (^{R Core Team, 2021}) to determine which morphological variables best explained the variability in the dataset.

The most important morphological characters selected by DFA and PCA were illustrated by violin and box plots. A one-way analysis of variance with a 95 % confidence interval was performed for each character to evaluate the null hypothesis of equality between the populations. After the null hypothesis was rejected, a post-hoc Tukey test was performed on Statistica 5.0 software (StatSoft, Inc., Tulsa) to verify whether the groups significantly differed (^{Turan, 1999}).

Results

DFA was performed on 120 specimens of C. eremus, females and males, from wild and refuge populations. 17 of the 41 morphological and meristic characters among females significantly distinguished the two populations (Table 3). The overall Wilks’ lambda (λ) was 0.06404 (p < 0.0001), indicating a high degree of discrimination between the two female populations. A significant difference was observed for eight variables (Table 3). According to PCA for wild and refuge females, principal components 1 and 2 combined explained 38.656% of the total variance, the PC1 and PC2 explained 25.329 % and 13.327 %, respectively (Supplementary Table 1).

The scatterplot shows segregation between wild and refuge females, mostly along PC1 (Figure 2A). The variables that most contributed to PC1 were body depth, base of the last dorsal fin ray to pelvic fin origin, dorsal fin origin to base of the last anal fin ray, dorsal fin origin to isthmus, body width, anterior depth of caudal peduncle, depth of caudal peduncle, and occiput to pelvic fin origin, among others (Figure 2B; Supplementary Table 2). Of these, only body depth was selected by the DFA, and it was slightly non-significant (p = 0.0565, Table 3). In addition, the DFA selected the length of dorsal fin base, occiput to isthmus, length of pectoral fin base, length of the upper jaw, and width of gape to significantly (p < 0.05) discriminate between groups (Table 3); these variables also contributed to PC1 in the PCA (Figure 2B; Supplementary Table 2). Conversely, head width, caudal fin rays, and pectoral fin rays were significant (p < 0.05) in the discriminant function (Table 3) but contributed least to PC1 in the PCA (Supplementary Table 2).

Fig. 2 PCA for the two C. eremus females’ populations: (A) Scatterplots showing the position of the females along the first two PCs, the ellipses represent the 0.95 confidence intervals; (B) the correlation circle of the 15 variables that most contribute to these PCs (See supplementary table 2 for more information about the contribution of the variables).  

Fig. 2 ACP para las dos poblaciones de hembras de C. eremus: (A) Gráfico de dispersión mostrando la posición de las hembras en los dos primeros CPs, las elipses representan los intervalos de confianza de 0.95; (B) círculo de correlación de las 15 variables que más contribuyen a estos CPs (Ver tabla suplementaria 2 para más información acerca de la contribución de las variables). 

The DFA selected 25 morphological and meristic variables that best discriminated the two male populations (Table 4). As observed in females, a high degree of discrimination between both populations was observed based on the Wilks’ lambda (λ = 0.08582, p < 0.0001). Significant differences were observed in 13 variables (p < 0.05) (Table 4). The first two PCs in the PCA explained 31.87 % of the male variability, where PC1 explained 19.46 % and PC2 12.41% of the variance (Supplementary Table 3). The variables that contributed more to the PC1 were body depth, dorsal fin origin to the base of the last anal fin ray, occiput to isthmus, head length, depth of caudal peduncle, dorsal length of head, and occiput to pelvic fin origin (Figure 3B; Supplementary Table 4). However, the two C. eremus male populations were mostly differentiated along PC2 (Figure 3A), where the variables that most contributed to this PC were head width, the width of gape, lower lip to the center of the eye, base of the last dorsal fin ray to pelvic fin origin, base of the last anal fin ray to pelvic fin origin, eye diameter, length of upper jaw, preanal length, and anterior depth of caudal peduncle, among others (Figure 3B; Supplementary Table 4). Of these variables, the head width, base of the last dorsal fin ray to pelvic fin origin, lower lip to center of the eye, length of upper jaw and anterior depth of caudal peduncle were also significant (p < 0.05) in the DFA.

Fig. 3 PCA for the two C. eremus males’ populations: (A) Scatterplots showing the position of the males along the first two PCs, the ellipses represent the 0.95 confidence intervals; (B) the correlation circle of the 15 variables that most contribute to these PCs (See supplementary table 4 for more information about the contribution of the variables).  

Fig 3 ACP para las dos poblaciones de machos de C. eremus: (A) Gráfico de dispersión mostrando la posición de los machos en los dos primeros CPs, las elipses representan los intervalos de confianza de 0.95; (B) círculo de correlación de las 15 variables que más contribuyen a estos CPs (Ver tabla suplementaria 4 para más información acerca de la contribución de las variables). 

Other key variables in both the DFA (p < 0.05) and PC2 of the PCA were ventral length of head, head length, length of depressed dorsal fin and length of depressed anal fin (Table 4; Supplementary Table 4). However, the body depth, base of the last dorsal fin ray to ventral base of the caudal fin, scales from the dorsal fin origin to anal fin origin, and dorsal fin rays were important in the DFA (p < 0.05) (Table 4) but contributed less to PC2 of the PCA (Supplementary Table 4).

Finally, the most significant morphometric variables in the DFA that contributed most to the PCA in both female and male populations, were plotted using violin and box plots (Figure 4). There were 19 divergent morphometric characters between wild and refuge C. eremus populations. Both sexes of the refuge population exhibited a greater length of the upper jaw but a shorter body width, length of the depressed dorsal fin, dorsal fin origin to the base of the last anal fin ray, anterior depth of the caudal peduncle, dorsal fin origin to the isthmus, and base of the last dorsal fin ray to pelvic fin origin compared with the wild population (Figure 4A-G ).

Fig 4 Violin and box plots of the most notable morphological characters for differentiating wild and refuge populations of Cyprinodon eremus. The letters above the bars represent the significant differences among groups according to the Tukey test (p < 0.05), and the Y-axis values are in millimeters.  

Fig 4 Gráficos de violín y diagramas de cajas de los caracteres morfológicos más destacados para diferenciar poblaciones silvestres y de refugio de Cyprinodon eremus. Las letras sobre las barras representan las diferencias significativas entre grupos de acuerdo con la prueba de Tukey (p < 0.05) y los valores en el eje Y están dados en milímetros. 

The females of the refuge population had a reduced width of gape, lower lengths from the occiput to pelvic fin origin and from the occiput to the isthmus, lower depths of body and caudal peduncle, and shorter length of the dorsal and pectoral fin bases compared with those of the wild C. eremus population (Figure 4H-N, respectively). Conversely, the males of the refuge population showed greater widths of the gape (Figure 4H H) and of the head, higher ventral length of head, and preanal length, but a shorter length from the base of the last anal fin ray to pelvic fin origin and shorter length of depressed anal fin compared with those of the wild C. eremus population (Figure 4O-S).

Discussion

The present study showed evidence of morphological variations in the refuge population of C. eremus, after 29 years of isolation in an artificial pond with a homogeneous environment distinct from its wild habitat. Upon comparing the CES refuge population with the wild population originally collected in a natural stream habitat from the Sonoyta River during the establishment of the refuge, we detected morphotypes associated with habitat type. Changes were observed in the mouth, head, body, and caudal peduncle regions. The refuge population had a longer upper jaw and varied width of gape. ^{McGee et al. (2013)} suggested that jaw traits affect feeding kinematics in fishes. Changes in head orientation and upward repositioning have been detected in Devil’s Hole pupfish Cyprinodon diabolis (^{Wilcox and Martin, 2006}) and Cyprinodon bovinus (^{Black et al., 2017}), which may be related to foraging behavior (^{Black et al., 2017}). Furthermore, changes in body depth and width were observed in the refuge C. eremus population. Similar variations observed in Cyprinodon pecosensis have been related to the size of their intestine due to the different types of food available in their habitat (^{Collyer et al., 2015}). The distributions and types of food available in the water column likely differ between the wild and refuge habitats of C. eremus. Ultimately, these variables may contribute to morphological changes observed in the wild and refuge populations.

The refuge C. eremus males had wider heads. Previously, C. pecosensis found in lentic sinkhole populations exhibited longer heads, which was attributed to a larger gill size adapted to prevent hypoxia in a low dissolved oxygen environment (^{Collyer et al., 2015}). Here, the C. eremus refuge population inhabits an artificial pond with limited water circulation and the presence of algae. These factors may contribute to a wider head, that allows more gill space, thereby reducing the risk of hypoxia when the dissolved oxygen levels in the refuge pond decrease.

Modifications in the pectoral fin attachment have been associated with enhanced maneuverability in the water column (^{Black et al., 2017}). Moreover, the reduction in the length of the pectoral fin base in C. eremus refuge females, depressed dorsal fin in refuge individuals, depressed anal fin in refuge males, and base of the dorsal fin in refuge females may be associated with a lower requirement for stability in the artificial pond, an environment without running water. Also, the anterior depth of the caudal peduncle was lower in both females and males from the refuge population; in the refuge females, the length from the base of the last dorsal fin ray to pelvic fin origin, was shorter, and the depth of the caudal peduncle was lower accounting for slender caudal regions. Variations in the caudal region of Cyprinodon have been associated with water flow and the presence of predators (^{Tobler and Carson, 2010}; ^{Collyer et al., 2015}). Thus, individuals from the refuge population, especially females, had slender caudal regions, potentially because they did not need to swim against the watercourse or move between habitats in the lentic pond environment.

Although morphological differences were observed between the wild and refuge C. eremus males, more disparities were observed between females between these populations. In C. diabolis, wild males have been shown to be more aggressive than two populations stocked in two artificial ponds in defending their respective territories (^{Wilcox and Martin, 2006}). Similarly, C. eremus males fight each other to protect their territory and reproduce with receptive females (^{Cox, 1966}). Unlike in C. diabolis males, within the mechanisms operating in wild C. eremus males for intimidating opponents and courting females, the male body shape plays an essential role, which is retained in the refuge males. Refuge males that maintained a body shape similar to wild males likely exhibited better fitness if the selective pressures of the environment were not strong enough to determine survival. Similar results were found in Cyprinodon tularosa, wherein males showed a positive association between body depth and size, likely related to the territorial defense, and females showed a decreased association between body depth and size (^{Collyer et al., 2005}). Therefore, the morphological variations in C. eremus females were more pronounced than in males, because morphological character selection in females may be regulated by environmental conditions and not by sexual selection.

^{Rodríguez-Ramírez et al. (2023)} recently studied the genetic variability of the CES and other two refuge populations and wild C. eremus from the Sonoyta River using seven microsatellite loci. The CES population showed less genetic variability compared to the others. This lower genetic variation in CES refuge is more related to the time of isolation in contrast to the others analyzed (^{Rodríguez-Ramírez et al., 2023}). As mentioned by ^{Koike et al. (2008)} and ^{Finger et al. (2013)}, typical long established pupfish refuge populations showed low diversity and significant divergence in allele frequencies.

Notably, lower genetic variability may reduce the survivability of the CES refuge C. eremus population and its ability to reproduce in its native environmental, as observed for other Cyprinodon spp. (^{Wilcox and Martin, 2006}; ^{Collyer et al., 2011}). This could hinder attempts to re-establish or increase native pupfish populations in the Sonoyta River, as has been the case for other Cyprinodontidae species (^{Hendrickson and Brooks, 1991}; ^{Black et al., 2017}). Therefore, it is necessary to evaluate the phenotypic and genetic diversities of the remaining wild populations and the rest of the refuges, including CES, to detect isolation-induced morphological and genetic variations. Consequently, a management plan is necessary for the conservation of C. eremus, considering the information on morphology and genetic variation. Measures must be taken to avoid morphological changes by increasing the heterogeneity of artificial habitats conditions rendering them more similar to wild habitats conditions (^{Black et al., 2017}). Wild individuals should also be translocated to the refuge population to increase genetic variability and therefore reduce the decline in fitness (^{Wilcox and Martin, 2006}; ^{Araki et al., 2007}; ^{Frankham, 2008}; ^{Black et al., 2017}; ^{Rodríguez-Ramírez et al., 2023}). We recommend an extensive survey to obtain samples from all localities and habitat types of the entire C. eremus distribution to account for the phenotypic variability in each habitat.