LC and LZ are located in Central Mexico, within the Lerma-Chapala basin (Fig. 1). According to a morphostructural geological analysis, both lakes were jointed at the Pliocene-Pleistocene boundary by wide channels related to morphotectonic structures and then separated during the Pleistocene-Holocene period [36]. LC, the second largest natural lake in the country, is a shallow tropical system (375 km${}^2$) located between 20º04’34” and 19º53’15”N and 101º46’45” and 101º47’25” W, with one main tributary (Grande de Morelia River) and a manmade effluent (Yuriria Channel) [37]. In an east-west direction, this lake is divided into three zones according to the horizontal physical, chemical, and biological patchiness; it has frequent volume fluctuations and the western zone dries up in years of severe droughts [37]. Several impacts affect this lake, including watershed degradation by deforestation and soil erosion, water withdrawal from its main tributary, water pollution from the state capital city, industry, and agricultural activities. As for the most common threats within the ecosystems are the introduction of non-native species, increasing fishers and nets, organizational capacity loss, and overfishing [37].

LZ is a small remnant (0.335 km${}^2$) of a larger wetland (261 km${}^2$) located between 19º49’26” and 19º49’40” N and 101º46’45” and 101º47’25” W. The main tributaries to the lake are two groups of springs (La Angostura and La Zarcita), and the natural effluent is the Angulo River [33]. This small spring-fed lake is threatened by the proximity to the urban area, the agricultural use of the southern shoreline, the water subtraction from the springs that feed the lake, and the introduction of exotic species [38].

To relate species densities to environmental variables, a data set of H. turneri in LC was used, including 448 individuals collected from five sites in four months during 1979 (February, April, August, and November). This is the last report of the species in the area after comprehensive ichthyological surveys [13]. The environmental data available for this lake consisted of annual average values at each site (Table 1). The data were collected at the surface during the day simultaneously with the fishes’ catches. In LZ, 2,135 individuals were caught and analyzed in four months of 1995 (January, May, July, and October). For this lake, water quality information was also recorded simultaneously with fish capture at the surface at four sites during the day and nighttime (Table 1). Also for LZ, environmental data and fish samples were collected during 2019 and 2020 at the same sites and months during the day (274 individuals). The sites in both lakes were selected in such a way to represent the different habitat characteristics (water inflow from springs and rivers and the effluents, marshy areas and zones with distinct submerged macrophytes species, and nearshore influence from agricultural activities and urban development). Fishes were collected in 1979 and 1995 using a 50 m long trawl net with a 4 m $\times$ 1.8 m rectangular mouth and 5 mm $\times$ 5 mm mesh wings with a cod end; the sample area at each site was approximately 165 m${}^2$. The fish density was adjusted for the net area, and the results were expressed in the number of fishes per square meter. In 2019 and 2020 we used a smaller long trawl net of 25 m length with a 4 m $\times$ 1.8 m rectangular mouth and 5 mm $\times$ 5 mm mesh wings with a cod end. Due to the small number of individuals captured, we also used five cylindrical minnow traps set for one hour per site (stainless steel, 42 cm long and 19 cm in diameter, stretch mesh 0.5 cm, with two 2.5 cm holes with inverted cone inlets) and electrofishing (Power ~200 W, peak voltage ~250 V, peak current ~10 amp). All sites in both lakes were located nearshore due to their accessibility and according to the fishing gears used for sampling.

The population density and frequency of H. turneri were analyzed and compared with the information of other fish species, aiming at describing the status of the species and identifying potential interactions, mainly with the introduced organisms in LC. We used Spearman’s correlation coefficients to measure the strength of the spatial association between species. In addition, we implemented a cluster analysis, with the sites as the grouping units, to recognize the species distribution affinities in multidimensional space. We applied the Bray-Curtis dissimilarity measure and Ward’s minimum variance method to link similar points [39]. The correlation coefficients and significant values were calculated with the ‘Hmisc’ package (v. 4.5-0) [40], and clusters were obtained with the ‘vegan’ package (v. 2.5-7) [41]. Both methods were performed in the R statistical language (v. 4.1.2, Vienna, Austria) [42].

The relationships between environmental aspects and H. turneri abundance were explored with a generalized additive mixed model (GAMM) in 1995 in LZ. Environmental variables were analyzed to search for outliers and screened for collinearity via Cleveland dotplots and by applying Pearson’s correlation coefficients ($>$0.75), respectively [43]. Bacteria’s counts were log-transformed, and nitrites square root transformed to avoid large values. Correlations were found between temperature and DO (r${}^{\mathrm{2}}$ = –0.76, p 0.001) and temperature and conductivity (r${}^{\mathrm{2}}$ = –0.76, p 0.001), as well as between nitrite and hardness (r${}^{\mathrm{2}}$ = 0.89, p 0.001). We removed different covariates according to collinearity (e.g., temperature and nitrite), low average spatial variation within optimal values for H. turneri (e.g., pH, conductivity), and to avoid skewness because several parameters represent similar water properties (e.g., turbidity, total bacteria).

A GAMM was applied, assuming a Poisson distribution because species counts include a low proportion of zeros (9.4%). Five covariates were used as smoothers (spline functions), assuming no strong relationships between covariates and response variables. We used random effects to model the correlation in space and time [43] because spatial and temporal autocorrelation was expected, since distance among sites is small (average 375.5 m) and samples were obtained every three months. We included a random intercept of Month (impose correlations between observations in the same sampling event) and Site (nested in Month) as the random intercept to allow for correlations between observations made in the same month (small scale spatial and temporal correlation). The final GAMM construction for H. turneri count data followed the equation:

Hubbsina${}_{i\mathit{}j\mathit{}k}$ ~ Poisson($\mu$${}_{i\mathit{}j\mathit{}k}$) log(Hubbsina${}_{i\mathit{}j\mathit{}k}$) = $\alpha$ + Total phosphorus${}_{i\mathit{}j\mathit{}k}$ + fDissolved oxygen${}_{i\mathit{}j\mathit{}k}$ + fDepth${}_{i\mathit{}j\mathit{}k}$ + fNH${}_{\mathrm{3}\mathit{}i\mathit{}j\mathit{}k}$ + fSuspended Solids${}_{i\mathit{}j\mathit{}k}$ + flogTotal coliforms${}_{i\mathit{}j\mathit{}k}$+ Month${}_i$ + Site${}_{i\mathit{}j}$

The model establishes that the counts of H. turneri at month i, site j, and observation k follow a Poisson distribution with mean $\mu$${}_{ijk}$. The log link function is implemented with the expectation that the values of H. turneri and $\mu$${}_{ijk}$ are a function of the covariates and avoid non-negative fitted values. We assessed the overdispersion of the model by dividing the sum of squared Pearson residuals by the difference of the number of sampling units minus the degree of freedom from each component in the model [43]. The GAMM was fitted using the gamm4 package (v. 0.2-6) [44] within the R language (v. 4.1.2, Vienna, Austria) [42].

The optimal linear estimation method proposed by Solow and Roberts [45] was applied to estimate a calendar date interval when H. turneri might be considered extirpated or locally extinct from LC with 95% statistical confidence. Collection records from several scientific collections were used to implement the model, including records from the University of Michigan Museum of Zoology, National Museum of Natural History-Smithsonian Institution, Tulane University, and the Colección Nacional de Peces Dulceacuícolas Mexicanos de la Escuela Nacional de Ciencias Biológicas, IPN (Table 2, Ref. [13]). This approach is based on the shape parameter of the Weibull distribution for the interval between successive dates. The model is described as:

(1)$$v=\frac{1}{k-1}\sum _{i=1}^{k-2}\log \frac{T_1-T_k}{T_1-T_{i+1}}$$

where T${}_{\mathrm{1}}$$\mathrm{>}$T${}_{\mathrm{2}}$$\mathrm{>}$…T${}_k$ are the k most recent sighting times of a species, ordered from most recent to least recent. The confidence interval is given by:

(2)$$T_1+\frac{T_1-T_k}{S_L-1},T_1+\frac{T_1-T_k}{S_U-1}$$

where the lower (S${}_L$) and upper (S${}_U$) limits are calculated as:

(3)$$S_L={\left(-\log \frac{1-\frac{\alpha }{2}}{k}\right)}^{-v}$$

(4)$$S_U={\left(-\log \frac{\frac{\alpha }{2}}{k}\right)}^{-v}$$

were confidence intervals were set at $\alpha$ = 0.05.

Hubbsina turneri had the sixth-highest abundance and cohabited with 12 species of fishes in LC, and the dominant species were the native atherinopsid (Chirostoma jordani) and goodeids (Goodea atripinnis, Xenotoca variata, and Zoogoneticus quitzeoensis). Three introduced species were reported, two cyprinids (Cyprinus carpio and Carassius auratus) and one cichlid (Oreochromis niloticus). The spatial distribution of Hubbsina turneri showed smaller values at sites related to the affluent and effluent (10 and 12 individuals at sites 4 and 5, respectively) and higher values at sites within the lake (64, 212, and 150 individuals at sites 1, 2 and 3, respectively). The sites with higher abundances were related to the presence of submerged macrophytes (Potamogeton filiformis and Potamogeton pectinatus). Temporally, more than 50% of the organisms (289 individuals) were captured in April, and only 3.6% (16 individuals) were captured in February.

In the fish community in LZ during 1995, H. turneri cohabited with 13 species, and the dominant species were the native cyprinid (Notropis grandis), goodeids (G. atripinnis and Skiffia lermae), and atherinopsid (Chirostoma humboldtianum). Hubbsina turneri had the fifth-highest abundance. As in LC, H. turneri exhibited smaller abundance values in the main affluent (214 at site 2), but with higher values in the more stagnant site and the effluent (876 and 748 individuals at sites 1 and 4, respectively). The occurrence of H. turneri was also related to the presence of submerged macrophytes (P. pectinatus and Ceratophylum demersum). Temporally, 35.7% of the organisms (763 individuals) were captured in October and only 17% (362 individuals) were captured in January. During 2019 and 2020, from the trawl net information, the dominant fish changed and C. humboldtianum and the goodeid X. variata had the higher values. Hubbsina turneri only occurred in 4% of the samplings and occupied the eighth position of the nine species collected. The few places where the species was captured were the same where it was dominant in 1995. With the electrofishing and the minnow traps, we captured more individuals of H. turneri (218), and at more sites, than with the trawl net (56).

In LC, H. turneri was highly correlated with Alloophorus robustus (r${}^{\mathrm{2}}$ = 1, p 0.001) and Poeciliopsis infans (r${}^{\mathrm{2}}$ = 0.97, p = 0.005), and the species was more related to C. auratus (r${}^{\mathrm{2}}$ = 0.8, p = 0.1) than to other introduced species (O. niloticus: r${}^{\mathrm{2}}$ = 0.7, p = 0.19 and C. carpio: r${}^{\mathrm{2}}$ = 0.37, p = 0.54). Three species groups were identified in the multidimensional space analysis: group 1 included the dominant fish, group 2 included the two introduced cyprinids and one native species, and group 3 included the intermediate dominant species including H. turneri, two other goodeids, and the introduced cichlid but the poecilid (Fig. 2). In 1995 Hubbsina turneri in LZ was more related to the goodeids Zoogoneticus quitzeoensis (r${}^{\mathrm{2}}$ = 0.74, p 0.001), A. robustus (r${}^{\mathrm{2}}$ = 1, p 0.001) and X. variata (r${}^{\mathrm{2}}$ = 1, p 0.001) than other species. About the introduced species H. turneri was more related to Ptenopharingodon idella, which was captured at site 2 (r${}^{\mathrm{2}}$ = 0.45, p = 0.01). Two groups were differentiated from the cluster analysis: (1) fishes with high densities including H. turneri and (2) species with low abundance values (Fig. 2). In 2019 and 2020, we found no relation of H. turneri with the other species because of the low individual number. Two species groups were identified in the multidimensional space analysis: (1) The intermediate and dominant species and (2) the species with lower abundance, where H. turneri is now located along with the predator A. robustus (Fig. 2).

Both lakes showed contrasting environmental characteristics. On the one hand, LC was shallow and saline with higher temperatures and higher values of different chemical parameters (conductivity, hardness, salinity, alkalinity, and pH) than LZ. On the other, LZ reached deeper areas, with very high dissolved oxygen (DO) values and low turbidity, conductivity, and hardness (Table 1). The correlation analysis revealed a close relationship between H. turneri with depth (r${}^{\mathrm{2}}$ = –0.9, p = 0.04) and DO (r${}^{\mathrm{2}}$ = 0.9, p = 0.04) in LC. In LZ, H. turneri was negatively correlated with total coliform bacterias (r${}^{\mathrm{2}}$ = –0.55, p = 0.03) and DO (r${}^{\mathrm{2}}$ = –0.49, p = 0.06).

When we checked the overdispersion value in the GAMM, the result was sufficiently close to 1 (1.14), meaning that the model had low overdispersion and an adjusted R-squared of 0.54. The GAMM outputs, and the H. turneri distribution according to the general water quality parameters in both lakes, indicated the characteristics of the sites where it would be more likely to find higher densities of the species (Table 3). Shallow to medium depth (up to 3–4 m); high oxygen concentrations, where LZ had oversaturation and the species was found in medium values (from 13 mg L${}^{-1}$), meanwhile in LC only in the higher values ($>$6.5 mg L${}^{-1}$); suspended solids (between 8.5 to 9 mg L${}^{-1}$); soft to moderate hardness water (35–40 in LZ, and 100–140 mg CaCO${}_3$ L${}^{-1}$ in LC); a pH from neutral to basic (7.0–8.7); alkalinity values lower than 30 mg CaCO${}_3$L${}^{-1}$. According to the nutrient and bacteria information obtained only in LZ, the species were mainly found in sites with higher NH${}_3$ (above 0.25 mg L${}^{-1}$) and total phosphorus values, and lower values of total coliform bacteria (Fig. 3).

The optimal linear estimation analysis indicated that the time interval within which H. turneri may be regarded as extirpated from LC is between the years 2013 and 2018; in other words, the species may have already been extirpated. After this time, in statistical terms, the species may be considered locally extinct.

The survival rate of a species is a historical indicator of the differential anthropogenic impacts on communities among aquatic systems. In goodeids, this value is 0.57 for LC and 0.94 for LZ [46]. A decline in the number of fish species in LC was described between 1970 and 1980 [13]. Of the 15 native species originally documented, only 4 occurred in the 1990s [13]. The fishery showed an important production in the lake at the beginning of the 1990s, followed by a drop in the second half of the decade related to a severe drought, which was then followed by a continuous decline in the 2000s (3250 t in 1994, 730 t in 1999, and 300 t in 2006). It is unclear what aspects contributed specifically to the disappearance of the species in the lake. Although the probable main causes are noted, we assume that the severity of the effects depends on the magnitude, frequency, duration, and combined impacts of the causes [47], which created circumstances that did not allow the population to recover, some of them are:

First, habitat reduction negatively affects population size and viability [48, 49]. There were significant surface fluctuations in LC due to hydrometeorological variation. The last report of the species in the lake coincided with highly dry conditions (1980–2000) [50]. Two-thirds of the lake remained dry in 1988, and the maximum depth was one meter [51]. Volume drop also relates to water use because the state capital city is located within the lake basin. The population almost doubled between the 1970 to 1980 decades [52], and the water demand increased likewise (18 $\times$ 10${}^3$ m${}^3$ in 1970 and 33 $\times$ 10${}^3$ m${}^3$ in 1980) [53]. The substantial exploitation and management of water resources (e.g., dams, water diversion structures, and over-extraction of aquifers) had caused severe ecological damage and had imposed real threats to imperiled species in basins everywhere [54, 55, 56]. Lake volume fluctuation was associated with mineralization changes, which in LC was reflected in alkalinity (lake average in mg CaCO3 L${}^{-1}$: 856 in 1980, 1407 in 1985, and 403 in 2002) and conductivity (lake average mS cm${}^{-1}$: 3 in 1980, 9 in 1985, and 1 in 2002) [57]. Ion concentration could exceed the organisms’ physiological adaptability, leading to the cessation of growth and increased mortality [58], an aspect that could squeeze H. turneri into marginal conditions. This pattern has been described temporally with environmental variables in other ecosystems [58].

Second, eutrophication is one of the most threatening effects in lakes because a drastic increase in nutrient concentration affects water quality and promotes changes in biodiversity and biogeochemical processes [11, 13]. The distribution area of H. turneri in LC (eastern zone) presented the highest Chlorophyll-a and nitrate concentrations associated with inputs from the Grande de Morelia River [59]. This river crosses the capital state city, where the treatment plant wasn’t constructed until 2007 [59]. In addition, the river flows and is channeled in the Morelia-Queréndaro irrigation district, which has a surface area of 200 km${}^2$ [60]. In the 1970s, a significant change related to the green revolution occurred in Mexican agriculture, when mechanization and agrochemicals use significantly increased. In this region, between 1970 and 1980, the cultivar structure switched from cereal grains for human consumption to forage species [61]. From the 1980s records, a consistently increasing trend in phosphates has been reported (lake average in mg L${}^{-1}$: 0.28 in 1980, 0.23 in 1985, and 72.4 in 2002) [57]. In addition, gully erosion increased from 9 km${}^2$ in 1975 to 23 km${}^2$ in 2000 [44].

Third, biological interactions with native as well as exotic species at different niche aspects, which impose direct impacts like competitive exclusion, niche displacement, predation, and indirect impacts, such as habitat alteration and the spread of emerging diseases, can be particularly severe in imperiled populations [62, 63]. Volume drop could reduce survival rates of vulnerable fish in lake zones due to limiting similarity processes like competition for food [64, 65]. Hubbsina turneri and Z. quitzeoensis are both species that inhabit similar habitats and consume analogous prey (chironomids). However, Z. quitzeoensis was more abundant (1.03 org m${}^{-2}$ and 0.53 org m${}^{-2}$, respectively) and consumed other common prey (the insect Notonecta). Predators influence prey’s population structure and dynamics; H. turneri is prey for A. robustus and may have been more vulnerable to predation if macrophyte-covered habitats no longer concealed it. This assumption has been corroborated in other predator-prey interactions, such as the higher success of largemouth bass (Micropterus salmoides salmoides) predation on bluegill (Lepomis macrochirus) in low macrophyte stem densities [66]. In the small, conserved lake La Mintzita, H. turneri disappears, and Z. quitzeoensis and A. robustus persist with higher abundances [67]. The interaction for space and food could also occur with introduced species [68]. In general, introduced species are strong competitors, and some studies have described the competitive interaction between O. niloticus and native fishes that results in the displacement from the structured habitats occupied by the native fishes [69].

Fourth, fisheries could contribute to the disappearance of populations and species in some distribution areas. The continuous decline of several species, including those captured as bycatch, is an indirect effect of using non-selective fishing gear, which is an aspect that is discussed in different environments [70, 71]. In LC, seine nets have been used to capture the more valuable and locally appreciated silverside Chirostoma jordani. Because of the mesh size of these nets, it is common to capture other juveniles and adults of native and introduced species.

Official records described a 150 km${}^2$ lacustrine zone in the Zacapu wetland with depths up to eight meters and several springs that were mainly along the south shore [72]. However, in 1896, the federal government approved the desiccation of approximately 123 km${}^2$, and the region turned into the main grain producer in the state [73]. This action fragmented the H. turneri populations, which are now restricted to LZ, and it was also present in the years 2006 and 2007 in the spring of the Naranja de Tapia town (12 km SE from LZ), and a small creek in the town of Jesus María (33.42 km NE from LZ) [33]. The water and habitat quality and biotic integrity evaluations in LZ fluctuate between regular and good, and the fish composition has remained the same [38]. However, the number of native species in the Naranja de Tapia spring has declined from five reported in 1995 to only one in 2006, and it presents poor biotic integrity; H. turneri and other species only persisted in a channel leading away from the system (Fig. 4) [38]. The situation of the species is also precarious in the Jesús María creek because it is located behind a livestock barn and the habitat is continually altered by cattle trampling (Fig. 4). Livestock grazing and trampling alter the complexity of the habitat in stream channels and the structure, pattern, and processes of riparian vegetation and, ultimately, the survival of other aquatic species like fishes [74].

A synthetic fiber factory was installed in 1948 on the shore of the LZ effluent [75]. This factory accelerated the urban growth, and the population of the area doubled between 1940 and 1950 (from 6 169 to 14 346 inhabitants). In 1995 the population grew 3.5 times more (48 307 inhabitants) and the last population census in 2020 reported 55 287 inhabitants. This growth exponentially increased the demand for basic services as well as the production of wastewater [75]. The city of Zacapu had a complete wastewater network and treatment plant with a treatment capacity of 120 L s${}^{-1}$ installed until 2000 [15]. However, there were some reports of wastewater discharges in LZ from the city meat processing facility in 2010 [76]. Most environmental variables remain similar at the different decades but NO${}_3$, which had an important increase in the recent samples (average $\pm$ sd in 1995: 1.1 $\pm$ 0.3 mg L${}^{-1}$; in 2019–2020: 7.3 $\pm$ 0.8 mg L${}^{-1}$). The natural level in the surface water of nitrates and other forms of nitrogen is typically less than 1 mg L${}^{-1}$, but the excess of nitrates is related to eutrophication processes [77]. The habitat structure in lakes is important for supporting biota and it is related to the presence of springs, riparian vegetation and littoral cover complexity, high diversity of underwater plant species, substrate diversity, and anthropogenic disturbance [78]. Lake Zacapu had a spatial heterogeneity general gradient in an east-west direction from low to high values [79], and something that has been perceived between decades, but not quantified, is a greater modification in the south shore due to the establishment of agricultural activity and removal of submerged vegetation in the west zone. This modification in the habitat structure may be related to the reduction in the population of some species, influencing the change in the dominance from littoral to more limnetic fish. Another imminent threat found in the most recent sampling on the springs of La Angostura, the main affluent of LZ, was the presence of the species Pseudoxiphophorus bimaculatus and Xiphophorus hellerii, because negative interactions with native fishes have been reported for both species in other aquatic ecosystems [80].

The current distribution of H. turneri is of special concern because fishes with restricted distributions are generally more prone to extinction than widely distributed species [81, 82]. To devise conservation and restoration measures for this species, it is important to understand the status of the lakes the species has been inhabiting. Lake Cuitzeo (as part of an extensive region) and LZ (by itself) are recognized as priority hydrological regions by the National Commission for the Knowledge and Use of Biodiversity [83]. Specifically, LZ was declared as a natural protected area in 2003 [84] and is a Ramsar-listed wetland (No. 1465, designation date 05-06-2005) [85].

The central issue regarding the reintroduction of H. turneri in LC is that the system requires a major intervention, including habitat restoration at multiple spatial scales [86]. Several plans have been elaborated by the government [84, 87] and research institutions [50, 88], and they have already answered some critical questions, including those related to the disturbances associated with land use, the amount of natural land cover in the basin, the impact and treatment of residual solids, and the implementation of wetlands to reduce nutrient inputs [1, 89]. Unfortunately, the misperception of the lake as a decaying, senescent, and hypereutrophic water body has overshadowed the understanding of the multiple ecosystem services it provides, which has reduced the support for restoration efforts [46].

Lake Zacapu is important to protect because it provides many ecosystem services, including urban, industrial, recreational, and agricultural water use services to the town and wetlands of Zacapu. The lake has an artisanal fishery where ten local families harvest silverside (Chirostoma humboldtianum) and common carp (Cyprinus carpio) and extract clams (Anodonta grandis grandis) and crustaceans (Cambarellus montezumae) [33]. In addition, LZ is an important habitat for endemic fishes (Notropis grandis and Allotoca zacapuensis) and amphibian species (Ambystoma andersoni) [90]. This lake represents a groundwater-dependent ecosystem that is influenced by the chemical, ecological and hydrological characteristics of infiltration and water transport, and it is threatened by different land-use activities, pollution, and climate change [91]. In this context, higher elevation areas and adjacent recharge areas must be protected, particularly at the level of landscape units related to volcanic processes with pine-oak forests [92].

We suggest the implementation of three actions to preserve the LZ ecosystem: (1) protect spring zones and regulate the water extraction for the conservation of aquatic biodiversity, (2) restoration of the littoral areas, and (3) adapt decision-making initiatives of the lake use and urban development to long-term objectives to preserve the ecosystem structure and function [5, 93]. Exotic species may represent a threat to native species, and non-native fish control measures must be implemented to reduce the already established populations (i.e., reducing restocking and increasing harvesting). New introductions must be avoided, and a well-implemented educational campaign is critical to prevent the release of organisms by aquarists into the lake [1]. Hubbsina turneri and goodeid fishes lack substantial economic fishery value, but the family has drawn attention from conservationists worldwide (i.e., Goodeid Working Group). Several populations are kept in captivity as a source of species data and as a gene bank by hobbyists, universities, public aquaria, and zoos. The data on this species include a detailed protocol on population treatment for restoration and reintroduction projects [94].