Introduction

Taxonomic distribution patterns arise from interactions among historical, ecological, and stochastic processes (^{Halffter & Morrone, 2017}; ^{Haydon et al., 1994}; ^{Wiens & Donoghue, 2004}). The analysis of these events within a temporal framework may explain how geological and climatic changes promote diversification patterns (^{Bryson et al., 2012}; ^{Martínez-Méndez et al., 2019}). The integration of both historical and ecological components is required to explain species richness patterns and community composition at local and regional scales (^{Holloway, 2003}; ^{Wiens & Donoghue, 2004}).

In historical biogeography, the shared patterns among organisms with different dispersal capabilities allow to evaluate the relationships between different geographic units through identifying generalized tracks (GTs), which arise from the geographic congruence of different individual tracks (ITs). GTs suggest the existence of ancestral biotic components that were affected by vicariant events (^{Morrone, 2001}). Furthermore, panbiogeographic nodes (PNs) that arise as a product of the convergence of individual or generalized tracks represent regions where different biogeographic histories overlap, indicating areas of high species richness (i.e., hotspots) that could lead to the establishment of conservation policies (^{Morrone, 2001}). Analyses of spatial diversity patterns at different scales (e.g., alpha diversity, or species richness within a community; gamma diversity, or species richness within a landscape; and beta diversity, or variation in species composition among communities within a landscape) are fundamental for understanding community patterns and processes (^{Hawkins et al., 2003}). Such analyses also link local community structure and composition to regional and continental assemblages (^{Cornell & Lawton, 1992}).

Studies in Mexico have evaluated diversity patterns using many different criteria (^{Espinosa-Organista et al., 2000}; ^{Morrone et al., 2017}; ^{Olson et al., 2001}). The existence of general distribution patterns can be determined from geological history, climatic changes, GTs, and PNs (^{Aguilar-Aguilar et al., 2008}; ^{Contreras-Medina et al., 2007}; ^{Corona et al., 2009}; ^{Escalante et al., 2018}; ^{Huidobro et al., 2006}). On the other hand, by using grid-cells in analyses of spatial diversity, researchers have recognized that most species occupy small areas; finally, other studies highlighted the role of the Mexican Transition Zone and the occurrence of diversity hotspots (^{Aguilar-Aguilar et al., 2008}; ^{Koleff et al., 2008}; ^{Ochoa-Ochoa et al., 2012}).

Globally, continental gastropods are one of the most diverse animal groups; their low vagility and dispersal potential, combined with their high endemicity, make them an ideal model for biogeographic studies (^{Hylander et al., 2005}). However, beyond a few regional inventories, no studies have addressed the diversity patterns of continental gastropods in Mexico (^{Correa-Sandoval et al., 2007}; ^{Naranjo-García, 2003}; ^{Naranjo-García & Fahy, 2010}; ^{Van Devender et al., 2012}). Diversity patterns in continental gastropods are different compared to other terrestrial taxa such as vertebrates, besides, the patterns differ between land and freshwater snails, in the case of land snails, most of the diversity is concentrated in islands and limestone mountain ranges (^{Solem, 1984}); on the other hand, freshwater snails richness is concentrated in springs and groundwaters, large rivers and their tributaries, ancient oligotrophic lakes and monsoonal wetlands (^{Strong et al., 2008}). In this study, a comprehensive biogeographic approach was conducted to identify GTs and PNs and evaluate spatial diversity patterns in Mexican continental gastropods as a whole; additionally, we evaluate land and freshwater snails regardless to evaluate if there are similarities in their distribution patterns despite their different environments.

Materials and methods

Records of continental gastropods in Mexico were gathered from a checklist by ^{Thompson (2011)}, an updated checklist of freshwater snails by ^{Czaja et al. (2020)}, and the GBIF database (http://www.gbif.org/). In addition, a search in the available literature was carried out and new species described from 2008 to the present were incorporated. Collection records with missing coordinates were georeferenced using maps and gazettes. Introduced species recognized by ^{Thompson (2011)} and new introduced species reported by ^{Naranjo-García and Castillo-Rodríguez (2017)}, as well as ambiguous and duplicated geographic records, were not considered. The taxonomic status of each species was verified using Molluscabase (http://www.molluscabase.org/) and original descriptions. The below mentioned analyses were performed considering 3 sets of data, the complete continental gastropods fauna and land snails and freshwater snails separately.

The majority of gastropod species in Mexico have restricted distributions, and most of them are only known from their type locality. Nevertheless, a panbiogeographic analysis involve the comparative study of several unrelated taxa and the integration of these with the tectonic history (^{Heads, 2005}). For this reason, panbiogeographic track analysis was performed at the genus level as has been used in other studies where the number of records diminish considerably if used at species level (^{Palma-Ramírez et al., 2014}). Individual tracks were obtained using minimum spanning trees constructed from the geographic distances among collection records and the Prim algorithm in the ‘fossil’ package (^{Vavrek, 2011}) in R (^{R Development Core Team, 2017}); this analysis was based on the methods proposed by ^{Escalante et al. (2017)}. Generalized tracks (GTs) were obtained using a parsimony analysis of endemicity with progressive character elimination (PAE- PCE; ^{Echeverry & Morrone, 2010}). The 14 biogeographic provinces defined by ^{Morrone et al. (2017)} were used as geographic units. ITs were codified as multistate non additive characters (absent = 0, present = 1, and crossing through the province = 2). Phylogenetic analyses were performed in TNT version 1.1 (^{Goloboff et al., 2008}) using a traditional search based on 1,000 replicates of random sequence addition and branch exchange via tree bisection and reconnection (TBR). Ten trees were retained in each run, and a final strict consensus tree was generated from the most parsimonious trees. GTs in the tree were identified as nodes supported by at least 2 geographic synapomorphies (ITs) shared by a biogeographic province. This procedure was iterated until no more GTs were recovered. Finally, GTs and panbiogeographic nodes (PNs) were plotted using the extension Trazos2004 (^{Rojas-Parra, 2007}) distinguishing between individual nodes (IN) that represent the intersection of at least 2 individual tracks and generalized nodes (GN) that represent the intersection of 2 generalized tracks (^{Morrone, 2015}).

The geographic area of Mexico was divided into 184 grid-cells of 1° latitude by 1° longitude to control for the effect of area size on diversity values (^{Koleff et al., 2008}). One-degree cells were used instead of the 30-second grid- cells used in other studies to avoid a large number of empty cells where no gastropod records are available. A presence-absence matrix (PAM) was constructed from the 184 grid-cells and the species recorded in this study. Alpha diversity and average proportional range size were generated (Koleff et al., 2008) and edited in QGIS version 3.14.0 (^{Team, 2016}), meanwhile, the Shannon index was calculated in DIVA-GIS version 7.5.0 (^{Hijmans et al., 2001}). Cells with proportion of species higher than 10% were recognized as diversity hotspots for the continental gastropods, meanwhile, for the land and freshwater snails datasets the threshold value was fixed at 7% due to lower number of species recorded. After it, the cells were intersected with the terrestrial ecoregions of the world and the freshwater ecoregions of the world to have a more precise location of the hotspots (^{Abell et al., 2008}; ^{Olson et al., 2001}).

Finally, to compare the spatial distribution of continental gastropods, as well as land and freshwater snails among the Nearctic and Neotropical regions and the Mexican Transition Zone (^{Morrone et al., 2017}), we recorded the number of total and endemic species in each region.

Results

A total of 4,379 records for 1,075 continental gastropod species were recovered, of them, 890 were terrestrial species and 185 were freshwater species. The species were from 187 genera in 63 families (Supplementary material Table S1, available at: https://github.com/HOmarMejiaG/diversity-patterns).

The panbiogeographic analysis for the continental gastropods included 159 ITs from 187 recorded genera; 28 of the genera were restricted to a single biogeographic province. Three PAE-PCE runs were performed until no more GTs were recovered. The first PAE-PCE run produced 1 tree with a length (L) of 786, a consistency index (CI) of 0.40, and a retention index (RI) of 0.47. The second run produced 1 tree with a L of 713, a CI of 0.41, and a RI of 0.48. The last run produced 1 tree with a L of 905, a CI of 0.30, and a RI of 0.18. From these 3 cladograms, 5 GTs (GT1-5) were generated (Table 1, Fig. 1). In the Nearctic region, GT1 included the Sonora biogeographic province and was supported by the occurrence of the genera Eremarionta and Helmintoglypta (Table 1, Fig. 1A). GT3 was located between the California and Baja California biogeographic provinces and was supported by the occurrence of 5 genera (Table 1, Fig. 1B). GT2 and GT4 were present throughout the Nearctic and Neotropical regions, as well as the Mexican Transition Zone. GT2 crossed 10 biogeographic provinces and was supported by 7 genera (Table 1, Fig. 1A); GT4 was supported by 4 genera and was present in 12 of the 14 biogeographic provinces, absent only in the Balsas Basin and Sierra Madre del Sur provinces (Table 1, Fig. 1B). GT5 was supported by 8 genera, being the only track that is supported by the distribution of freshwater gastropods (Aroapyrgus and Balconorbis) and was located near the Atlantic slope; it included the Tamaulipas, Veracruzan, Sierra Madre Oriental, and Transmexican Volcanic Belt provinces (Table 1, Fig. 1C). A total of 86 individual nodes (IN) were identified: 13 in the Nearctic zone, 36 in the Mexican Transition Zone, and 37 in the Neotropical zone, on the other hand, a single generalized node (GN) was recovered by the intersection of the Generalized Tracks 1 and 4 in northwestern Mexico (Fig. 1D). The analysis for the land snails was performed with 149 genera and the results were similar to that obtained with the complete dataset (data not shown), finally, for the analysis of freshwater taxa a total of 38 genera were used, a single tree was recovered for the PAE-PCE analysis with a length (L) of 137, a consistency index (CI) of 0.44, and a retention index (RI) of 0.47. Two generalized tracks were recovered, the first extends from southern Chiapas towards Tabasco and Veracruz and is supported by the distribution of 3 genera, meanwhile, the second track extends from northern Mexico through the Atlantic coast to the Yucatán Peninsula and is supported by the distribution of 12 genera, only 2 individual nodes were recovered (Supplementary material Table S2, Fig. S1, available at: https://github.com/HOmarMejiaG/diversity-patterns).

Figure 1 Panbiogeographic generalized tracks and nodes recovered in this study for the analysis of Mexican continental gastropods: A) GT1 and GT2, B) GT3 and GT4, C) GT5, and D) the only GN (generalized node) recovered by overlapping of the GT1 and GT4 and the other 86 IN (individual nodes). 

Of the 189 grid-cells, 13 (7%) of them had no species records. Most of the other grid-cells were characterized by low diversity values (<6% of total species). However, 3 grid-cells contained over 10% of continental gastropod species and could be considered hotspots. The first hotspot contained 108 species (10% of the total) and was located in the Veracruz moist forest and the Oaxacan montane forest ecoregions. The second hotspot is located in Veracruz moist forest and Veracruz montane forest ecoregions containing 130 species (12% of the total). The third hotspot was located in the Oaxaca montane forest, Sierra Madre de Oaxaca oak-pine forest and Petén Veracruz moist forest ecoregions and include 229 species (21%), a fifth of Mexico’s continental gastropod diversity (Fig. 2A). For the land snails, 4 grid cells could be considered as hotspots, the number of species ranged from 7.1% to 13.4%, the location of the grid cells match the previously described hotspots, but an additional one is located in Tehuacán Valley matorral and Sierra Madre de Oaxaca pine-oak forest ecoregions (Supplementary material Fig. S2, available at: https://github.com/HOmarMejiaG/ diversity-patterns). In the freshwater snails, 3 hotspots were recognized, the first one contained 18 species (9.7 % of the total) and was located in Pánuco tropical and subtropical coastal rivers freshwater ecoregion, the second hotspot contained 15 species (8.1%) and was located in the Río Conchos xeric freshwater and endorheic basins freshwater ecoregion, finally, the third hotspot comprised 13 species (7.1%) and is located in the Grijalva-Usumacinta tropical and subtropical coastal rivers freshwater ecoregion (Supplementary material Fig. S3, available at: https://github.com/HOmarMejiaG/diversity-patterns)

The Shannon diversity values tended to be more geographically homogeneous, but this analysis indicates grid-cells with high diversity values (3.6-4.4) in Hidalgo, Puebla, Queretaro, Nuevo León, Coahuila, Oaxaca, and Veracruz for the complete dataset and the land snails. On the other hand, for the freshwater snails the highest diversity cells (2.1-2.6) are located in Coahuila, San Luis Potosí and Veracruz (Supplementary material Fig. S4, available at: https://github.com/HOmarMejiaG/diversity-patterns)

Of the 1,075 continental gastropod species studied, 814 (75.7%) are endemic to Mexico, but only 15 species (1.4%) (all terrestrial) had widespread distributions, occupying 10 or more grid-cells: Drymaeus multilineatus (10 grid-cells), Praticolella berlandieriana (10), Bulimulus unicolor, Allopeas beckianum (11), Choanopoma largillierti, Rabdotus schiedeanus (11), Gastrocopta pellucida (11), Drymaeus emeus (12), Helicina chrysocheila (15), Neocyclotus dysoni (15), Orthalicus princeps (17), Thysanophora horni (17), Succinea luteola (18), Hawaiia minuscula (19) and Rabdotus alternatus (21). Meanwhile, 433 (39.8%) of the species occupied less than 10 grid-cells, and the majority of species (586, or 58.7%) are restricted to only 1 grid-cell (Fig. 2B; Supplementary material Table S1, available at: https://github.com/HOmarMejiaG/diversity-patterns)

Figure 2 Spatial diversity patterns of Mexican continental gastropods: A) proportion of the 1,075 gastropods species present in each of the 189 grid-cells that comprise Mexico, the numbers on cells represent the 3 hotspots recovered; B) average proportional range size of Mexican gastropods for each of the 189 grid-cells. 

The Nearctic region comprises 5 biogeographic provinces with an estimated area of 959,205 km² and includes 280 species of continental gastropods of which 204 are endemics; the Neotropical region comprises 4 biogeographic provinces with an estimated area of 433,558 km² and include 452 species of continental gastropods of which 326 are endemics, finally, the Mexican Transition Zone comprises 5 biogeographic provinces with an estimated area of 591,901 km² and include 343 species of continental gastropods of which 284 are endemics. For the land snails, the Nearctic region included 250 species (156 endemic), the Neotropical region 460 species (281 endemic) and the MTZ 185 species (175 endemic); meanwhile, for the freshwater snails, the Nearctic region included 50 species (45 endemic), the Neotropical region 90 (44 endemic) and the MTZ 40 species (13 endemic).

Discussion

There are several critiques for PAE as a panbiogeographic method for recovering historical biogeographic patterns, such as the method does not take into account the phylogenetic relationships (^{Humphries, 1989}), the use of a posteriori optimality criterion (^{Szumik, 2002}) and the use of artificial delimited areas (^{Peterson, 2008}), among others. Nevertheless, most of the critiques arise from the use of PAE as a method for recovery of historic relationships among areas and not as a method for the recognition of areas of endemism (^{Nihei, 2006}). Areas of endemism are strong evidence of the historical structure of biotic components and represent a basic step in evolutionary biogeography, where the PAE PCE method appears as the best approach for the identification of generalized tracks (^{Morrone, 2014}).

The biogeographic patterns recovered in this study were similar to those previously established in other biological groups with the same methodological approach of PAE PCE (^{Contreras-Medina et al., 2007}). Beyond agreement with previous research, there are also other pieces of evidence that support the patterns recovered in this study.

First, the GTs of northwestern Mexico (GT1 and GT3; Fig. 1, Table 1) were supported by 7 genera of Nearctic affinity restricted to the arid ecosystems of the Sonora, California, and Baja California provinces (^{Naranjo-García & Fahy, 2010}) and are similar to GTs recovered for other taxa (^{Escalante et al., 2018}; ^{Rojas-Soto et al., 2003}). The heterogeneous landscape, as well as geological changes that occurred in the Baja California Peninsula 6-15 MYA can explain these patterns (^{Escalante et al., 2018}; ^{Naranjo-García, 1988}; ^{Riddle et al., 2000}).

Second, the geographic distribution of GT2 and GT4 (Fig. 1) are congruent with the typical Neotropical pattern proposed by ^{Halffter (1987)}, the lowland Mesoamerican track proposed by ^{Morrone and Márquez (2001)}, and the patterns recovered in other groups (^{Escalante et al., 2018}; ^{Gonzalez-Ávila et al., 2017}; ^{Luna-Vega et al., 2001}; ^{Márquez & Morrone, 2003}; ^{Rosas-Valdez & Pérez- Ponce de León, 2008}). The Neotropical affinity groups may extend beyond the coastal regions and low altitudes, as evidenced by the presence of Neotropical affinity species above 1,400 m such as the genera Allopeas and Drymaeus in GT2 that could reach up to 2,059 and 3,698 m, respectively, and the genera Echinix and Polygyra in GT4 that occur up to 1,485 and 2643 m, respectively. This distribution pattern in highland regions has also been found in other groups, such as volant mammals, and is characteristic of a transitional fauna that may have arisen either by the dominance of endemic or widely distributed species or by a balance between the number of Neartic and Neotropical species in grid cells (^{Ortega & Arita, 1998}).

Furthermore, the support of GT2 by the genera Drymaeus, Euglandina, and Orthalicus demonstrate the role of Sierra Madre del Sur, the Sierra Madre Occidental, and the Transmexican Volcanic Belt as important biogeographic barriers to dispersal for several terrestrial gastropod families (^{Mejía & Zúñiga, 2007}; ^{Naranjo-García & Fahy, 2010}).

The panbiogeographic tracks recovered for the freshwater snails were similar to the tracks recovered for the land snails; although an incomplete species sampling over the whole country could lead to misleading historical patterns (^{Brooks & Van Veller, 2003}; ^{Nihei, 2006}), our results support that a panbiogeographic analysis is useful as a preliminary approach to state primary biogeographic homology hypothesis in groups where phylogenetic analysis are scarce.

The species richness in tropical regions is explained by both ecological and historical factors. The number of species recovered in this paper in the Neotropical region (5 times greater than the Nearctic region and 2.5 times higher than the MTZ) were also recovered in other biological groups such as different phyla of freshwater organisms (^{Quiroz-Martínez et al., 2014}) and mammals (^{Gómez-Ortiz et al., 2017}), supporting that the Neotropical region comprises the higher diversity in Mexican fauna. Additionally, the results of this study highlight the complexity of the MTZ (^{Corona et al., 2009}; ^{Escalante et al., 2018}; ^{Halffter & Morrone, 2017}) that has been considered either as a high richness species area (^{Aguilar-Méndez et al., 2021}) or as an area of high phylogenetic diversity, (^{Gómez-Ortiz et al., 2017}) that are reflected in the high number of total (343) and endemic (284) species of continental gastropods. Thus, the Mexican Transvolcanic Belt provide a set of factors that are related with high speciation rates in land snails such as a topographic mosaic, rain shadow effects, heat budget differences, vegetation cover and effective refugia from floods that serve as shelter areas (^{Solem, 1984}). The diversity hotspots of freshwater snails recovered in this study were similar to the previously reported by ^{Czaja et al. (2020)}, for Cuatro Cienégas-Sabinas system and the Río Verde-Guayalejo, with the exception of the Río Conchos- Guzmán system hotspot that has not been recovered in this work, although we were able to recovered a grid in the Grijalva that has been postulated as a potential hotspot (^{Czaja et al., 2020}). The main difference between both approaches is that ^{Czaja et al. (2020)} used the complete freshwater ecoregions meanwhile in this work we used grid cells that reduce the size area effect allowing a more precise comparison.

In particular, the nodes gathered in this study were geographically similar to the ones proposed by ^{Morrone and Márquez (2008)}, who defined 15 nodes in the boundaries of some Mexican biogeographic provinces; they are also similar to 3 of the 5 freshwater snails endemicity hotspots suggested by ^{Czaja et al. (2020}) and to the highest species diversity areas for terrestrial gastropods recovered in previous regional studies (^{Correa-Sandoval et al., 2007}; ^{Naranjo-García & Fahy, 2010}).

In this study, a single GN was recovered in Northwestern Mexico and only 2 INs were recovered for the Mexican Pacific slope. In contrast, 2 regions along the Gulf of Mexico basin had very high node concentrations. The intersection of the Tamaulipas, Sierra Madre Oriental, Veracruzan, and Transmexican Volcanic Belt biogeographic provinces had up to 60 INs; this coincides with the proposal by ^{Correa-Sandoval (2003)}, who suggested the Sierra Madre Oriental malacological province for its distinctive fauna. The nodes identified in the central area of the Veracruzan province are geographically similar to the endemicity hotspots for freshwater gastropods in the Río Verde, Río Guayalejo, and Río Papaloapan regions, as previously defined by ^{Czaja et al. (2020)}.

The results of spatial richness patterns obtained in this study are similar to those recovered in other taxa (^{Aguilar-Aguilar et al., 2008}; ^{Koleff et al., 2008}; ^{Ochoa-Ochoa et al., 2012}). However, only 3 grid-cells contained more than 10% of the species diversity, which is in contrast to previous research on vertebrates, where there were a large number of grid-cells with at least these diversity values (^{Koleff et al., 2008}). This could be the product of environmental factors, which have a considerable effect on species richness patterns for taxa with low dispersal rates because they have geographically restricted distributions and high levels of endemism (^{Clements et al., 2008}). Finally, the INs and highest species diversity values detected in this study coincide with conservation priority areas previously defined either by species richness or by the occurrence of biotic elements of different origins (^{Corona et al., 2009}).