Dna Barcodes, a Powerful Tool for a Rapid Construction of a Baseline and Conservation of Aquatic Ecosystems in Sian Ka’an Reserve (Quintana Roo State, Mexico) And Adjacent Areas

This study is focused to the aquatic environments of the Sian Ka’an reserve, a World Heritage Site. We applied protocols recently developed for the rapid assessment of most animal taxa inhabiting any freshwater system, by using light traps, and DNA barcodes, represented by the mitochondrial gene Cytochrome Oxidase I (COI). We DNA barcoded 1037 specimens of mites, crustaceans, insects, and fish larvae from 13 aquatic environments close or inside the reserve, with a success rate of 99.8%. In total, 167 Molecular Operational Taxonomic Units (MOTU’s) were detected. From them we identified 43 species. All others remain as a MOTU. For analyzing the adult fish communities, we applied the non-invasive method of environmental DNA (eDNA), and identified the sequences obtained with the Barcode of Life Database (BOLD). We found 25 fish species, and other terrestrial vertebrates from this region. No alien species was found. After a comparison of the MOTU’s from all systems, we found that each water body was unique respect the communities observed. The reference library presented here represents the first step for future programs to detect any change in these ecosystems, including invasive species, or improve knowledge of freshwater zooplankton, because most of the MOTU’s are possibly new species to science.


Introduction
Freshwater is the most valued resource globally, representing around 1% of all water in the biosphere. However, pollution, excess use, and alterations due to global warming and dry seasons are changing its availability, affecting all terrestrial life forms that depend on it for their survival. As a result, all the animal communities dwelling in these freshwater ecosystems are highly vulnerable, especially to the invasions of alien species. Mexico is not an exception [1] to this situation. Hence, we chose to study an important system and its surroundings, located on the east coast of the Yucatan Peninsula: the Sian Ka'an biosphere reserve (Fig. 1).
First and foremost, it is essential to point out that the inventories of freshwater species are still in an incipient phase of knowledge since few groups of specialists in the matter are scattered throughout the world. None can identify all individuals present in any particular system to the species level, making it even more relevant to join efforts. For example, back in 2005, an estimation of the diversity of branchiopods, a crustacean group dominated by freshwater species, was c. 1180 species in the world; this figure was also reported to be at least two times more [2]. Since then, about 156 species of branchiopods have been described formally (Web of Science, May 3, 2021; Search: "new species" and "branchiopod"); Adamowicz and Purvis concluded that the neotropics are one of the regions where more species remain to be discovered [2] suggesting the importance of its study. Some recent developments helped to understand the biodiversity in these particular environments through DNA barcoding, making it a powerful tool to study aquatic organisms from invertebrates [3] to vertebrates [4]. Therefore, the construction of databases with this information is fundamental, being the Barcode of Life Database (BOLD, boldsystems.org) one of the best to use because it focuses on biodiversity exclusively [5].
With the popularity of metabarcoding after the development of new generation sequencers, many labs focused on searching targeted species in aquatic environments, as aliens [6] or endangered [7]. Others focused on a whole community, like fish, but these latter studies rely upon a public database of DNA barcodes [8]. The relevance of the database has been highlighted for many groups as the insects before the metabarcoding, even though the detailed curation of names is incomplete [9].
Hence, the main problem for any aquatic community study is the lack of knowledge of the species (the so-called taxonomic impediment) and the consequently limited information on the DNA analysis.
For several years now, our group has been working on constructing a public database with the most accurate information for DNA barcodes. This effort has been an example for others [10]; nonetheless, a proper assessment of the biodiversity within Mexico remains distant, mainly because of the global taxonomy crisis associated with the fact that many new barcodes have not been matched with their corresponding species. To overpass this taxonomic impediment has particular relevance in the neotropics, the most speciesrich region in the world [11]. Consequently, we consider the best practice is to document the specimens studied with vouchers deposited in museums. If their DNA is barcoded, there are several approaches to establish them as molecular taxonomic units (MOTUs) that can be a hypothesis to be tested against a real species identification [12]. BOLD provides a useful calculation of these MOTUs, the Barcode Index Numbers (BINs) [13] were created as part of a rapid solution to this taxonomic impediment.
In the case of the Mexican continental territory, the Yucatan Peninsula, set in the neotropical region, has a complex system of underground freshwater that emerges to the surface as sinkholes (locally known as "cenotes"), lagoons, and "aguadas," mostly shallow rounded water bodies [14]. Its east coast hosts the second largest coral reef in the world, known as the Mesoamerican Reef. Unfortunately, tourism has been growing with no order along this coastline, from Cancun to Chetumal city (near the border with Belize), including the Sian Ka'an biosphere reserve. This expansion shows no respect to natural resources.
As a result, the freshwater systems in the region present pollution of hydrocarbons going from Cancun (in the north) to Bacalar Lake and Milagros lagoons, located 10 km from Chetumal city [15]. This situation is relevant because the physical and biological interactions between mainland freshwater and the sea are not well understood due to the lack of studies. For example, just recently, thanks to the finding of larvae and juveniles of the fish Cyprinodon artifrons (Hubbs, 1936), which adult stages are found around the mesoamerican reef, were discovered to had migrated to breed into the freshwater system of Bacalar Lake, located more than 70 km from the ocean [16]. We do not know how these larvae reached the lake because adults have not been found here.
Sian Ka'an, which means "Origin of the Sky" in the Mayan language, is a Natural Protected Area with nearly 4,000 km2 of land surface and continental waters (Fig. 1). It has an outstanding conservation status on its hydrology and ecosystems, including tropical forests, mangroves, wetlands, beaches, sinkholes, and marshes (Fig. 2). It was designated as a biosphere reserve of international importance by UNESCO in 1986 and recognized as a World Heritage Site in 1987 [17]. In 2003, it was recognized as a wetland of importance by the Ramsar Convention [18]. Nevertheless, notwithstanding its importance, there is surprisingly limited information regarding the freshwater species inhabiting the aquatic habitats here. The existing records include more than 40 fish species, with only two alien species: the Mozambique tilapia, Oreochromis mossambicus, reported long ago from an "aguada" near the reserve limits [19] and the Nile tilapia, O. niloticus. Among aquatic invertebrates, recent surveys in these environments allowed the description of one calanoid copepod, Mastigodiaptomus siankaanensis [20], and the harpacticoid Remaneicaris siankaan [21]. Nevertheless, after these studies, there is no other formal description of the aquatic biota from Sian Ka'an reserve, except non-peer-reviewed lists provided by Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) [22], and Naturalista [23]. There is another unverified list published by Comisión Nacional de Áreas Naturales Protegidas (CONANP, Mexico) [24]. Finally, a report from United Nations Development Programme [25] listed 80 operational taxonomic units (OTUs), but the authors warned that this report was preliminary.
Recently, Elías-Gutiérrez et al. [16] and Montes-Ortiz & Elías-Gutiérrez [26] showed the usefulness of light traps combined with DNA barcoding to recognize the zooplankton biodiversity despite the little taxonomical knowledge in some temperate and tropical water systems on the Yucatan Peninsula.
We present here how to establish a new and rapid baseline of all possible BINs or potential freshwater species found, including invasive alien species, if present, around Sian Ka'an reserve and nearby aquatic systems. It will be based on sequencing of fragment of the COI gene (for zooplankton) and metabarcoding methodologies [8], using the eDNA for fish as an alternative to collect the specimens directly. Finally, we will discuss how both results could be used in the future as a biomonitoring tool.

Sampling
Samples were collected from 13 water bodies near and around the buffer zone of the Sian Ka'an reserve (Table 1, Fig. 1) on August 21-26, 2019. We decided to work within the reserve, including other adjacent water bodies, because we know that all these systems are not isolated. They can be connected by subterranean waters or on the surface after the floods, a common phenomenon in this region.
Laguna Muyil and Chunyaxché were the only systems sampled in two sites. Light traps with 50 µm mesh [26] were deployed overnight in the limnetic and (or only) littoral zones (the latter designed with a number with the same name of the locality), depending on the size of the water body. The year 2019 was uncommonly dry, during August (supposedly part of the rain season, occurring from May to October). This phenomenon possibly triggered forest fires. Consequently, these fires limited the possibility of visiting more systems (also unusually dry) [26].  After collection, water from all the samples was extracted with a 50 µm sieve by washing them with cold (4°C) 96% ethanol. Specimens were then transferred into jars with approximately 1/3 of the sample and 2/3 of ethanol [16]. Next, the sample jars were placed in a container with ice for being transferred to the laboratory, where they were stored in a freezer at -18°C for at least one week. After this period, samples were kept at room temperature.
For metabarcoding from each system, we collected 1 L of water from the sub-surface in a sterilized CIVEQ bottle as suggested for the nearby Bacalar Lake [20]. Water samples were filtered in a field lab in Carrillo Puerto town, to minimize eDNA degradation. We filtered at least 250 ml of water from each point in 0.22 µm filters. All filters were stored in cold gels till further analysis (see 2.4. Metabarcoding and eDNA section for more details).

DNA extraction and amplification of individuals
All specimens collected were sorted out under a stereomicroscope from the zooplankton samples; representatives of each morphologically distinct taxon were then photographed using the same microscope. After being photographed, specimens were again placed in ethanol and stored at room temperature.
Small specimens (less than one mm) were destructively analyzed. DNA from larger specimens was obtained from an antenna, eggs or embryos, or a piece of the abdomen, preserving the remaining parts, like the head and the end of the abdomen (e.g., chironomid pupae). In the case of water mites, voucher specimens were recovered as suggested for Collembola [27] and preserved in 96% ethanol with a drop of glycerol. Finally, for fish larvae, one eye of the right side or a piece of muscle was used. The vouchers (specimens not lost during extraction) were deposited in the Reference Collection at El Colegio de la Frontera Sur (ECOSUR), Unidad Chetumal. All specimens were identified to the finest possible taxonomical level by using general identification keys [28] or detailed descriptions, as M. siankaanensis, for example [20], or by comparison with previous DNA barcodes in the Barcode of Life Database (BOLD, bodsystems.org).

Sequencing and data analysis
PCR products were cycle sequenced using a modified [30] BigDye© Terminator v.3.1 Cycle Sequencing Kit (Applied Biosystems, Inc.) and sequenced bi-directionally in a sequencing facility (Eurofins) using M13F and M13R primers. Sequences were edited using CodonCode v. 9.0.1 (CodonCode Corporation, Dedham, Massachusetts) and uploaded to BOLD, and are available in the dataset Baseline Sian Ka'an I (DS-BASKAAN; dx.doi.org/10.5883/DS-BASKAAN). All data were analyzed with the quality tools on BOLD, and all sequences were examined for the presence of stop codons and indels as a check against NUMTS.
A Neighbor-Joining (NJ) tree was calculated by using the Maximum Likelihood method based on the Kimura two-parameter (K2P) distance model [31] with 500 Bootstrap replications including all taxa found in the four major groups: Arachnida, Crustacea, Insecta, and Actinopterygii found in all systems, using the MEGA 7 software [32]. Each tree was simplified with the Compress feature of Mega 7. Additionally, we prepared a general ID tree with all specimens, with the analytical tools provided by BOLD (Supplementary File 1). We selected this method because it allows the rapid analysis of large data sets of specimens [33]. MOTUs as a proxy to species were delimited with the Barcode Index Number (BIN) [13] that has proven more than 80% of effectiveness [16] and has been widely accepted [9]. Based on the BINs (=MOTUs), we prepared a list to the finest possible identification of all species found. Mexico has one of the most extensive datasets in public databases of freshwater species in the world [10].

Metabarcoding and eDNA
Water filters from each sampling point were sent to be processed in the Centre for Biodiversity Genomics in the University of Guelph (Canada). The interval between filtration and DNA extraction was less than 48 hrs.
Each sample tube was kept at -20°C before processing. Before DNA extraction, all lab surfaces and pipettors were sterilized using 10% bleach, followed by 70% ethanol [8]. Prior to extraction, each filter was placed in DNEasy PowerWater Bead Tube (forceps were sterilized with 20% bleach, followed by 100% ethanol and triple-flame sterilization between samples). DNA was extracted with minor modifications to previous publications [34]: a volume of 900 µL of ILB buffer with 100 µL Proteinase K was added to each tube, tubes were incubated at 56°C for 30 min and vortexed on Genie 2 vortex at max speed for 5 min, followed by 1.5-hour incubation at 56°C. Tubes were centrifuged at 2,000 g for 2 min and ~700 µL of lysate transferred to a clean tube containing 1.4 ml of 5M GuSCN buffer; all resulting volume was applied to Epoch Biolabs column in three subsequent centrifugation steps at 6,000 g (each 700 µL transfer). Silica membrane was washed 1x with 400 µL of PWB, 2x with 700 µL of WB, dried at 56 °C for 20 min. DNA was eluted in 100 µL of EB buffer at 11,000 g. DNA was transferred into a 96-well plate in 3 replicates per sample. After extraction, we followed a two-step PCR approach with the first round employing conventional primers, while the diluted PCR product served as a template for a second round of PCR with fusion primers containing sequencing adapters and UMI-tags. A 184-187 bp segment of the barcode region of COI was amplified with two primer sets (Aq-uaF2/C_FishR1, AquaF3/C_FishR1). The PCR reactions employed the master mix described previously [35] and Platinum Taq. The first round of PCR employed the following thermal regime: 94˚C for 2 min, 40 cycles of 94˚C for 40 s, annealing at 51˚C for AquaF2 or 50˚C for AquaF3 for 1 min, and 72˚C for 1 min, with a final extension at 72˚C for 5 min. PCR´s resulting products were diluted 2x and used for second round PCR with fusion primers for 20 cycles to create IonXpress MID-tag labeled libraries. The PCR regime for the second round consisted of 94˚C for 2 min, 20 cycles of 94˚C for 40 s, annealing at 51˚C for 1 min, and 72˚C for 1 min, with a final extension at 72˚C for 5 min. PCR products were visualized on an E-Gel96 (Invitrogen, Thermo Fisher Scientific).
The library for each unique primer pair was pooled without normalization and purified using magnetic beads with a 0.5 to 1 bead to product ratio for the uppercut and 0.8 to 1 ratio for the lower cut. Each library was diluted to 26 pM and mixed in a 1:1 ratio for the S5 run using Ion 510/520/530 chef kit. Each DNA replicate for each primer combination was processed under a separate IonXpress MID-tag.
The following procedure to process the raw NGS reads: Cutadapt (v1.8.1) was used to trim primer sequences; Sickle (v1.33) was used for size filtering (sequences from 150 -250 bp were retained), while Uclust (v1.2.22q) served to recognize OTUs based on the >98% identity and a minimum read depth of 2 thresholds. The Local Blast 2.2.29+ algorithm was then used to compare each OTU to the reference sequences in five datasets: public fish data from BOLD filtered to genus and species ID level (130,357 sequences), and public BOLD data for Amphibia, Aves, Mammalia, and Reptilia represented by the following datasets: DS-EBACAMPH (11,018 sequences), DS-EBACAVES (28,914 sequences), DS-EBACMAMM (39,890 sequences), DS-EBACREPT (5,424 sequences). Raw Blast output results were parsed using custom-built Python scripts. Processed results in tab-delimited format were imported to MS Excel, filtered by a minimum score of 250 and 97-100% percent identity range. Blast search results were parsed and concatenated using custom-built Python scripts, exported to Excel, and visualized using Tableau software.

Statistical analyses
Since data are absence-presence, variation species composition between sampling places was performed using the Sørensen index of dissimilarity [37][38][39]. All analyses were performed with Betapart Package in R software [36].

Results and discussion
All aquatic ecosystems studied here, except Tres Reyes 1, are of transparent blue water (see Fig. 2), accordingly with most of the aquatic systems in the Yucatan Peninsula, which are oligotrophic, with an average Secchi disk of 7.6 m [37] and exhibit a predominant presence of diatoms [38].

DNA barcoding baseline
1,037 specimens were processed in total for DNA barcodes. We had success for 1035, corresponding to a 99.8%. These results are remarkable because we used for all groups presented here only a single pair of primers, Zplank forward and reverse [39]. We followed all recommendations by Elías-Gutiérrez et al. [16] about the fixation of the material with cold ethanol and preservation in cold storage the first week. We consider these procedures contributed to the high success obtained here. All sequences and chromatograms passed the quality controls of BOLD, and none was contaminated, or with stop codons, only some few sequences were shorter than 500 bp (see dataset DS_BASKAAN, Baseline Sian Ka'an dx.doi.org/10.5883/DS-BASKAAN in BOLD).
A total of 43 species were identified to species level, and they are coincident with the BINs assigned by BOLD (Figs. 3-6), with less than 2% of intraspecific divergence. All others remained as a MOTU (=BIN), a putative species.
The list of species and potential species, represented by the BINs presented here, is the largest published freshwater biota of Sian Ka'an reserve with 167 (Supplementary File 1). It includes many groups not considered as "true" zooplankton, such as the water mites or chironomids. However, Montes-Ortiz & Elías-Gutiérrez [25] discussed the presence of these groups in the zooplankton community.
We found that each system has a unique assemblage of species, and Muyil and Chunyaxché have a unique distribution of the species, which will be discussed later.
Arthropods were represented by water mites (Arachnida) with 209 specimens comprising 11 families and two orders, with 40 BINs.
The high diversity of mites in this region was noticed in two previous studies [16,40]. Their relevance as indicators of water quality is well known [41]. Their taxonomy is yet unresolved, and they have a limited distribution in many cases, with possible endemics. Some of them as representatives of Limnesia (ACY7380), Unionicola (AEB1594), or some Pionidae (AEA4808), and Eylaidae (AEA5669 and AEA4696) were present in only one system (Fig. 3) [40]. . Simplified ID tree for Arachnida represented in the samples. All of them belong to Trombidiformes. Some of the sequences presented here have been published before (18). Numbers represent the aquatic system, as presented in Table  1. The last number corresponds to the BIN assigned in BOLD. Support bootstrap values are on each branch. Crustaceans were represented by 436 specimens comprising five classes, five orders, and 13 families, corresponding to 50 BIN's, 20 of them allowing identification to species. Other seven were previously identified as possible comparable species but possibly different (cf.). Finally, seven more were identified to genus, and all remaining can be considered putative species after the BIN.
The Cladocerans were represented by 30% of the total crustaceans and usually were not dominant in the samples. Although scarce, Sminov and Elías-Gutiérrez [38] found eight cladoceran genera after a survey in sediments from 25 systems in Yucatan. For example, recent studies on 56 systems from Mexico and Central America [42,43] found only three taxa in five out of 19 systems located in the Yucatan Peninsula. These numbers indicate, as Smirnov & Elías-Gutiérrez [38] suggested, that cladocerans´ preference for littoral areas and their rapid decomposition in the bottom mud difficult further preservation. From the 15 branchiopod BINs found, only four could be identified with certainty, confirming the need to improve the taxonomy for this group [3,44].
It was surprising to find that the sinkhole Tres Reyes 2, with blue and clear water (number 8 in Fig. 4), is highly dominated by a species of the daphniid Ceriodaphnia, with more than 90% of the total zooplankton. Named here Ceriodaphnia cf. rigaudi, it is possibly an undescribed species with previous records in Yucatan [3]. All the other systems dominated copepods (calanoids and cyclopoids). Tres Reyes 2 is quite remarkable in a region were cladocera seem to have vanished as the dominant group. We cannot explain the reason for this phenomenon, but the system where Ceriodaphnia cf. rigaudi dominates seems to be a deep oligotrophic sinkhole with vertical walls. This particular system requires further studies, but access to it is pretty difficult.  Table  1. The last number corresponds to the BIN assigned in BOLD. Bootstrap support values are shown in each branch.
The ostracods were found in ten systems (Fig. 4). They did not appear in all systems, probably due that samplings were conducted only with light traps. It is necessary to compare results using different devices as plankton tows to clarify if this group was underrepresented here. Although recently several publications included this group [45][46][47], we could identify only Darwinula stevensoni to species level. Ostracods are an important group in Yucatan Peninsula, more abundant in sediments than cladocerans [38], with great value as paleobioindicators [46]. Although the study of the diversity of this group is just starting here, a new genus was described recently here [45].
The most common and dominant members of the planktonic community in most systems were the copepods, represented by two groups: Calanoida and Cyclopoida. Among the identified species, only the copepod Mesocyclops edax is distributed from Yucatan to Canada. We cannot establish if it was introduced or not. Mexican previous records of this species are only from Yucatan [28]. Even though some of the taxa registered here are restricted to Yucatan or the south of Mexico, others have a more global range, such as Tropocyclops prasinus and Macrocyclops albidus are found from the southern lowlands to highlands of the central plateau of Mexico [28]. However, these two species seem to be a complex of a sibling or translocated species [28,48], requiring a detailed study.
The analyses of these animals in Yucatan allowed the discovery of new species, possibly endemic in sinkholes [49] where cyclopoids have been barcoded since 2008 [3], allowing the identification of ten species of the 24 BIN's found. In the case of calanoids, we identified the four BINs found to species. However, two could be cryptic, named Arctodiaptomus cf. dorsalis and A. cf. dorsalis 1 [3]. Calanoids were absent in five systems (Km 48, Tres Reyes 2, Santa Teresa, El Toro, and Pucté 2), but cyclopoids were present in all. When calanoids were present, they were usually the most common group in the samples.
From the six BINs representing Decapoda, we identified only Palaemonetes intermedius. This caridean shrimp that can penetrate from the sea to freshwater systems, widely distributed on the Western Atlantic coast of the continent [50]. Most decapods, mainly represented by zoea larvae, were found in specific sites of Muyil and Chunyaxché lagoons, both with direct connection to the Caribbean, located 10 km away in a straight line, through a series of channels.
In the case of insects, from the 50 BINS found, Chironomids were the most important group present in all samples (Fig. 5). We also found three other three orders and four families. In total, 44 BIN's were chironomids, and all of them need further studies since most of this fauna is not well known in the study area. Three BINS were shared in three clusters, but the number of sequenced specimens is still low (two of them are singletons), and this phenomenon has been noticed before in aquatic insects [51]. Chironomids were represented basically by pupae, but some adults emerged inside the trap. Recently, in a study about subfossil Chironomidae present in sediments from Yucatan Peninsula, the maximum resolution was to groups of species, and the authors concluded that due to poor sedimentation and preservation of remains, cenotes have limited potential for palaeolimnological studies [37]. This conclusion highlights the urgency to bio-assess this community at the present day. Our material was sequenced, but taxonomical parts (head and tail) are preserved in the zooplankton reference collection at ECOSUR Chetumal, allowing future identification of the BINs. The only detailed recent study on this fauna is the description of six new species from several systems of Yucatan [52]. Previously, the same author pointed out that this region remains unknown regarding chironomid fauna [53]. In the list presented by these latter authors, they identified only 13 species (most of them represented by adults) from a total list of 86 potential taxa. Another small report lists 42 genera from Calakmul reserve in the south of the Yucatan Peninsula [54], with no comments. Other insects less common but important in some systems were the larvae of Chaoboridae, mostly restricted to the three southernmost areas (Fig. 5). These were represented by two BINs. It is remarkable the presence of two Ceratopogonidae BINs because most of the larvae for these tiny biting midges are known to be associated with wetlands but have not been reported as part of the "true" zooplankton from lakes or sinkholes. In this case, the larvae swim into the trap, attracted by the light, although they do not seem to be common. The only findings of them were at Muyil Lake and El Toro sinkhole. From the whole group, only two BINs could be identified to species level: the common chironomid Cladopelma forcipes and the odonate Argia translata. All others need further studies to assign the correct species name.   Table 1. The last number corresponds to the BIN assigned in BOLD. Bootstrap support is shown in each branch. Figure 6. Simplified tree for Actinopterygii. Numbers represent the aquatic system, as presented in Table 1. The last number corresponds to the BIN assigned in BOLD. Bootstrap support is shown in each branch.

Preprints
Finally, chordates collected with light traps were represented by 94 fish (Fig. 6), mostly larvae. They included six orders, seven families, 11 genus, and 12 species with a total of 12 BINs. All fish species collected with the light trap were also detected with metabarcoding. The latter method added 15 species more ( Table 2).
The fish have a complete taxonomy work, and the database for them in Yucatan Peninsula started to be built in 2005 [4]. All BINs from the light traps could be identified to species. Only one of them, with a clear, unique haplotype (Atherinella), still needs clarification about its identity. This species was found in three sinkholes (Pucté Cafetal 1, Pucté 2, and Sijil Noh Há), but records from BOLD indicate a distribution of this BIN (AAI4788) in the southeast of Mexico. All other species are well known in this region. The rest of the fish species we found have been previously reported in Yucatan Peninsula.

Metabarcoding and eDNA
Metabarcoding was focused only on vertebrates, particularly fish. Negative controls did not produce any visible PCR products and meaningful data. Human DNA was detected and filtered from the final results. Only high-quality reads assigned to correct Ion Express MID tags were used for the analysis. In total, 172,244 valid sequences were obtained from a total of 17,944,836 reads. They are summarized in Table 2. In total, 25 fish species, two reptiles, five mammals, and three birds were found. All of them matched our previous baseline. From the total, one slider (Trachemys sp.) and two fishes (Bramocharax-Astyanax and Cyprinodon beltrani-simus complexes) had a low interspecific resolution that has been highlighted previously [8]. The site with more species was Muyil 1 with 12 species. Meanwhile, Sijil Noh Há had the lowest number (two species). All systems were positive for eDNA, and all these species had previous records here or in the nearby systems [55]. Not any alien species were found among the fish. A previous old record of the alien tilapia (Oreochromis mossambicus) [19] in Chancah Veracruz (=Yodzonot in the old record) was not found. However, this species was found in a previous eDNA survey of Bacalar Lake, located about 70 km to the south of the Sian Ka'an reserve [8]. All species of larvae found in the light traps were confirmed with metabarcoding.

Species composition comparison
Sørensen index of dissimilarity indicates high differences among the BINs assemblages of sampled places. They reached a value of 0.87  0.006 SD. Based on the pair-topair comparison, the cluster analyses reflect a high dissimilarity in the BINs assemblages among sampled sites (Fig. 7). Sørensen dissimilarity plot indicates lower dissimilarities between Chunyaxché 1 and Chunyaxché 2 (0.43), between Chunyaxché 2 and Muyil 2 (0.49), and between Muyil 1 and Chunyaxché 1 (0.55), the remaining dissimilarity values between pairs are higher than 0.60 (Fig. 7). Chunyaxché and El Muyil have well-differentiated zones, with a different community in the eastern localities (Muyil 1, Chunyaxché 2, see Fig. 8). However, these two systems are permanently connected by an ancient channel possibly built by the Mayan culture. They also erected an important ceremony center near the shore of the lagoon with the same name currently. Then the two lagoons keep a similar assemblage of species, as we can observe after Sørensen (Fig. 7). Minicenote, in the southern limit of the reserve, hosts a unique assemblage of species. The peculiarity of this system [14] allowed describing a particular pattern of migration of the Mastigodiaptomus nesus population here to the walls [56]. A later, more detailed study of Minicenote concluded it as an oligotrophic system inhabited by 18 taxa, 13 of them rotifers and five crustaceans, and not any chironomids [57]. Our study confirmed the presence of a bosminid in this sinkhole, named by them Bosmina hagmanni, and that DNA barcodes allowed to distinguish it as a confusing morpho of Bosmina tubicen present in several systems of the Yucatan Peninsula [44]. The presence of M. nesus and Thermocyclops inversus was also confirmed, but not any rotifer. It seems most rotifers are not attracted by the light of the traps, although some have been reported [16].
Though several BIN's are shared among all systems, the Sørensen analysis shows that each system has a unique assemblage of species. This singularity for each aquatic system makes evident the fragility of these ecosystems in the whole region and confirms previous analyses with mites [40]. We do not know why each system has a unique assemblage of species if they eventually can be connected in the surface after floods or through the complex underground system of rivers. For example, M. siankaanensis was found only in two sinkholes (Tres Reyes 2, co-existing with the dominant C. cf. rigaudi and Pucté Cafetal) meanwhile, Mastigodiaptomus nesus was more widespread and common in the samples. The original localities of the description for M. siankanensis near Vigia Chico (Aguada Vigia Chico, Savannah 2 and Aguada Limite de la Reserva, the type locality) [20] were dry during our survey, as was previously mentioned. We do not know anything about the biology of any zooplankter in this region. However, it was recently detected an apparent genetic structure in the haplotypes of M. siankaanensis from the center to the south of Yucatan [58], being different haplotypes and possibly in the process of differentiation.
From the origin of the water and possible water flows derived from remote sensing and fault structure, it is suggested that there is a fault from NE-SW communicating systems from Sian Ka'an to the Hondo River (the border between Mexico and Belize) [59,60]. Nevertheless, the biological communities of zooplankton seem to respond more to local variables that could be predation or morphometry of the systems [14] and other unknown factors. This uniqueness of each system or group of systems, as Muyil and Chunyaxché, highlights their vulnerability to environmental issues and the urgency of a permanent biomonitoring system based on the knowledge of most species dwelling in each water body. This study could be considered the basis to understand any posterior change in these water systems.

General remarks and future biomonitoring
With the baseline provided here, we have a reliable source of information to apply next-generation techniques connected to eDNA. We demonstrated the usefulness of the DNA barcodes database to identify all species with non-invasive methods with the fish.
In the different freshwater zooplankton groups, there is no laboratory with the ability to identify all the species found in any given aquatic ecosystem. This situation is even worse in megadiverse countries such as Mexico that is the fourth country with the highest biological diversity in the world [61]. Although most of the BINs found with this study have no formal name to species, all material was deposited and photographed. With time, some specialists will collaborate to identify this material.
Though the work to document Mexican freshwater zooplankton with DNA barcodes has been an example for several years [10], most species remain not well known. The use of new collecting methods as the light traps increased dramatically the species found in neotropical systems of Mexico [16,26]. With the DNA barcode baselines, we can apply the metabarcoding immediately. These methods will allow the detection of any change in the aquatic community if a periodic biomonitoring system is established. These methods should have the same efficiency with bulk samples of the zooplankton (Elías- Gutiérrez, unpublished).
The subsequent phases should be continuing with the baseline to expand it to more water systems and seasons of the year. The knowledge of the aquatic diversity will help calibrate metabarcoding methodologies to detect any environmental threat from alien species introduction (carrying changes in predation, for example) to pollution. It has been demonstrated that zooplankton responds rapidly to these stressors [62], making it an excellent element to assess any change in the environment [63].

Conclusions
We currently can construct a fast baseline of the species in any neotropical aquatic environment despite the taxonomic impediment. If the species are unnamed, they can be kept in a scientific collection; meanwhile, the genetic data for their identification should be stored in a public database (as BOLD). With time these BINs will receive a name if the institutions where the collections are deposited keep the policy to continue the taxonomical studies and encourage the specimens' curation.
This previous step will allow the implementation and continuity of non-invasive biomonitoring with all new tools for eDNA currently available. However, this second step will depend on one side of society committed to conservation and sustainable development. On the other side, it will depend on the institutions with the equipment and knowledge to perform these complex but efficient techniques.
The economy of Quintana Roo state in the Yucatan Peninsula depends almost only on tourism. All the visitors are looking for the natural attractions that this region offers [64]. The recent presence of sargassum in the shoreline of the Mexican Caribbean caused that many visitors moved to inland freshwater sinkholes and archaeological sites. As a result, Bacalar (located south of Sian Ka'an) doubled the number of visitors from 2017 to 2018. Currently, it is noticeable an effect of these activities on the microbialites from this site that are unique in the world [65].
Even though Sian Ka'an reserve is a natural protected area, the water systems that support all the equilibrium in this region are pretty vulnerable to the indirect or direct effect of pollution of the underground or surface water, the introduction of exotic species, and all effects derived from mismanaged tourism. We consider it urgent to consider a proposal to monitor, in a fast and reliable way, all future changes in this environment, as suggested here. Funding: This study was financed by the Global Environment Fund through the United Nations Development Programme (UNDP, Mexico), Comisión Nacional para el Conocimiento y Uso de la Biodiversidad (CONABIO) and Comisión Nacional de Áreas Naturales Protegidas (CONANP) as part of the investigation called: Programa de detección temprana piloto de especies acuáticas invasoras a través de los métodos de código de barras de la vida y análisis de ADN ambiental en la Reserva de la Biosfera Sian Ka'an within Project 00089333 "Aumentar las capacidades de México para manejar especies exóticas invasoras a través de la implementación de la Estrategia Nacional de Especies Invasoras". Acknowledgments: Special thanks to José Angel Cohuo Colli, Adrian Emmanuel Uh Navarrete, Ivan Canul Palma, and Georgina Alexandra Prisco Pastrana for the field assistance. Alma Estrella García Morales from the Mexican Barcode of Life (MEXBOL) node Chetumal assisted with DNA extraction, PCR reactions, and sequence edition of all material presented here. We appreciate to Felipe Ángel Omar Ortiz Moreno, Director of the Sian Ka'an reserve, his support and facilities to realize this study. Jorge Vidal Quijada, Heliseo Torres Otero, Gilberto Paredes Pat, Isaías Tun Santiago, Alejandro Tuyub Tuyub from Limones town gave us guidance and advice to move inside and the surroundings of the Sian Ka'an Reserve. Edilberto Sosa Martín from the reserve´s main station supported as well our sampling effort. All ejidal commissars supported our study in their communities. Sarah Dolynskyj and Natalia Ivanova from the Biodiversity Institute of Ontario performed the laboratory processes for metabarcoding. Finally, two anonymous referees gave us valuable advice to improve this manuscript.

Conflicts of Interest:
The authors declare no conflicts of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.