Subterranean Waters of Yucat á n Peninsula, Mexico Reveal Epigean Species Dominance and Intraspeciﬁc Variability in Freshwater Ostracodes (Crustacea: Ostracoda)

: The Yucat á n Peninsula is a karstic region, rich in subterranean environments with a diverse crustacean stygobiont fauna. In order to gain insights into the biological evolution of the subterranean environments of this region, we evaluated the ostracode species composition of caves and cenotes in ﬁve independent sampling campaigns (2008, 2013, 2017–2019). Using morphometric analyses, we evaluated inter-population morphological variability; using molecular analysis based on mitochondrial COI and nuclear 18S rDNA, we evaluated genetic differentiation in selected species. The observed fauna is composed of 20 (epigean) species, presenting a lack of strict stygobionts. Morphometric analyses discriminated up to three morphotypes in each of the three most abundant species: Cytheridella ilosvayi , Alicenula sp. and Cypridopsis vidua . High intraspeciﬁc morphological variability was found either in shape or size. Phylogenetic analysis based on COI demonstrated the existence of three lineages on C. ilosvayi , with high support (>0.9). The 18S rDNA sequences were identical among individuals of different populations. A lack of congruence between the genetic markers precluded us from postulating speciation in subterranean environments. It is likely that Late Pleistocene—Early Holocene climate variation related to sea level and precipitation was forcing agent for epigean ostracode dominance in subterranean environments of the Peninsula.


Introduction
Subterranean waters are critical environments for biodiversity. The combination of unique abiotic features such as stable and relatively low temperatures, geologically driven water chemistry, space restriction and low energy and food availability, challenge the success of most taxa [1][2][3][4]. Aquatic species have been able to colonize subterranean environments through morphological, behavioral and physiological adaptations [4][5][6], but given the highly selective environment, the number of species per site (i.e. α diversity) and their ecological interactions are usually limited [1,[7][8][9]. Species inhabiting subterranean waters are usually isolated, with a limited capacity of dispersal. These conditions may drive large differences in species composition between regional sites (high β diversity) [10,11]. Speciation and endemicity are therefore common in these environments and are promoted by events such as vicariance and spatial ecological partition of the environment [11][12][13].
Global patterns of biodiversity and endemicity in subterranean waters are largely varied among regions and are influenced by factors such as subterranean water availability, the geological and climatic history of the region and the capacity of the organisms to succeed in these environments [2,[13][14][15]. In tropical regions for example, the biodiversity of obligate aquatic cave dwellers is overall lower than that in temperate regions [16,17], but large tropical regions remain unexplored.
The Yucatán Peninsula is a calcareous platform located in southern Mexico; it is partially surrounded by the Gulf of Mexico and the Caribbean Sea. The Peninsula has one of the largest underwater cave systems in the world. This system is characterized by extensive interconnected caves throughout the Peninsula [18]. Physical and chemical properties of the subterranean system, such as water level and ionic water composition are strongly influenced by the proximity and connectivity with marine environments, particularly the Caribbean Sea [19][20][21].
In the Yucatán Peninsula groundwaters become frequently exposed by the continuous dissolution of limestone [22][23][24] and underground waterways become exposed. The resulting karstic phenomenons are locally named "cenotes" ("dzonot" in mayan language, meaning "hollow with water"). Currently, more than 8000 cenotes are estimated to be sparsely distributed in this region (Secretaría de Desarrollo Sustentable, Yucatán census). In the northern portion of the Peninsula, cenotes are almost the unique source of freshwater and represent an important interaction between groundwater (hypogean) and surface environments (epigean) [25,26]. Given their geomorphological features, subterranean habitats of the region can be caves (partial or completely flooded) or cenotes (deep wells), sometimes with extensive networks of cave passages [27]; in some of them, marine waters can underlie the freshwater table [28].
Biotic diversity in subterranean waters of the Peninsula is dominated by crustaceans and the region is recognized as one of the foremost sites for this group with more than 40 obligate subterranean species from fresh and anchialine waters, also known as "stygobiont" species [27].
Freshwater ostracodes are micro crustaceans with a maximum length of 7.5 mm [29] (but mostly <5 mm), characterized by having their body enclosed between two calcareous valves [30]. Ostracodes are typical of almost all aquatic systems, including ephemeral water deposits, subterranean and interstitial waters, semiterrestrial environments and a wide variety of epigean water systems [31]. In subterranean environments, ostracodes have been documented worldwide [32][33][34][35][36][37] and are characterized by high levels of endemism and a number of morphological adaptations including reduction of size, pigmentation loss, elongation of sensorial structures and reduction of setae [33]. Globally, the study of subterranean ostracodes is still under development, but evidence from karstic systems in Europe [32,33,38] and the arid regions of Australia [39] have demonstrated a high level of diversity. In Western Australia (Pilbara region), subterranean ostracode fauna can be more diverse than that in surficial environments [39][40][41].
Studies of ostracode and crustacean species from cenotes of the Yucatán Peninsula have revealed high species diversity, with a variety of ecological interactions between epigean, hypogean, and anchialine fauna [26,27,[42][43][44]. For ostracodes, a sampling campaign undertaken in cenotes in 1936, revealed the presence of 23 freshwater species in Yucatán, out of which 13 of them were described as new to science [45]. More recently, sampling campaigns undertaken in 1989 [46] and 1998 [47] in cenotes located along the Caribbean Sea coast (<5 km separated from the coast) revealed the presence of two anchialine species in the water column, below the halocline. On the Yucatán Peninsula, therefore, integrated studies aimed to understand species compositions are required to further test ecological and evolutionary hypotheses, such as epigean and hypogean species interactions, colonization, speciation events and endemicity. In this study, (i) we evaluate ostracode faunal composition from subterranean habitats (cavern type and cenotes) of the Yucatán Peninsula; (ii) we evaluate morphological variability in valves and using molecular analyses, we assess the genetic differentiation occurring among subterranean ostracodes, and (iii) we discuss the role of past climate variations (since the Late Pleistocene) in structuring the composition of subterranean ostracode fauna of the Peninsula.

Materials and Methods
We sampled a total of 32 cenotes in 25 settlements of the Yucatán Peninsula (Yucatán, Campeche and Quintana Roo federal states), Mexico, during independent sampling campaigns performed in 2008 (eight sites) and 2013 (ten sites) and between 2017 and 2019 (five sites in 2017, three sites in 2018 and six sites in 2019) ( Figure 1, Table S1). Each site was visited only once, as each campaign took place in different regions of the Peninsula.
3 composition from subterranean habitats (cavern type and cenotes) of the Yucatán Peninsula; (ii) we evaluate morphological variability in valves and using molecular analyses, we assess the genetic differentiation occurring among subterranean ostracodes, and (iii) we discuss the role of past climate variations (since the Late Pleistocene) in structuring the composition of subterranean ostracode fauna of the Peninsula.

Materials and Methods
We sampled a total of 32 cenotes in 25 settlements of the Yucatán Peninsula (Yucatán, Campeche and Quintana Roo federal states), Mexico, during independent sampling campaigns performed in 2008 (eight sites) and 2013 (ten sites) and between 2017 and 2019 (five sites in 2017, three sites in 2018 and six sites in 2019) (Figure 1, Supplementary Table S1). Each site was visited only once, as each campaign took place in different regions of the Peninsula.  Table S1. The light area within the Yucatán state represents the Geohydrological Reserve in Yucatán, while the dark area depicts the urban extension of the city of Merida.

Sampling, Valve Extraction, Counting and Identification
Biological samplings were carried out using two different methods. During the years 2008 and 2013, sampling was undertaken by the Instituto Tecnológico de Chetumal (ITCH), Chetumal, Quintana Roo, Mexico and El Colegio de la Frontera Sur (ECOSUR), Unidad Chetumal, Chetumal, Quintana Roo, Mexico. During these campaigns, biological samples were collected from the littoral zone, water column and deepest bottom of cenotes. At littoral areas, samples were taken from submerged vegetation using a hand net of 250 µm mesh size. In the water column, samples were collected with vertical tows and horizontal trawls with a 150 µm mesh and a 20 cm wide net. Sediment samples were collected with an Ekman grabfrom the deepest bottom; only the uppermost 2 cm of sediment of each grab were collected to avoid fossil or subfossil shells. Samples fixation was done  Table S1. The light area within the Yucatán state represents the Geohydrological Reserve in Yucatán, while the dark area depicts the urban extension of the city of Merida.

Sampling, Valve Extraction, Counting and Identification
Biological samplings were carried out using two different methods. During the years 2008 and 2013, sampling was undertaken by the Instituto Tecnológico de Chetumal (ITCH), Chetumal, Quintana Roo, Mexico and El Colegio de la Frontera Sur (ECOSUR), Unidad Chetumal, Chetumal, Quintana Roo, Mexico. During these campaigns, biological samples were collected from the littoral zone, water column and deepest bottom of cenotes. At littoral areas, samples were taken from submerged vegetation using a hand net of 250 µm mesh size. In the water column, samples were collected with vertical tows and horizontal trawls with a 150 µm mesh and a 20 cm wide net. Sediment samples were collected with an Ekman grabfrom the deepest bottom; only the uppermost 2 cm of sediment of each grab were collected to avoid fossil or subfossil shells. Samples fixation was done in situ with 96% ethanol and then preserved on ice. Thereafter, samples were kept in refrigeration at −4 • C for seven days to preserve DNA from degradation. Samples were deposited Ostracode specimens sorting and counting were carried out in 50 cm 3 of wet sediment, using a stereomicroscope. Only organisms with complete soft parts or well-preserved valves were extracted. For species identification, organisms with complete soft parts were dissected and mounted in a mixture of glycerol and formaldehyde (1:1). Preparations were sealed with Entellan mounting media. Shells were stored in micropaleontological slides. Valves and appendages were measured with a stage micrometer adapted to an Olympus light microscope. Non-dissected material was preserved in plastic tubes in 70% ethanol. All studied material was deposited at the zooplankton reference collection from ECOSUR, Unidad Chetumal.
Samplings between 2017 and 2019 were conducted by Cenoteando, a research group team from the Unidad Multidisciplinaria de Docencia e Investigación Sisal, Mérida, Yucatán, México (UMDI-Sisal), Facultad de Ciencias, Universidad Nacional Autónoma de México, Mérida, Yucatán, México (UNAM). Sediment samples from different depths of the cenotes and associated submerged cave passages were collected during scientific SCUBA cave dives in 50 mL sample tubes and in 1000 mL plastic bags. The sediments collected corresponded to the uppermost layers of the cave floors. Sample fixation was done in situ, by replacing all water with 96% ethanol; then, the sediment samples were preserved in ice. Large samples were filtered in situ using a 150 µm mesh size net. Thereafter, samples were kept in refrigeration at −4 • C. Species extraction was carried out at the UNAM UMDI-Sisal facilities. Small portions of sediment, diluted with ethanol, was examined using a Nikon SM Z800 stereomicroscope.
Ostracodes and other minute invertebrates such as thermosbaenaceans, stygiomysids, mysids, and amphipods were selected. Only ostracodes with well-preserved valves were extracted. The individuals selected were preserved in labelled sample tubes in 96% ethanol and were sent for further analyses to the Instituto de Geología of the UNAM in Mexico City and to the ITCH. The remaining sediment samples were deposited at the Laboratorio de Ecología y Manejo de Costas y Mares of UMDI-Sisal, UNAM.
In the case of the non-submerged caves X'tacumbilxuna'an, Yumku, and Peba, sediment samples were collected in dry cave passages from freshwater ponds. Samples were processed and deposited, as described previously.
Scanning electron microscope (SEM) analyses were conducted for selected and representative species with a Jeol Jsm-6010 plus/LA (Tokyo, Japan) scanning electron microscope from ECOSUR, Unidad Chetumal.
An Olympus BX51 microscope (Tokyo, Japan) was used for light photographs of specimens using the stacking method.
Species identification was conducted with species keys provided by Karanovic [31] and original species descriptions. Taxonomic classification follows Cohuo et al. [48].

Morphometric Analysis of Ostracode Valves
Morphological variability between populations and subterranean habitats was evaluated with geometric morphometrics. Valve shape analysis was conducted in the three species most frequently found in our sampling sites: Cytheridella ilosvayi Daday, 1905, Alicenula sp., and Cypridopsis vidua (O.F. Müller, 1776). For C. ilosvayi, we evaluated a total of 13 individuals from seven populations (caves Peba and Yumku and cenotes Huul K'in, Yokdzonot, Azul (in Campeche), Azul (in Quintana Roo) and Xoch). We decided to use female right valves for this analysis as they were the most abundant and in better preservation state. Although the left valve is more stable in ostracodes, for C. ilosvayi, right valve has demonstrated morphological stability and accurate results on inter-and intra population comparisons [49,50]. For Alicenula sp., specimens were mostly identified , and X'batun (San Antonio Mulix) were used. Only adult well-preserved valves were used in the analysis. Given that Alicenula sp. and C. vidua lack enough landmarks (type I and II) in their valve surface, we used for all three species an outline approach based on the elliptic Fourier method (EFA). The outlines of the specimens were created as two curves, starting at the maximum anterior and posterior curvature of the valves. Considering the differences in the valve shapes of the target species, 104 semilandmarks were drawn, equidistantly spaced for C. ilosvayi, 106 for Alicenula sp., and 98 for C. vidua. Semilandmarks were created manually in the tpsutil32 and tpsDig2 software packages ( [51]; https://life.bio.sunysb.edu/morph/).
A set of 20 harmonics that best describe the shape of the valves were calculated with the EFA3D software (3D elliptic Fourier; https://life.bio.sunysb.edu/morph/softoutlines.html). The harmonics were invariant to the conditions of the valves, such as size, rotation, and position. A principal component analysis (PCA) with the harmonic scores was conducted to visualize the patterns of shape variations, and EFA coefficients were used for valve shape reconstruction. Significant differences in shape morphologies between morphotypes were assessed with MANOVA and canonical variates analysis (CVA). For MANOVA, morphotypes represented by two or more specimens were included in the analysis, and Bonferroni correction was applied to the resulting p-values. CVA was performed using the PCA scores, and, as groups of interest, we used the associations resulting from the harmonics scores in the PCA morphospace. As CVA interpretations can be strongly affected by the number of variables (when variables exceed the number of specimens per group) [52], we used only those components that describe 90% of data variability. PCA, MANOVA and CVA were performed with the PAST version 3 software (Oslo, Norway) [53].
The resulting geometric morphometrics analysis did not uncover differences in individual size. A length-height analysis was then performed to provide a graphical representation of size differences among populations. Valve size analysis was conducted using 17 specimens of C. ilosvayi from eight populations, 19 specimens of Alicenula sp. from eight populations, and 30 specimens of C. vidua from seven populations. For this analysis, we used adult well-preserved and partially broken valves, as far as the length and height were clearly observable. We measured the total length (L) and height at the middle of the valve (H) of each specimen in each of the three groups. The L/H ratio was then calculated to identify changes in the lateral shape of the valves. The normality of the measurements was assessed by Shapiro-Wilk tests. As the data sets did not fit a normal distribution, differences between populations in length, height, and L/H ratio were assessed with non-parametric Kruskal-Wallis tests. In order to provide a graphic representation of the relationship between length and height for all species, scatter plots were generated. Descriptive statistics such as mean, median, standard deviation, and maximum and minimum measures of the length and height of each group were calculated.

Molecular Analysis: DNA Extraction, Amplification and Sequencing
To test genetic differentiation in subterranean waters (cenotes and caves), we performed phylogenetic analyses, using two independent genetic markers-the mitochondrial COI and nuclear 18S rDNA gene. Molecular analysis was conducted exclusively on C. ilosvayi, as it is the unique species that we were able to collect complete specimens (valves + soft parts) of. Given the practical absence of genetic sequences for this species, and to provide a reference framework for species delimitation, we included individuals of six epigean populations. Subterranean habitats corresponded to Cave Yumku, Cave Peba, Cenote Azul (in Campeche), Cenote Azul (in Quintana Roo) and Cenote El Padre. Epigean environments are represented by lakes characterized by different trophic conditions. The following lakes are oligotrophic: Bacalar and Señor, while the remaining ones-Tres Garantías, Campeche, Caobas and Perdida-are eutrophic lakes.
From each of these systems, three individuals were used for molecular analysis. The valves of each organism were removed and stored in micropaleontological slides, and the soft parts were further used for DNA extraction.
A volume of 4 µL of DNA was used as a template for amplifications. PCR products were visualized in 1.5% agarose gel, and sequencing was done with a BigDye ® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Foster City, CA, USA) on an ABI 3130 automated DNA sequencer. Chromatographs were manually checked, and sequences were edited using CodonCode Aligner (v. 9.01, Codon Code Corporation, Centerville, MA, USA) where necessary. All new sequences were deposited in GenBank (https://www.ncbi.nlm. nih.gov/genbank) with the following accession numbers: MW018845-MW018862 for COI and MW018806-MW018816 for 18S rDNA (Table S2).

Molecular Analysis: Sequence Alignment and Phylogenetic Analysis
COI and 18S sequences were aligned in MEGA software version 10.1.8 [58], the Clustal W algorithm was used for COI [59], while non-coding sequences from 18S were aligned with the MUSCLE algorithm. Sequences were trimmed manually when necessary. To avoid stop codons, COI sequences were translated into amino acids using the ORF finder from the NCBI page. All COI and 18S sequences are newly obtained in this study. With the purpose of comparison and in order to provide robustness to the analysis, we also used five sequences from closely related species of C. ilosvayi: two COI sequences are from Paracythereis opesta (Brehm, 1939), which are newly obtained in this study, two sequences are from a putative Cytheridella species from GenBank (accession numbers: Podocopida sp.1 and 2-MG449867 and MG449978) and an 18S sequence from Metacypris digitiformis Smith amd Hiruta, 2004 (accession numbers: AB674964). Two ostracode species where used as the outgroups for the COI and 18S data sets, respectively (accession numbers: Limnocythere inopinata (Baird, 1843)-AJ534412 and Limnocythere sp.-AB076635).
For phylogenetic analysis, the best-fitting model of molecular evolution and the partition scheme, under the corrected Akaike information criterion (AICc) [60] was estimated for the COI and COI-18S concatenated database using Partition Finder 2 on XSEDE [61].
A Bayesian approach was implemented using MrBayes on XSEDE [62]. According to the best-fitting model, the following partitions were included: SYM + I for the first COI codon position, GTR for the second COI position, K80+G for the third COI position, and GTR for the 18S gene segment. Two runs of one million generations and four chains sampled for every 1000 generations were used. The first 25% of the resulting trees corresponding to the initial phase of the Markov chain were discarded as burn-in. All analyses were implemented in in CIPRES Science Gateway [63].

Ostracode Species Composition in Subterranean Waters of the Yucatán Peninsula
A total of 20 ostracode species, belonging to four families, were recovered from subterranean habitats (caves and cenotes) of the Yucatán Peninsula ( Figure 2, Table 1). We found  Table 1 summarizes the taxonomic classification and distribution of the species in the study areas. The species C. vidua and C. ilosvayi were the most widely distributed with 15 and 13 records, respectively. For all other species we collected mostly empty shells, mainly from bottom sediments. As a result, we were unable to confirm the taxonomic classification for at least eight species, which remain in open nomenclature.
we collected mostly empty shells, mainly from bottom sediments. As a result, we were unable to confirm the taxonomic classification for at least eight species, which remain in open nomenclature.   the influence of the marine environment on species composition. Species richness in other cenotes varied from one to three species. In cave environments such as Yumku, Peba and, X'tacumbilxuna'an, between one and two species were observed.
All taxa recorded in the sampling sites are epigean fauna. The species composition is integrated with species widely distributed in tropical America or cosmopolites such as Diaphanocypris meridana (Furtos, 1936) (Furtos, 1936). The latter species (K. xanabanica) is the only species currently known to be distributed exclusively in cenotes of the northern Yucatán Peninsula.

Shape and Size Morphometric Analysis of Selected Ostracode Species from Caves and Cenotes
The principal component analysis based on the first and second components of elliptic Fourier coefficients, explain more than 59% of total data variation for C. ilosvayi ( Figure 3A), Alicenula sp. (Figure 3B), and C. vidua ( Figure 3C). For C. ilosvayi, the first and second principal components of the PCA explained 85% of total data variation. Patterns in morphospace suggested the presence of three morphological forms that seem not to be randomly distributed but, rather, associated with habitat types ( Figure 3A). Cave systems Yumku and Peba were associated with the first morphotype (type CIL I), in which the outline describes elongated organisms with a well-developed brood pouch, exceeding the dorsal margin, and with the anterior margin slightly pointed. The second morphotype (type CIL II; Figure 3A) is related to cenotes, with the water surface directly in contact with the exterior. These organisms are narrow, elongated, and square-shaped, with an almost straight dorsal margin (without overlapping of the brood pouch) and a broadly rounded anterior margin. The third morphotype (type CIL III; Figure 3A) is considered a valve deformation, as it was observed in a single specimen. It was characterized by a narrowly developed and totally compressed brood pouch area. This specimen belongs to cenote Yokdzonot and contrasted markedly with the other forms. Due to the low number of specimens recovered, we cannot provide a reliable conclusion on the nature of this form. A MANOVA test showed significant differences (Wilks' lambda test: Λ Wilks = 0.0014, F = 151.7, p < 0.006) between CIL I and CIL II morphotypes. Jackknifed rate classification was 83% correctly assigned for this species. Positions of CIL I and CIL II morphotypes on the CVA plot were clearly differentiated and no overlap was observed ( Figure 4A).
The first and second components of the PCA explained 94% of the total shape variation of the ostracode Alicenula sp. Two forms were recovered from the morphospace, mainly associated with PC1 ( Figure 3B). Similarly, to C. ilosvayi, the forms were not randomly associated but clearly corresponded with well-defined water systems. The cenotes Dzonbacal and Dzonotila were related to the morphotype AI, whereas all other cenotes were associated with the second morphotype AII ( Figure 3B). An observed difference in the forms was the anterior margin, which is rounded in the AI form and slightly pointed in the AII form. The MANOVA test revealed significant differences (Wilks' lambda test: Λ Wilks = 0.0001, F = 2164, p < 0.0004) between the AI and AII morphotypes. Jackknifed classification rate was 100% correctly assigned. On the CVA plot, morphotypes did not overlap ( Figure 4B).
In C. vidua, the first and second components explained only 59% of shape data variability, but a notable trend was depicted in the PCA morphospace, with three morphotypes ( Figure 3C). Two of them (type CIV I and CIV II) were distributed in cenotes with waters directly in contact with the exterior, whereas the third morphotype CIV III, inhabits in the cave environment of X'tacumbilxuna'an. A MANOVA test showed significant differences between C. vidua morphotypes (Wilks' lambda test: ΛWilks = 0.0003, F = 33.92, p < 0.003). Classification rates were 100% correctly assigned in the confusion matrix. The CVA plot clearly discriminates three morphotypes, as no overlaps were observed ( Figure 4C).  Table S1. Abbreviation BP = brood pouch. Black arrows represent anterior the position of valves.   Table S1.
As elliptic Fourier analysis did not recover variability in valve size, a morphome analysis based on length and height was conducted. Data dispersion in the scatterplo C. ilosvayi showed overlapping positions between the populations of different hab types (cenotes and caves; Figure 5A). Comparisons of length, height and L/H r showed no significant differences between the size of the eight populations evalua   Table S1.
As elliptic Fourier analysis did not recover variability in valve size, a morphometric analysis based on length and height was conducted. Data dispersion in the scatterplot for C. ilosvayi showed overlapping positions between the populations of different habitat types (cenotes and caves; Figure 5A Table S1). Dashed circles indicate cave populations. .

Molecular Analysis of Cytheridella ilosvayi Epigean and Hypogean Populations
For COI alignment after excluding the primer region and editing the sequence, the consensus length of sequences was 442 bp (base pairs) and no insertions, deletions, or stop codons were observed. For C. ilosvayi, the final COI data set was integrated by 16 new  Table S1). Dashed circles indicate cave populations.

Molecular Analysis of Cytheridella ilosvayi Epigean and Hypogean Populations
For COI alignment after excluding the primer region and editing the sequence, the consensus length of sequences was 442 bp (base pairs) and no insertions, deletions, or stop codons were observed. For C. ilosvayi, the final COI data set was integrated by 16 new sequences of three populations from subterranean habitats and six populations from epigean environments and also 2 sequences from GenBank (Podocopida).
Phylogenetic analysis based on COI and Bayesian Inference (BI) resulted in a tree topology that discriminates three lineages with high support (>0.9; Figure 6). Populations of freshwater and oligotrophic lakes and cenotes located relatively close to the Caribbean Sea coast, constituted the first lineage. Populations of lakes from southern Yucatán Peninsula, characterized by extensive littoral areas and eutrophic conditions, are the second lineage. The third lineage is composed of populations from cave environments Yumku and Peba, both in the northern part of the Peninsula.
sequences of three populations from subterranean habitats and six populations from epigean environments and also 2 sequences from GenBank (Podocopida).
Phylogenetic analysis based on COI and Bayesian Inference (BI) resulted in a tree topology that discriminates three lineages with high support (> 0.9; Figure 6). Populations of freshwater and oligotrophic lakes and cenotes located relatively close to the Caribbean Sea coast, constituted the first lineage. Populations of lakes from southern Yucatán Peninsula, characterized by extensive littoral areas and eutrophic conditions, are the second lineage. The third lineage is composed of populations from cave environments Yumku and Peba, both in the northern part of the Peninsula.  Table S2).
The 18S data set comprised 11 sequences of 1078 to 1170 pb length, and they show extremely low variation between populations; most sequences were identical among species. This made it impossible to differentiate intraspecific groups, based on this gene. The concatenated COI-18S data set consisted of 16 sequences of 2238 pb length, and the resulting tree topology was highly concordant with COI tree topology. The lineages discriminated by independent COI analysis were also uncovered by the COI-18S tree ( Figure S1). This is expected as the COI sequences were responsible for the variability of COI-18S data set.

Low Diversity and Epigean Species Composition of Ostracode Assemblages in Subterranean Environments of the Yucatán Peninsula
The Yucatán Peninsula is a region where the biological diversity of freshwater ostracodes in epigean environments is high (>30 spp.), compared to other tropical regions of the American continent such as Central America (⁓25 spp.), the Caribbean Islands (⁓40 spp.) and the transition zone between the Neotropics and the Nearctic (Central Mexico; ⁓20 spp.) [45,48,54,64,65]. This region is indeed considered a site of (micro) endemism, Figure 6. Phylogenetic relationship in Cytheridella ilosvayi populations of oligotrophic and eutrophic lakes, cenotes and caves in the Yucatán Peninsula. COI tree topology is based on Bayesian Inference (BI). Numbers above branches represent Bayesian posterior probabilities. Abbreviations on tree topology correspond to GenBank species identification code (see Table S2).
The 18S data set comprised 11 sequences of 1078 to 1170 pb length, and they show extremely low variation between populations; most sequences were identical among species. This made it impossible to differentiate intraspecific groups, based on this gene. The concatenated COI-18S data set consisted of 16 sequences of 2238 pb length, and the resulting tree topology was highly concordant with COI tree topology. The lineages discriminated by independent COI analysis were also uncovered by the COI-18S tree ( Figure S1). This is expected as the COI sequences were responsible for the variability of COI-18S data set.

Low Diversity and Epigean Species Dominance on Ostracode Assemblages of Subterranean Environments of the Yucatán Peninsula
The Yucatán Peninsula is a region where the biological diversity of freshwater ostracodes in epigean environments is high (>30 spp.), compared to other tropical regions of the American continent such as Central America (~25 spp.), the Caribbean Islands (~40 spp.) and the transition zone between the Neotropics and the Nearctic (Central Mexico;~20 spp.) [45,48,54,64,65]. This region is indeed considered a site of (micro) endemism, with more than 70% of the species restricted to single lakes (e.g. Cypretta spinosa Cohuo et al. 2013) [61] or limnological regions (e.g. S. intrepida, C. petenensis) [66]. In subterranean waters, however, the observed αdiversity is relatively low, with an average of two ostracode species per site and with a total of 13 genera and 20 species in the 32 systems evaluated. The low biodiversity in ostracodes of subterranean waters of the Yucatán Peninsula, at this stage, is coincident with the global trend observed for subterranean assemblages, in which species richness is higher in temperate regions of higher latitude, compared to the tropics [16,67]. In groundwaters of the karstic regions of Europe such as Romania, a total of 19 genera and 42 ostracode species have been recorded [68] in an area of a similar extent to Yucatán. The number of ostracode species found in our study, may, however, represent a gross underestimation of the biological diversity of subterranean waters in the Yucatán Peninsula, as our sampling sites represent less than 1% of the total subterranean systems available in the region. Particularly, the presence of ostracodes in all systems evaluated (caves and cenotes and its subterranean passages) may be an indicative of a higher diversity in the region, as far as more systems are evaluated.
Surprisingly, the ostracode assemblage in the cenotes and caves of the Peninsula was constituted exclusively by epigean species, while stygobionts are lacking. All of the fauna collected has also been observed in lake environments and about 40% of these species display wide distributions in the tropical America or are cosmopolitan (generalists) [48,69]. At the regional scale, in the northern Neotropics (southern Mexico, Central America, and the Caribbean), subterranean ostracodes have been studied mainly in the Yucatán Peninsula [45][46][47] and the Antilles (including the Caribbean coast of Venezuela) [70][71][72][73]. These studies have demonstrated that the region harbors stygobiont fauna in subterranean waters in both fresh and anchialine environments. In freshwaters, three endemic genera have been described so far in the Caribbean: Caribecandona Broodbakker, 1983, Cubacandona Broodbakker, 1983and Danielocandona Broodbakker, 1983. The first two genera (with six species) have been found exclusively in Cuba and Hispaniola [73], reflecting endemicity in the older islands of the Caribbean. The genus Danielocandona was found only on the northern coast of Venezuela [73]. In the Yucatán Peninsula, stygobiont ostracode species in freshwaters have not been found yet, but in anchialine environments, two stygobiont species have been recorded, so far. These species belong to the marine order Halocypria and are likely endemic to the Peninsula. Humphreysella mexicana (Kornicker and Iliffe, 1989) was first found in cenote Mayan Blue, near Tulum in Quintana Roo [46], and was originally described as Danielopolina mexicana Kornicker & Iliffe, 1989 [46]. Similarly, Spelaeoecia mayan Kornicker and Iliffe, 1998 was first found in cenote Mayan Blue [47]. These species have been recently collected in other localities of the submerged cave system "Ponderosa" and cenotes Muknal, Bang, 27 steps, and Cristal, located along the Caribbean coast of the Yucatán Peninsula, always in marine water habitats, under the halocline [74].
The low numbers of obligated cave dwellers found in freshwaters in previous sampling campaigns and the present study on the Yucatán Peninsula, may reflect the scarcity of or difficulty in collecting these organisms with standard sampling methods, or with the current sampling effort. For instance, the two known stygobiont species of the region (H. mexicana and S. mayan), were collected during scientific cave diving expeditions using plankton nets [46,47], suggesting that this alternative collection method could increase the opportunities of finding new stygobiont species. However, freshwater ostracodes, with the difference of anchialine species, which can be planktonic, are mostly benthic or nektobenthic. Therefore, the intensification of bottom sediments and ecological niches explorations, particularly in cenotes and caves walls, may increase the opportunity of finding stygobionts. Furthermore, a detailed analysis of the soft parts of the currently collected organisms may reveal stygobiont species.

Morphological Variability of Freshwater Ostracodes in Subterranean Waters of Yucatán Peninsula
Geometric morphometrics is a technique widely used in freshwater ostracodes, in both shell and appendages, to evaluate morphological variability within species populations as an environmental effect [75,76] or to postulate new species on the basis of the integrative taxonomy paradigm [77,78]. Understanding the source of variability in ostracodes, is fundamental for their use as bioindicators of present and past environments and to gain insights into the biological evolution of a region.
In the Neotropical region, the use of geometric morphometrics in ostracodes is rare; and most studies have been used to clarify the taxonomic status of the species (e.g., Cypridopsis silvestrii (Daday, 1902)) [79] or to evaluate the influence of the environment in the phenotype (e.g., Argentocypris fontana (Graf, 1931) and Limnocythere rionegroensis Cusminsky and Whatley, 1996 [80,81].
Cytheridella ilosvayi is the neotropical species in which the use of geometric morphometrics has been more intensively applied. For this species, large-scale biogeographical patterns [49,50], taxonomic issues [50], and the influence of the environment, i.e., climate and hydrogeochemistry on the shell and appendages [82], have been tested. These studies have illustrated that C. ilosvayi displays morphological variability in both shell and appendages between populations of its distribution range (southern Mexico, Central America, the Caribbean, and South America) [48,49]. Outline is the most important valve shape attribute that varies most significantly between geographical regions [50]. Given the large and disjunct distribution areas of C. ilosvayi, valve variability has been attributed to biotic (type of reproduction, dispersal) and abiotic (water physics and chemistry) features of the host environments [49,82].
In the Yucatán Peninsula, the geometric morphometric approach effectively discriminates morphotypes. For instance, Wrozyna et al. [82], identified two morphotypes in epigean environments: one of them is attributed to C. ilosvayi, and the other to an undescribed species. In the subterranean waters of the Peninsula, similarly, our study reveals two well-defined morphotypes of C. ilosvayi, significantly differing in shape, but not in size. The CIL I morphotype, corresponding to cave environments, is distinctive as it displays poor ornamentation in the valve surface (the anterior first third of the valve is almost smooth) and a prominent brood chamber that exceeds the dorsal margin of the shell. These characteristics (attributed to their life in caves) differentiated them from the other forms of the Peninsula, and therefore, represents the third morphotype of C. ilosvayi in the region. This is relevant as at the continental scale, excluding Yucatán species, only two morphotypes have been recognized in adult specimens (Florida and Brazil-Colombia) [49,50,82]. Therefore, high morphological variability in Yucatán can be assumed. The presence of C. ilosvayi in the caves of Yucatán represents the first record of the species in subterranean environments, as previously it was documented in lakes, rivers, cenotes, springs, flood plains, and coastal areas [49,83]. For the CIL II morphotype which is widely distributed in the region, we could not establish a morphological relationship with the morphotypes identified by Wrozyna et al., [82], particularly because a full description of those morphotypes is not presented.
The low number of specimens evaluated (13 individuals) in C. ilosvayi may, however, represent a source of bias for geometric morphometrics, as outline analysis requires a number of specimens (independent data) to consider assumptions to be statistically robust. For ostracodes, however, geometric morphometric analysis has shown acceptable results with a relatively low number of specimens, such as 27 [84] and 29 [74], particularly in studies of paleoenvironments, where shell recovery can be a limiting factor. For cave and cenotes (mostly small-sized and oligotrophic), valve recovery of a single species can be complex and not always possible with standard methods. For instance, contrary to what was expected (as subterranean waters are low energy), the valves recovered in our samplings were not always in a good preservation state, and this may be related to ecological interaction as predation, which may affect valve integrity.
For the Alicenula sp. and C. vidua populations, valve shape variation between the morphotypes identified in the PCA plot, was significant. Furthermore, for both species, valve size analysis (length/height ratio), resulted in significant differences between the populations. This reveals that morphological variation is relevant in these species and its related to both valve size and shape. In comparison to their epigean relatives, the most prominent pattern is valve reduction. For Alicenula sp., epigean forms have a length range of 560 to 580 µm [54], and in our study forms of cenotes Xbatun, Dzonotila and Calcuch, the valve size of the organisms ranges from 420 to 450 µm. Specimens from cenote Kankirixché, however, were even larger than those in epigean environments, with a lengthof more than 600 µm. For darwinulid ostracods (such as Alicenula sp.), valve shape and size have been demonstrated to be highly conservative in epigean environments, even across large geographical areas [85]. In the Yucatán Peninsula, molecular evidence has highlighted that minimal changes in valve structure and size may indicate different species [54]. Currently, there are two Alicenula species described so far in Yucatán: A. yucatanensis and A. serricaudata [48]. Differentiation between these species is difficult based on morphological or morphometric attributes of valves only. Nevertheless, two well-defined groups were recovered with both the elliptic Fourier method and length-height analysis. The two groups identified were closely match, with the overall shorth and low A. yucatanensis (514 µm length, 198 µm height) and A. serricaudata (572 µm length, 210 µm height), but size discrepancies were observed in particular in A. yucatanensis. A definitive taxonomic position of Alicenula species inhabiting subterranean environments of Yucatán cannot be provided at this stage and the ecological significance of the presence/absence of a particular species, the valve size reduction or variability and the possible co-occurrence in a single system, must be carefully tested and should be accompanied by soft part analysis and studies of genetic divergences.
Similarly, for C. vidua, the variability of cenotes and cave populations with respect to epigean fauna is evident. The size of epigean individuals ranges from 550 to 650 µm [54], whereas, in subterranean waters, we recorded individuals less than 450 µm in cenote Xoch; most specimens have a length ranging between 480 and 530 µm. As C. vidua is a cosmopolitan species, the morphological features of shell and genetic diversity are variable, and the morphotypes may represent a complex interaction between genotype and environment. In the Yucatán Peninsula, three Cypridopsis species have been reported so far (Cypridopsis rhomboidea Furtos, 1936; C. niagranensis; C. vidua), with distinct characters mostly in appendages and reproductive structures [45,54]. Furthermore, for C. vidua cryptic speciation has been documented in the Peninsula with at least three genetic lineages, that cannot be fully discriminated by morphometric analysis [54]. For this species uncertainty in species identification remains, as we were unable to collect well-preserved and complete organisms and corroborate their taxonomic identity.
For Alicenula sp. and C. vidua, an integrative approach using morphological, morphometric, ecological, and molecular analysis is highly recommended to clarify the influence of the environment on shell structure and to evaluate diversity [54], particularly for their use as bioindicators in paleoenvironments.

Genetic Differentiation of Cytheridella ilosvayi in Subterranean Waters of the Yucatán Peninsula
In the Yucatán Peninsula, there are currently known two species belonging to the genus Cytheridella: C. ilosvayi and Cytheridella americana (Furtos, 1936) [45,48]. Based on morphometric differentiation, Wrozyna et al. [50] postulated the existence of a statistically differentiated morphotype of C. ilosvayi in Yucatán, but it remains uncertain whether the species is C. americana or an undescribed species. In ostracodes, valve morphology can be strongly influenced by the physical and chemical water properties of the systems. Most importantly, valve variations attributed to the environment are size, ornamentation, and shape [78,[86][87][88]. Appendages are also subject to phenotypic plasticity as a response to local or regional environmental factors [78,89]. For ostracode species, therefore, minimal changes in both hard and soft parts must be evaluated on the basis of the integrative taxonomy approach [90]. Molecular analysis has been demonstrated to be a fundamental technique to complement the morphological determinations, and the combination of them have proved to be reliable for species delimitation in the Yucatán Peninsula [54].
For the C. ilosvayi individuals of our study, it remains uncertain whether the morphological interpopulation variability is related to local environmental factors of caves and cenotes or derived from speciation processes.
Phylogenetic analyses undertaken in C. ilosvayi specimens of subterranean environments and epigean populations shows three main groupings on tree topology based on COI, and with high Bayesian support (0.9). This suggests the existence of three genetic lineages in the region. Furthermore, these lineages corresponded to well-defined habitat types, i.e., cave populations, cenotes/oligotrophic lakes, and lakes with eutrophic conditions. This specificity for habitat types may also indicate an ecological separation between lineages. COI has been widely used for aquatic species identification and has a high resolution for (cryptic) species delimitation [91,92]. This gene, therefore, suggests the existence of three independently evolved lineages in the Peninsula, associated to habitat types, that may represent different species. Differentiation between habitat types for C. ilosvayi must, however, be taken with caution as this species is considered broadly tolerant to environmental conditions, as it is widely distributed in the tropical America. A well-defined ecological differentiation of morphotypes, must be accompanied by controlled conditions and experiments such as transplantation or translocation.
The 18S gene is highly intraspecifically conservative and is also effective in specieslevel identification, even with single nucleotide polymorphisms [93,94]. Intraspecifically, 18S sequences similarities close to 100% have demonstrated to be an appropriate criterion for species delimitations in aquatic groups such as copepods [95]. For C. ilosvayi, the 18S gene sequences were identical between populations (100% similarities). For this gene, the lack of differentiation may further suggest that reproductive intermixing occurs between populations. Therefore, lineage discrimination is not supported by the 18S gene in C. ilosvayi.
The lack of concordance between the mitochondrial and nuclear markers in C. ilosvayi populations, represents a challenge to determine speciation in subterranean populations.
Concatenated COI-18S analysis revealed, however, Bayesian posterior probabilities of 0.7 for C. ilosvayi lineages; the evidence for genetic differentiation within the group is, therefore, weak. Given the lack of additional morphological evidence, particularly of soft parts, we consider that there is not enough evidence to support diversification in subterranean environments for C. ilosvayi in the Yucatán Peninsula.
Regarding the taxonomic identity of Cytheridella specimens of Yucatán subterranean waters, based purely on valve size, we assume that C. ilosvayi is the most appropriate name, as they average 1078 µm in length and 553 µm in height (n = 17), which is coincident with most records of this species in tropical America [83]. Conversely, C. americana valve sizes range 550 µm in length and 310 µm in height [45]. At this stage, we are unable to test the hypothesis of Wrozyna et al., [50], on the existence of an undescribed Cytheridella species in lakes (or cenotes) of Yucatán and to resolve this taxonomic issue. Genetic sequences of Cytheridella of other different regions of the American continent will be necessary.

Late Pleistocene and Early Holocene Climate as the Driver of Subterranean Biodiversity in the Yucatán Peninsula
Aquatic environments of the Yucatán Peninsula during the Late Pleistocene and Early Holocene have been influenced by regional and local events, mainly related to climate fluctuation. From a regional perspective, environments of the Peninsula and, particularly, subterranean systems have been strongly influenced by the relative position of the sea level. Evidences from multiple biological and non-biological indicators, such as corals and speleothems, suggest that during the Late Pleistocene and, particularly, the Last Glacial Maximum (LGM), the marine water level in the Caribbean decreased bỹ 120 m, relative to today's levels [96,97]. This suggests that the water level in cenotes, was much lower than today. Evidence from sedimentary sequences of Lake Petén Itzá, the largest and oldest system in the Yucatán Peninsula, supports that events such as Heinrich Stadials and deglacial were characterized by precipitation decreases and periods of severe droughts [98,99].
In terms of freshwater availability, Lake Petén Itzá decreased during the Heinrich stadial 1 by about 50 m with respect to their current water level, with a reduction of about 80% of the total area [100][101][102][103]. This estimation allows us to postulate that water systems in the entire Yucatán Peninsula were severely affected by similar magnitudes as Lake Petén Itzá [98,99,101]. Few systems must have remained filled during the dry spells of Late Pleistocene-Early Holocene, with water chemistry that was likely saturated with minerals as a consequence of evaporation rates.
For subterranean waters, the lowering of the sea level and the prevailing inland dry conditions, might reduce subterranean water availability and caused modifications of the chemical composition of the systems by rock weathering, i.e., increased conductivity by calcium and magnesium incorporation. These conditions may have triggered the reduction of groundwater sources and the millennial-lasting complete isolation of subterranean environments. Long-lasting systems isolation, together with environmental variability, is recognized as a driver for biological diversification worldwide [104,105]. Currently, in the Yucatán Peninsula, a high number of freshwater endemic and stygobiont species have been recorded so far, in groups such as Copepoda, Thermosbaenacea, Stygiomysida, Mysida, Isopoda, Amphipoda, and Decapoda [26,74,106,107], all with well-developed capabilities to inhabit temporarily or permanently the water column.
The low number of freshwater endemic and stygobiont species in Ostracoda of subterranean waters of Yucatán, can be related to the ecological behavior of the group and the environmental fluctuations of the late Pleistocene and Holocene [108]. For instance, ostracodes can be considered benthic species adapted to exploit sediments, and with limited capacities of swimming in the water column.
Sedimentary sequences of the Yucatán Peninsula have revealed that the sea-level rise was progressive after the drastic decrease of the Late Pleistocene. During the Holocene, bottom sediments were subject to strong environmental alterations given sea level fluctuations. Data from the currently submerged speleothem of the Yucatán Peninsula have revealed periods of either increased/decreased water level during 8.9, 8.6, 8.4, and 6.0 ka BP (before present) [109]. The frequency of bottom sediments instability in subterranean environments might have represented an important limiting factor for benthic species success, i.e., adaptation and diversification.
In cenote Aktun Ha, located~8 km away from the Caribbean sea coast, the fossil record of C. ilosvayi and D. stevensoni during the last~6.8 ka, has revealed phases of relatively high ostracode abundance during low water stands (relative to lower sea levels). The increase of water levels following a sea-level rise of about 4 m at 6.5 ka BP, drove the extinction of C. ilosvayi in Aktun Ha, likely due to the effect of marine water intrusion in the subterranean aquifer [110]. Conversely, D. stevensoni, remained present in the system despite hydrological alterations [110]. The ecological response of C. ilosvayi and D. stevensoni in cenote Aktun Ha may illustrate the regional effect of sea level variability in benthic assemblages on a regional scale, with some species more susceptible to subterranean environmental change and others with better capacity to tolerate them.
In isolated sites such as Yumku and Peba, which are systems not directly connected to the subterranean aquifer, C. ilosvayi demonstrated morphological variability in carapace structure, as well as a clear genetic differentiation compared to their relatives in cenotes and lakes. This may be indicating different evolutionary pathways for species of subterranean environments in which marine intrusion has shown none or little relevance in the geological history. These sites are, therefore, highly important for a comprehensive evaluation of the biological evolution of subterranean species associated with sediments in the region. The exclusive occurrence of epigean species in freshwaters in cenotes and caves of Yucatán may, therefore, be reflecting the ecological alteration of the Late Pleistocene and Early Holocene climate variability in the region [111].
Colonization of epigean species is a common trend in subterranean environments worldwide, and in most systems, a combination of epigean and hypogean fauna can be observed [8,68,112]. Evidence from different aquatic groups reveals that colonization of subterranean environments is explained by active and passive models [31,32,113]. Active models describe that: (i) the species found refugium from climate variations (the climaterelict model), or strong ecological competitive pressure; (ii) generalist species colonized subterranean environments even in the absence of climate or ecological stressors (the active colonization model), and (iii) marine species invaded land through subterranean connections [31,113]. Passive models suggest that colonization occurred by transportation, either by biotic or abiotic factors [31]. The evidence observed from the freshwater species assemblages of subterranean environments in the Yucatán Peninsula is, by far, not definitive; given the combination of a high percentage of generalist species and epigean regional endemics and the practical absence of microendemics and obligate cave-dwellers, we can infer that active and relatively recent colonizations (Holocene) may have taken place in the Yucatán Peninsula. Such colonization may be occurring at different time scales (with different lineages) and may be occurring polytopically, as it has occurred in subterranean zooplankton assemblages worldwide [32].

Conclusions
The Yucatán Peninsula is a karst region with high availability of subterranean waters and hypogean biodiversity. The dominance of epigean species and low freshwater ostracode species richness characterize the subterranean waters of this region. About 40% of the species are generalists widely distributed in the American continent, and the remaining 25% are regional endemics. All other species remain in open nomenclature, given the limitations in the taxonomic analysis based on soft part analysis. The absence of stygobiont species in freshwaters, can be related to the applied sampling methods or effort in the field. For instance, in this study, we used scuba diving, sediment analysis, and water filtration, but the increased use of sampling and analysis of bottom sediments and cenotes walls, may represent an appropriate method for hypogean species recovery, given the benthic habits of the ostracodes. Morphometric analyses of valve shape based on the elliptic Fourier analysis and valve size using length-height determination, allowed us to capture the morphological variation between populations. The three species evaluated demonstrated morphological differences between populations. Cave (epigean) populations display the highest morphological variation (shape and size) with respect to species of cenotes. Molecular analysis is a powerful tool for understanding interpopulation variability in C. ilosvayi. This species displayed three molecular lineages within subterranean and epigean habitats, but we lack conclusive evidence of species diversification in caves and cenotes. The climate history of the region is assumed to be the main driver for species composition. Subterranean marine intrusions, as well as severe droughts that occurred during the Late Pleistocene and Early Holocene, likely modified water and sediment chemistry, which may have had a direct effect on benthic species composition.
Supplementary Materials: The following are available online at https://www.mdpi.com/1424-2 818/13/2/44/s1, Table S1: Location data of the studied cenotes and caves, Table S2: GenBank identification codes and accession numbers for COI and 18S rDNA sequences of Cytheridella ilosvayi populations, Figure S1: Phylogenetic tree of concatenated COI-18S dataset on Cytheridella ilosvayi populations of cenotes, caves and lakes in the Yucatán Peninsula. COI tree topology based on Bayesian Inference (BI); numbers below branches represent Bayesian posterior probabilities; abbreviations on tree topology correspond to GenBank species identification codes (see Table S2).

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to their usage in the ongoing study. Sequences of this study are available from GenBank (accession numbers MW018845-MW018862 for COI and MW018806-MW018816 for 18S rDNA).