DNA barcodes and evidence of cryptic diversity of anthropophagous mosquitoes in Quintana Roo, Mexico

Abstract Culicidae mosquitoes are potential vectors of pathogens that affect human health. The correct species identification, as well as the discovery and description of cryptic species, is important in public health for the control and management of specific vectors. In the present study, the diversity of anthropophagous mosquitoes in Quintana Roo, at the border between Mexico and Belize, was evaluated using morphological and molecular data (COI‐DNA Barcoding). A total of 1,413 adult female specimens were collected, belonging to eight genera and 31 morphospecies. Most species formed well‐supported clades. Intraspecific Kimura 2 parameters (K2P) distance average was 0.75%, and a maximum distance of 4.40% was observed for Anopheles crucianss.l. ABGD method identified 28 entities, while 32 entities were identified with the BIN system. In Culex interrogator and Culex nigripalpus a low interspecific genetic distance of 0.1% was observed. One undescribed species belonging to the genus Aedes (Aedesn. sp.) was discovered, but no clear genetic divergence was found between this species and the closely related species Aedes angustivittatus. An intraspecific K2P distance greater than 2.7% was observed in Aedes serratus(3.9%), Anopheles crucianss.l. (4.4%), Culex taeniopus (3.7%), Haemagogus equinus (3.9%), Culex erraticus (5.0%), Psorophora ferox (4.5%), and in Anopheles apicimacula(8.10%); therefore, evidences of cryptic diversity are shown in these species. This study showed that DNA barcodes offer a reliable framework for mosquito species identification in Quintana Roo, except for some closely related species for which it is recommended to use additional nuclear genetic markers such as ITS2, in order to resolve these small discrepancies.

. The combined use of DNA and morphology allowed to analyze a greater number of species worldwide and to discover cryptic speciation in taxa of medical importance such as culicides (Ashfaq et al., 2014;Hebert & Gregory, 2005;Laurito, Oliveira, Almirón, & Sallum, 2013;Torres-Gutiérrez et al., 2016;Wang et al., 2012). Though, given the diversity of mosquitoes, it is possible to affirm that integrative taxonomy in this group is at its early stages (Beebe, 2018).
For several years, the Cytochrome Oxidase I subunit (COI) gene as a molecular marker for mosquitoes, has evidenced the presence of species complexes and cryptic species within the genera Aedes, Anopheles and Culex (Ashfaq et al., 2014;Wang et al., 2012). These discoveries have been interesting from an ecological point of view, as molecular analyses demonstrated a wide distribution of several genera, and predicted that speciation processes may be occurring at an unobservable rate due to the constant exposure to physical and biological environmental factors. Therefore, molecular information in addition to predicting the presence of a great number of undescribed species also allows analyzing the variation at the genetic level (Weeraratne et al., 2018).
The discovery of cryptic species in Culicidae also entails serious public health implications as it directly impacts the design of vector control and management programs (Bickford et al., 2007). The correct identification of the species is essential for the development of programs for the control and prevention of diseases transmitted by mosquitoes, as it allows focusing only on the control of those species that transmit certain diseases (Azari-Hamidian et al., 2010;Bueno-Marí, Corella-López, & Jiménez-Peydró, 2010;Erlank, Koekemoer, & Coetzee, 2018;Laboudi et al., 2011).
Unlike other groups of insects, Culicidae mosquitoes are widely studied due to their role in human health (Cardoso et al., 2010;Rozo-López & Mengual, 2015). Nevertheless, their knowledge from a taxonomic point of view is still uncomplete (Kumar, Rajavel, Natarajan, & Jambulingam, 2007). The difficulties in the correct identification of mosquito species derive from the impossibility to use some diagnostic characters, such as scales and setae, which get often damaged during field collection or storage of specimens, and the presence of other characters which are evident only during some stages of their growth (Kumar et al., 2007). Furthermore, isomorphic species can be often found in mosquitos, grouped in closely related species (sibling species) (Beebe, 2018). These issues cause delays in the correct identification of the species (Ruiz-López et al., 2012), besides the fact that a high taxonomic expertise in needed.
There are 3,554 species of mosquitoes recognized worldwide (Harbach, 2018), of these only 1,254 have been barcoded that is <40% (www.Boldsystems.org) (Ratnasingham & Hebert, 2007). In Mexico, around 250 species belonging to 20 genera are reported (Bond et al., 2014;Ibáñez-Bernal, Strickman, & Martínez-Campos, 1996;Ortega-Morales et al., 2015). Among the important genera reported for Quintana Roo is the genus Aedes sensu Wilkerson et al. (2015) potential vector of pathogens causing diseases such as yellow fever, dengue, zika, among other diseases that represent a threat to human and animal health (Ortega Morales et al., 2010;Salomón-Grajales et al., 2012). In the present study, a total of nine species of Aedes genus were identified, although the number of species belonging to this genus is present in the State is still unknown if exist more due to the lack of exhaustive studies involving morphology and DNA barcoding, for this and other species of medical importance.
The objectives of this study were to delimit anthropophagous mosquito species using DNA barcodes and traditional taxonomy in the southeastern region of Mexico, an area designated as malarial on the border between southern Quintana Roo and Belize; to find evidences of cryptic species in designated species with wide distribution and, finally, to contribute to GenBank and Bold systems databases with new DNA barcodes.

| Mosquitoes collection
Adult females were collected in three selected locations of

| Depository and morphological identification
Vouchers were deposited in the Entomology Collection of the Zoology Museum of ECOSUR-Chetumal Unit where they were assigned a reference code (ECO-CH-AR/DP_0285-1697), and lateral view photographs of each voucher were taken. All specimens were identified using different available literature such as books and identification keys, which included Díaz Nájera (1965), Berlin (1969), Arnell (1976), Darsie and Ward (1981), Sirivanakarn (1982), Clark-Gil and Darsie (1983), and Wilkerson and Strickman (1990). The collected data were deposited in the BOLD Systems ® database (www.boldsystems.org) (Ratnasingham & Hebert, 2007) under the project "Mosquitoes from Rio Hondo, Mexican Caribbean" (MRRO). The Culicidae classification system proposed by Wilkerson et al. (2015) and available in the Walter Reed Biosystematic Unit (www.wrbu.org) was used for this study.

| DNA isolation, PCR amplification, and sequencing
For the molecular analysis, a total of 272 specimens that represented 31 morphospecies, were selected. Five to thirty specimens per morphospecies were selected. All specimens were processed by removing one hind leg from each individual, which was deposited in a sterile tube of 0.2 ml with 30 µl of ethanol (Ashfaq et al., 2014).

| Sequences analysis
All sequences obtained in this study were compared with sequences available in the BOLD Systems ® using "Identification Engine" and with sequences available in GenBank with Basic Local Alignment Search Tool (BLAST) (blast.ncbi.nlm.nih.gov/Blast.cgi). MEGA v 6.06 program (Tamura, Stecher, Peterson, Filipski, & Kumar, 2013) was used for sequences alignment, for searching for stop codons and for the calculation of intra-and interspecific genetic distances using the Kimura two parameters distance model (K2P) (Kimura, 1980;Tamura et al., 2013). To provide a graphic representation of the grouping pattern between the different species a K2P distance similarity tree was constructed using the Neighbor-Joining model (NJ) with

| Sequences adding from databases
Seventy-one high-quality public sequences belonging to 13 species were selected from the BLAST and Identification Engine analysis and added for the construction of the tree; specimens were collected in other geographic sites and identified as some of the species reported in Supporting Information Appendix S1. Only sequences belonging to the BOLD Systems database were added to this study, as this database provides information such as collection data, specimen depository, and identifier. For Hg. equinus were adding five sequences generated by us from specimens collected in other locality, theses sequences are still unpublished but have high quality.

| Delimitating species
To delimit the Molecular Operational Taxonomic Units (MOTU's), the ABGD (Automatic Barcode Gap Discovery) program was used (Puillandre, Lambert, Brouillet, & Achaz, 2012) with the following parameters: a minimum intraspecific distance (P min ) of 0.001, a maximum intraspecific distance (P max ) that oscillated from 0.02 to 0.1, the barcode gap width parameter with the default configuration (1.5), and the K2P and Jukes-Cantor (JC) (Jukes & Cantor, 1969)  obtained for this study and those downloaded from BOLD systems.
The BIN system method was also used to delimit the MOTU's in our dataset (Ratnasingham & Hebert, 2013). Psorophora champerico sequence was excluded from the ABGD analysis as this is sensitive to singletons.  (Table 1), within these a new species was discovered Aedes n. sp. (Figure 2).

| Morphological identification
The genus Aedes with nine species was the better represented,

| DNA barcodes
From the 272 specimens processed, 243 COI sequences were obtained. 29 sequences were discarded due to their low quality but the 31 morphologically identified species are represented by their corresponding DNA barcode. The size of the sequences obtained was between 570 and 650 bp. The total number of specimens morphological and molecularly analyzed is reported in The sequence generated in this study for the species Ps. champerico is a new contribution to both databases.

| Genetic distance
In this study, the average K2P intraspecific distance was 0.75%.
The maximum observed K2P distance between conspecific specimens was for An. crucians s.l. with a value of 4.40% (Tables 1 and   3

| Neighbor-joining tree
Most species formed well-supported groups in the NJ tree with bootstrap values between 99% and 100%. (Figure 3).  (Figure 3).

Monophyletic clades were observed for Cx. interrogator and
Cx. nigripalpus, and for Ae. angustivittatus and Ae. n. sp., with high bootstrap values (99%-100%) (Figure 3). Ae. serratus formed welldefined paraphyletic clades (Figure 3) including both the sequences generated in this study (clade QROO) and the sequences from the database belonging to specimens collected by Talaga et al. (2017) in French Guiana (clade FG) (Supporting Information Appendix S1). TA B L E 3 K2P sequence divergence of COI barcode region among the mosquito species with >2 specimens, analysis in Culicidae family among the four genera with two or more species by a single sequence MRRO218-16. In Cx. erraticus, we also found two well-supported clades, QROO clade with six sequences generated in this study and US clade that include the sequences obtained from the database of specimens collected from Florida (Supporting Information Appendix S1) (Figure 3). For Ps. ferox two clades were observed, one consisting of sequences generated in this study mix with sequences downloaded from database and collected in French

Hameagogus equinus
Guiana, Argentina, and Mexico (Supporting Information Appendix S1) (Figure 3), and a second clade formed by seven sequences belonging to US specimens (Supporting Information Appendix S1); both clades with bootstrap values >60%. Finally, for An. apicimacula, were formed two clades, one formed by the sequences generated here and the sequences GBANO930-14 and GBANO931-14 from Colombian specimens, and another clade formed by other sequences from Colombian specimens (Supporting Information Appendix S1) ( Figure 3).

| Delimitation of entities
The ABGD method identified 27 MOTU's for the 31 species (

| Distribution of anthropofagous species in the border Mexico-Belize
In the locality of Sacxan, six species of Culicidae had already been reported (Ortega-Morales et al., 2015), and 17 more species are reported in this study (  (Vargas & Martínez-Palacios, 1950), and its collection in this study suggests that their distribution in Quintana Roo is wide, this species together with An. vestitipennis and An. albimanus, are considered the main malaria vector in Mexico (Hiwat & Bretas, 2011;Loyola, Arrendondo, Rodríguez, Browns, & Vasa Marin, 1991;Mijares, Pérez Pacheco, Tomás Martínez, Cantón, & Ambrosio, 1999). By another hand, Aedes aegypti is the main vector of Dengue, Chikungunya, and Zika virus in Quintana Roo (Sánchez-Rodríguez et al., 2014;Torres-Avendaño et al., 2015). Hence, the importance to continue with the entomological surveillance of these species in the region.

| Genetic distance
In the present study, the overlap observed between the intra-and  Although has been pointed that the taxonomic performance of RESL in BIN system is stronger that ABGD algorithm, showing better results between species identified and MOTU assigned (Ratnasingham & Hebert, 2013). In ABGD, the results improve when is selected a best partitioning scheme, in the analysis the final partitioning depends to a large extent on the user-defined prior upper limit to intraspecific distance (P) as pointed out by Puillandre et al. (2012). Also, it is important to mention that both analyses largely depend on the representation of the species (Pentinsaari, Vos, & Mutanen, 2017).
The results obtained with the ABGD method in this study are similar to those obtained by Versteirt et al. (2015) for the culicids of Belgium, where Cx. pipiens and Cx. torrentium, Ae. cantans and Ae.
annulipes, Ae. punctor and Ae. communis, and An. maculipennis s.s. and An. messeae were grouped in the same MOTU. In Aedes genera, the intraspecific divergence was between 3.6% and 3.9%, evidencing that the delimitation of species in this genera depends largely on the sampling (Beebe, 2018).

| Anopheles apicimacula
Our results supported the hypothesis of Gómez, Bickersmith, González, Conn, and Correa (2015), suggesting that An. apicimacula is a species complex. In this study, when analyzing the sequences generated here and those of specimens from Colombia, the average intraspecific distance (K2P) was 5.6%. Our sequences formed a clade with sequences of An. apicimacula from Antioquia Colombia, while sequences of An. apicimacula from Choco and Valle del Cauca Colombia formed another clade. Due to the proximity between our study area (<300 km) and the type locality (Livingston, Guatemala), it is likely that our collected specimens corresponded to An. apicimacula s.s. (Dyar & Knab, 1906)

| Anopheles crucians s.l.
In Mexico, two species of this complex have been reported using only morphology, An. crucians s.l. which is reported in North America (Massachusetts to New Mexico, USA), Central America (Mexico to Nicaragua) and in the Caribbean islands (Wilkerson, Reinert, & Li, 2004), and An. bradleyi which is distributed along the coasts of the Atlantic and Gulf of the USA up to the south of Nicaragua (Floore, Harrison, & Eldridge, 1976). It is necessary to continue studying An.
crucians species and contribute with information in order to facilitate the separation of all the species of the complex. The genetic distances of 4.4% reported among the specimens collected in this study, suggested the presence of cryptic speciation. This result was not surprising as An. crucians has been previously considered as a complex formed by seven species which are impossible to separate using only adult female characters, and therefore, the use of molecular markers could represent the only way to allow their separation (Wilkerson et al., 2004). However, with the results obtained in this study, we might probably report an eighth species in this complex.

| Aedes serratus
This species showed a greater intraspecific genetic distance than the one previously reported (3.50%); furthermore, the distance between specimens collected in this study (Ae. serratus QROO) and sequences of individuals collected in French Guiana (Ae. serratus FG) showed the existence of two lineages. So far, a wide distribution has been reported for this species from Mexico to Brazil (WRBU, 2018) but it is very likely that this wide distribution will be interrupted by the Andes mountain range that crosses Colombia and Venezuela which is responsible for propitiating the speciation (De-Silva et al., 2017).

| Hameagogus equinus
The three clades observed with an average intraspecific distance of 3.05% also suggest cryptic speciation within the genus. No morphological variation has yet been found in the analyzed specimens but we suggest including the observation of immature stages such as eggs, larvae, and pupae, as well as the genitalia of the adult males to find the trait(s) that separate these species (Versteirt et al., 2015). Recently, Mello, Santos-Mallet, Tátila-Ferreira, and Alencar (2018) found significant differences in the external morphology of the eggs, the exochorion ornamentation, in three populations (the USA, Trinidad and Tobago, and Brazil). These differences were pointed out by Linley and Chadee (1990) in the population of Southeast USA and Trinidad and Tobago; however, they were not considered by the separation of species. Nevertheless, we now show molecular evidences that suggest more than one lineage for Ps. ferox. These results imply the necessity to continue with studies including integrative taxonomy of Ps. ferox populations throughout the American continent, which would help delimiting new species and establishing their geographic range (Dayrat, 2005;Will, Mishler, & Wheeler, 2005). Mendenhall, Bahl, Blum, and Wesson (2012)

| Culex taeniopus
The usefulness of DNA barcodes for the identification of species of the subgenus Melanoconion has been reported, as well as the discovery of cryptic species or new species within the genus Culex (Laurito et al., 2013;Torres-Gutiérrez et al., 2016). Values of 3.70% are reported for specimens of Cx. taeniopus collected for this study. The specimens were collected in relatively close localities (Palmar and Ramonal), with a distance between localities of 12 km, and no apparent geographic barrier that can lead to speciation was present; therefore, the cryptic speciation could be the result of a sympatric speciation. To clarify the high values of intraspecific divergence, we suggest extend the collection of specimens to increase the number of sequences and represent a greater geographical area, at this moment Guatemala is only represented with one specimen, performing ecological studies and observation of structures such as genitalia (in males) and the morphology of larvae and eggs to find any character that could separate the species. It is in our interest to rule out the possibility of the presence of a complex of species and not only cryptic speciation.

| Morphological differences versus low genetic distances
The interspecific distance between Ae. angustivittatus and Ae. n. sp. of 1.1%, the low average is due that probably both species belong to the Scapularis group. However, morphological characters of Ae. n. sp.
did not match those of Ae. angustivittatus, nor those of Ae. trivittatus or the characters of the hybrid reported by Arnell (1976) (Laurito et al., 2013). The same case was reported between Culex minor and Culex spiculosus, which presented an average K2P interspecific distance of 1.86% (Wang et al., 2012). These results could be explained by an incomplete lineage sorting or introgression events (Beebe, 2018).
Cases like those here reported, where COI does not present enough divergence to separate close species because the lineages are not fully sorted into divergent clades, are not new, but our data increase the lineages record of recent divergence. In this sense, it is important to conduct further studies in order to look for evidence of reproductive isolation. Also, the use of a single gene region is an imperfect tool because will be overlooked recent species because of their low sequence divergence (Ratnasingham & Hebert, 2013). It is necessary to include at the analysis nuclear marker like ITS2 as well as ecological and biological data (Ajamma et al., 2016;Alquezar, Hemmerter, Cooper, & Beebe, 2010;Beebe, 2018;Mardulyn, Othmezouri, Mikhailov, & Pasteels, 2011;Versteirt et al., 2015;Wang et al., 2012).

| CON CLUS IONS
The results in this study contribute to the development of a reference DNA barcode library for Mexican culicids. In addition, our analysis reported the presence of a certain taxonomic differentiation that needs to be investigated in An. apicimacula, An. crucians s.l., Ae. serratus, Hg. equinus, Ps. ferox, Cx. erraticus, and Cx. taeniopus. For the new specie discovered, it is necessary to work on the description during all their life stages (larval exuvia, pupa, and adults), which allow the documentation of the link between morphological and molecular identification standards (Versteirt et al., 2015). Finally, the identification of species is an essential step for monitoring and vectors control. This information is valuable for the Ministry of Health provided taxonomic support on identification. All authors contributed to the text and approved the final manuscript.

DATA ACCE SS I B I LIT Y
The collected data are deposited in the database BOLD Systems ® (www.boldsystems.org) in the project "Mosquitoes from Rio Hondo, Mexican Caribbean" (MRRO). Molecular sequences have been submitted to GenBank on June 13.