DNA Barcoding of Mosquitoes from the Pantanos de Centla Biosphere Reserve, Southeastern Mexico

The family Culicidae is one of the most important arthropod groups of medical importance due to the large number of pathogens that some species can transmit to animals and humans. There are more than 3,500 species worldwide (Harbach 2013), while the number of recorded species in Mexico varies between 240 species (Rodríguez-Martínez et al. 2020) and 250 species (Bond et al. 2014, Chan-Chable et al. 2019). Historically, studies of mosquito taxonomy in Mexico have mainly focused on the tropical forests of the southeastern region, which includes the richest biodiversity of the country. Recently, studies on the taxonomy and distribution of mosquitoes have been undertaken in several states, but the taxonomic studies on the mosquito fauna that inhabit protected natural areas and biomes are limited worldwide (Dutta et al. 2010, Santos et al. 2015).

In Mexico, there are 44 protected natural areas declared as Biosphere Reserves, which makes Mexico the country with the largest number of biosphere reserves in Latin America (SEMARNAT 2018). Pantanos de Centla, a large wetland situated in Tabasco State, was declared as a Biosphere Reserve (PCBR) on August 6, 1992, and, at present, is the largest wetland in Mexico and one of the largest in North America (SEMARNAT 2016). Several vector-borne diseases are recorded in surrounding areas of PCBR since environmental conditions are favorable for the development and maintenance of large mosquito populations (Ulloa et al. 2003). For example, dengue fever and malaria are frequently reported in patients living within the PCBR, which is now considered an endemic area for these diseases. In addition, outbreaks of chikungunya virus and Zika virus (ZIKV) have also been reported in patients living in the area (SINAVE 2019). Dengue fever virus is the most important mosquito-borne pathogen in Tabasco State, with 1,630 cases during 2015–19, of which 1,135 were febrile patients and 495 patients developed hemorrhagic fever (SINAVE 2019). Chikungunya fever is an emerging disease, with 59 cases reported in Tabasco between 2015 and 2016, although there have not been any new cases reported since late 2016 (SINAVE 2019). Zika is another emergent disease, with 381 cases reported in Tabasco between 2015 and 2019, of which 286 were reported in pregnant women. Although no cases of microcephaly in fetuses or newborns have been reported, 3 cases of Guillain-Barré syndrome associated with ZIKV were reported in the state (SINAVE 2019). West Nile virus and Venezuelan equine encephalitis virus have also been reported in the state (Obreste et al. 1998, Estrada-Franco et al. 2003, Hidalgo-Martínez et al. 2008, Adams et al. 2012). Moreover, the PCBR is a tourist area and outdoor activities are frequently suspended due to the large number of biting mosquitoes, which negatively affects the economy of the region.

Adequate taxonomic identification of mosquito species is essential for the establishment of surveillance and control programs (Ruiz-Arrondo et al. 2019). However, morphological identification is difficult because of similar morphology between life stages in several species and the presence of species complexes (Cook et al. 2005, Chan et al. 2014, Batovska et al. 2016, Hernández-Triana et al. 2019). To overcome this taxonomic impediment, a small region (658 bp) of the mitochondrial cytochrome c oxidase subunit I (COI) gene has been proposed for species identification (DNA barcoding) (Hebert et al. 2003a, 2003b; Laurito et al. 2013).

The DNA barcoding methodology has been used extensively to identify mosquito species in different geographical regions worldwide, including Australia (Batovska et al. 2016), Belgium (Versteirt et al. 2015), Canada (Cywinska et al. 2006), China (Wang et al. 2012), Colombia (Hoyos-Lopez et al. 2015), Ecuador (Linton et al. 2013), French Guiana (Talaga et al. 2017), India (Kumar et al. 2007), Mexico (Adeniran et al. 2021), Pakistan (Ashfaq et al. 2014), Singapore (Chan et al. 2014), Spain (Ruiz-Arrondo et al. 2019), Sri Lanka (Weeraratne et al. 2017), Sweden (Engdahl et al. 2014), Turkey (Gunay et al. 2015), and the United Kingdom (Hernández-Triana et al. 2019). However, there is insufficient information on the use of this technique within the Culicidae in Mexico. The DNA barcoding has been used to confirm the presence of Psorophora albipes (Theobald) and Anopheles veruslanei Vargas in Quintana Roo State (Chan-Chable et al. 2016, 2018a), to support the presence of cryptic diversity in Aedes taeniorhynchus (Wiedemann), as well as to support the morphological identification and the presence of cryptic diversity within 7 species in Quintana Roo (Chan-Chable et al. 2018b, 2019). In the present study, the DNA barcoding approach was used to support the identification of the local mosquito fauna in the PCBR. In addition, the DNA barcode variability was assessed using genetic distance methods to detect cryptic diversity across different mosquito species.

The PCBR is located in Tabasco State, Mexico (18°20′00″N, 92°30′00″W). It has an area of 3,027 km², and borders the Gulf of Mexico in the north, the Bitzal and Grijalva rivers in the south, Campeche State in the east, and the Las Porfirias river and city of Villahermosa in the west. The reserve belongs to the physiographic region of the Coastal Plain of the Southern Gulf, subregion of Tabasco, Plains and Marshes, and includes the municipalities of Centla, Jonuta, and Macuspana (Fig. 1). The climate is warm-humid and tropical subhumid, with rains during the summer with an annual average rainfall of 1,400–1,800 mm (INEGI 2018).

Adult mosquitoes were collected from locations inside and within close proximity to the PCBR (Table 1) during July and November of 2016, which included the dry and rainy seasons. Collections were carried out using 10 Centers for Disease Control and Prevention (CDC) light traps located between 1 and 10.5 m from ground level, baited with octenol, and operated at night (1800–2200 h); 2 Shannon traps with humans employed as baits were used at night (2000–0300 h); and mosquitoes were also collected from resting places using 2 Insectzookas (BioQuip No. 2888A, Rancho Domínguez, CA) during the day between 0900 and 1700 h. Immature stages were collected directly from aquatic habitats. Larvae and pupae were stored in individual tubes to obtain the adult stages and associated exuviae. Adults were killed using lethal chambers with triethylamine vapors, stored in vials, and preserved in liquid nitrogen during transportation to the laboratory. All specimens were transported to the Molecular Biology Laboratory of the Parasitology Department of the Universidad Autónoma Agraria Antonio Narro, Unidad Laguna, for mounting and morphological identification. Adult mosquitoes were mounted on insect pins and identified using taxonomic keys (Sirivanakarn 1982, Clark-Gil and Darsie 1983). The classification system proposed by Wilkerson et al. (2015) was used for the tribe Aedini, while for the other tribes and Anophelinae the classification system of Knight and Stone (1977) was used.

A modified Hotshot technique (Montero-Pau et al. 2008, Hernández-Triana et al. 2019) was used for DNA extraction. In brief, 1 or 2 legs from individual specimens were put directly into 50 μl of alkaline lysis buffer in a 96-well plate, and sonicated in a water bath for 20 min. Subsequently, the plate was incubated in a polymerase chain reaction (PCR) block machine for 30 min at 94°C, and allowed to cool for 5 min at 4°C, after which 50 μl of the neutralizing buffer was added. The plate was stored at −80°C until processing the following day.

The PCR amplification was carried out using the primers LCO1490 and HCO2198 (Folmer et al. 1994), which are considered the standard for the amplification of the 658-bp barcode region located at the 5′ end of the COI gene (Hebert et al. 2003a, 2003b). The PCR products were obtained using the protocol of Hernández-Triana et al. (2017). A 1.5% agarose gel was used to visualize the PCR products, and samples showing the correct band size were sequenced in both directions using the ABI PRISM® BigDye® Terminator v3.1 Cycle Sequencing Kit (Applied Biosystems, Waltham, MA) at the Sequencing Unit, Animal and Plant Health Agency (Surrey, United Kingdom).

All bidirectional sequences were combined to produce a single consensus sequence, the full 658-bp barcode. The full data set was analyzed in MEGA v.6 (Tamura et al. 2013), and genetic relationships between species were analyzed using the neighbor joining (NJ) and maximum likelihood (ML) default values; Bayesian inference (BI) analysis was conducted in MrBayes 3.2 (Ronquist et al. 2012). The NJ and ML were undertaken using the Kimura 2-parameter (K2P) distance metric to represent their distribution pattern in the trees. For species that were not found in the PCBR but have been recorded for the region, we used COI barcode sequences for those species from Quintana Roo (Chan-Chable et al. 2019), as well as sequences available in Barcode of Life Data System (BOLD) (see Table 1). The robustness of the NJ and ML trees was calculated using the bootstrap methodology employing 1,000 pseudoreplicates; only groups with more than 80% support were mapped in the NJ tree (Tamura et al. 2013, Hernández-Triana et al. 2019). The BI analysis was conducted with 2,000,000 generations running and sampling sequences of 100 generations (Ronquist et al. 2012). Barcodes longer than 500 bp were allocated a Barcode Index Number (BIN) (Ratnasingham and Hebert 2013). Each BIN was then mapped onto the NJ tree to examine its distribution among every morphologically identified species.

A total of 1,402 females, 24 males, and 150 larvae were collected (1,576 specimens in total). Based on morphology, specimens were identified as belonging to the 2 subfamilies present in Mexico (Anophelinae and Culicinae), 6 tribes (Aedeomyiini, Aedini, Culicini, Mansoniini, Sabethini, and Uranotaeniini), 10 genera (Anopheles Meigen, Aedeomyia Theobald, Aedes Meigen, Psorophora Robineau-Desvoidy, Culex Linnaeus, Deinocerites Theobald, Coquillettidia Dyar, Mansonia Blanchard, Wyeomyia Theobald, and Uranotaenia Lynch Arribálzaga), and 35 species (Table 1). The genus Culex with 11 species was the most represented, while Aedeomyia, Deinocerites, and Wyeomyia were represented by only a single species each (Table 1). The most abundant species were An. vestitipennis Dyar and Knab (n = 800), An. albimanus Wiedemann (n = 121), and Ma. dyari Belkin, Heinemann, and Page (n = 89). Within the 35 species collected, 15 species were identified that are of medical and veterinary importance (Table 2).

In total, 225 COI DNA barcode sequences for 30 species were analyzed, which represents 85.7% of the species recorded at the PCBR (Table 1). The sequences generated in this study for De. pseudes Dyar and Knab and Ma. dyari are new additions to BOLD and National Center of Biotechnology Information (NCBI). The genetic diversity within our data set was analyzed using the NJ, ML, and BI methods. In general, the trees obtained showed a similar specimens topology with strong bootstrap support values; therefore, we illustrated only the NJ tree (Fig. 2). Only 2 sequences identified as Ae. scapularis (Rondoni) (CBMA045-12 and CBMA016-12) separated from the remainder of the specimens identified as this species in the ML analysis (tree not shown). Similarly, the BI analysis showed that 5 sequences generated from specimens identified as Ae. scapularis (94CBMA016-12, 102CBMA105-12, 152MF172267, 153MF172266, and 154MF172265) grouped separate from other specimens.

As a whole, most species formed well-supported groups with bootstrap values higher than 92% (Fig. 2); however, we found 2 well-supported groups within An. crucians s.l., one from Mexico (BOLD:AAC8253), supported by 2 sequences generated in this study; and one group from Florida, USA (BOLD:AAA5102), supported by sequences obtained from BOLD (Table 1 and Fig. 2). Two lineages were observed within Cx. erraticus (Dyar and Knab): one from PCBR, Mexico (BOLD:AAG3848), supported by 1 sequence generated in this study; and one group from Florida, USA (BOLD:ADR1028), supported by 5 sequences and obtained from BOLD (Table 1 and Fig. 2). Finally, 2 groups were obtained in Cx. conspirator, one group supported by 2 sequences (BOLD:ACU5030 and BOLD:ACU5302) from Colombia. Regrettably, since no males were collected, it was not possible to reliably distinguish the females of Wy. abebela Dyar and Knab and Wy. melanopus Dyar based on morphological traits; thus, we identified these specimens as Wy. abebela/melanopus.

The levels of sequence divergence were variable across the taxa, with conspecific individuals collected from a single site often exhibiting 0 (De. pseudes) or between 0.06–1.60% divergence values, while other specimens showed higher percentages (Table 1). In this study, the average K2P intraspecific distance was 1.52%. The maximum observed average K2P distance above 2% between conspecific specimens was 6.09% for female of An. crucians s.l., while Cx. erraticus and Cx. conspirator showed an average intraspecific distance of 1.5% (maximum average 2.50%) and 1.67% (maximum average 3.51%), respectively (Table 1). The interspecific divergence ranged between 2.29% and 24.3% (Table 1). As estimated, the smallest values of genetic divergence were identified among species in the same genus: for example, Aedes scapularis/Ae. euplocamus (Dyar and Knab) (2.69%), and Psorophora albipes/Ps. confinnis (Lynch Arribálzaga) (8.8%).

The 225 barcodes analyzed in this study produced 33 BINs, which were representative of 30 morphologically identified species. Of the 33 BINs, 25 were concordant with existing BINs and 3 were singletons (Cx. conspirator [BOLD:ACU5302], Cx. erraticus [BOLD:AAG3848], and De. pseudes [BOLD:ADE5088]). Furthermore, 2 species showed 2 BINs each: An. crucians (BOLD:AAA5102 and BOLD:AAC8253) and Cx. erraticus (BOLD:ADR1028 and BOLD:AAG3848). In addition, 2 species, Ae. euplocamus and Ae. scapularis, shared the same BIN (BOLD:AAH9007) (Fig. 2).

These data provide a faunistic survey of 35 mosquito species collected in the PCBR. Although 35 species records are documented in this data set, there are likely to be more species present in this region that were not collected in this study. Uncollected species might be due to the collection methods that were used, which were mainly aimed at collecting species attracted by CDC light traps and Shannon traps only, and in a few cases, immature stages from larval habitats. In addition, access to PCBR is challenging, and as a result, it was not possible to obtain specimens from all localities. All 35 species collected in this study were previously recorded in the survey undertaken by Ortega-Morales et al. (2019a) in Tabasco State.

Based on the ecology of mosquito species identified in this study, 4 ecological categories may also be recognized. The 1st category includes species that develop in swamps and temporal ponds at ground level with emergent vegetation. Taxa in this category include the majority of species recorded in the PCBR (n = 26, 74.2%), especially all members of the genera Aedeomyia, Psorophora, Coquillettidia, Mansonia, Uranotaenia, and some species of Anopheles, Aedes, and Culex. A 2nd category includes generalist species that develop in both swamps and temporal ponds at ground level, and inside natural and/or artificial containers. This category includes all species of the subgenus Culex (Culex), and Cx. (Phenacomyia) corniger Theobald (n = 4, 11.4%). The 3rd category includes species that develop mainly in either natural and/or artificial containers. This group includes species of Aedes (Stegomyia) and Wy. abebela/melanopus (n = 3, 8.5%). Both Ae. aegypti (L.) and Ae. albopictus (Skuse) are commonly discovered in artificial containers inside PCBR, while Wy. abebela/melanopus has been collected only in bromeliad axils. Finally, the 4th category comprises species that develop in crab holes, mainly De. pseudes (n = 2, 2.8%). Crab holes are common larval habitats in the PCBR. These habitats can be exploited by several mosquito species, mostly during the dry season when the water level of the swamps decreases and a large number of islets emerge, and the crab holes fill with fresh water from surrounding rivers and streams.

The values for the intraspecific and interspecific genetic divergences obtained in this study are within the values obtained by other studies in the family Culicidae, for example, Hoyos-Lopez et al. (2015), Talaga et al. (2017), Chan-Chable et al. (2019), and Hernández-Triana et al. (2019). The morphologically identified specimens of the same species formed well-defined groups in the NJ, ML, and BI analysis (Fig. 2), supporting the use of DNA barcoding in combination with morphological characters as a suitable approach for species identification. Nonetheless, out of the 30 species that showed a correspondence between morphological and molecular data, 3 species showed deep splits in the NJ tree (Fig. 2): Anopheles crucians s.l., Cx. erraticus, and Cx. conspirator.

The Crucians Complex (Floore et al. 1976, Wilkerson et al. 2004) of Anopheles (Anopheles) currently includes 7 species: 3 nominal species: An. crucians s.s., An. bradleyi King, and An. georgianus King; and 4 unnamed sibling species, which are apparently only identifiable using molecular approaches (Wilkerson et al. 2004). Anopheles crucians s.l. is distributed from the USA to Costa Rica and the Caribbean islands. The only species that have been reported in Mexico are An. crucians s.l. and An. bradleyi (Chan-Chable et al. 2019). Anopheles crucians s.l. has been reported in the northern states of Nuevo León and Tamaulipas (Vargas and Martínez-Palacios 1956, Ortega-Morales et al. 2019b), in the southern states of Veracruz, Yucatán, Quintana Roo, Tabasco, and Chiapas (Vargas and Martínez-Palacios 1956; Ulloa et al. 2009; Ortega-Morales et al. 2010, 2019a; Beltrán-Aguilar et al. 2011; Chan-Chable et al. 2018a), and in the central region of the country, San Luis Potosí and Hidalgo (Vargas and Martínez-Palacios 1956). San Luis Potosí and Hidalgo are located inland; thus, it is likely that the species belonging to the An. crucians Complex from those states are different from those reported in the northern and southern coastal regions. In the PCBR, members of the An. crucians Complex were collected using a Shannon trap placed on an islet during the dry season at night, and in the same location using CDC light traps with octenol as an attractant. The average genetic divergence of 6.09% observed between COI barcode sequences of members of An. crucians s.l. from PCBR and USA are similar to those from Quintana Roo where genetic distances of 4.40% have been reported (Chan-Chable et al. 2019). These high values of genetic distance support other evidence that An. crucians is a species complex, which is also supported by Wilkerson et al. (2004).

The subgenus Melanoconion of Culex includes 25 species in Mexico, 16 of which occur in Tabasco and 7 in the PCBR. This subgenus is divided into 2 species-groups, the Melanoconion and Spissipes Sections; both Sections occur in Mexico (Clark-Gil and Darsie 1983). Adult females in the Erraticus Group of the Melanoconion Section can be identified by the presence of a few to several scales forming a distinct patch on the upper corner of the mesokatepisternum, with the median surface of the mesanepimeron bearing a broad whitish scale-patch (Sirivanakarn 1982, Clark-Gil and Darsie 1983). Melanoconion mosquitoes with a distinctive patch of scales on the median surface of mesanepimeron are usually collected in the tropical and subtropical regions of Mexico. The typical form of Cx. erraticus is presumably the only species of the Melanoconion Section in Mexico with a scale-patch on the mesanepimeron. Culex erraticus is a Neotropical species that occurs from the USA to Middle and South America and some Caribbean islands; in Mexico, this species has been recorded in the states of Campeche, Guerrero, Michoacán, Morelos, Nuevo León, Oaxaca, Quintana Roo, Sinaloa, Sonora, Tabasco, Tamaulipas, Veracruz, and Yucatán (Díaz-Nájera and Vargas 1973). Chan-Chable et al. (2019) stated that COI barcode sequences analysis of Cx. erraticus specimens collected from Quintana Roo and Florida showed a maximum intraspecific distance of 5% (average 4.1%, n = 17), which implied the presence of cryptic diversity in this species. In this analysis, the average genetic distance between Cx. erraticus from the PCBR and Florida was 1.52% (n = 6), and a maximum intraspecific distance of 2.5% was observed; the deep division in the NJ tree (Fig. 2) suggests the presence of 2 lineages. The latter also corroborated the findings in Mendenhall et al. (2012) while analyzing sequences of the ITS2 and NADH dehydrogenase genetic markers in Cx. erraticus.

Culex conspirator has a Neotropical distribution, where it has been reported from Belize, Colombia, Costa Rica, Ecuador, El Salvador, Guatemala, Honduras, Mexico, Nicaragua, Panama, Trinidad and Tobago, and Venezuela (Gaffigan et al. 2017). The sequences analyzed originate from Colombia for Cx. conspirator showed genetic diversity of 1.67%, and a maximum intraspecific distance of 3.51%. As in the previous case, the deep split in the NJ tree might indicate the presence of different lineages within this species. However, the latter hypothesis will require further investigation employing more sequences across the distribution range of both Cx. erraticus and Cx. conspirator (Fig. 2).

Aedes scapularis is morphologically similar to Ae. euplocamus (Chan-Chable et al. 2019); thus, Basic Local Alignment Search Tool (BLAST) searches were undertaken for these sequences in the NCBI database to reconfirm their identification. The BLAST searches retrieved matches to Ae. scapularis with 98.7–99.6% (accession no. MF172267) support for their correct identification.

The present work brings together the available information regarding the utility of DNA barcoding to support the identification of the mosquito fauna in the PCBR. In addition, this study reports the presence of cryptic diversity in An. crucians s.l., the presence of potential different lineages in Cx. erraticus and Cx. conspirator, which will need to be investigated further. This study adds COI sequence data to the BOLD and NCBI databases for De. pseudes and Ma. dyari, and highlights both the medical and veterinary importance of key species that might be involved in pathogen transmission in the area. For poorly represented genera and species, further collecting efforts should be supported to increase the number of specimens and barcode sequences in public databases. This study supports the need for continuing research combining the use of molecular methodologies with morphological characters for mosquito species identification in Mexico.

We are grateful to Salvador Morales-Avitia and Ramón Méndez-López for help during the field collection trips, to the personnel of the biological station “Tres Brazos” for the facilities provided during our collection trips, and to Carlos E. Ramírez-Contreras for providing vehicles to transport equipment to the PCBR. We thank the Secretaría de Educación Pública (PRODEP, Mexico) (Grant No. 13-30-8257-7260 CUAC1414) and Consejo Nacional de Ciencia y Tecnología, Mexico (MEXBOL, Grants No. 251085 and 271108), for funding. L.M. Hernández-Triana thanks the United Kingdom Department for Environment Food and Rural Affairs (DEFRA), Scottish Government, and Welsh Government through grant SV3045, and the European Union Framework Horizon 2020 Innovation Grant, European Virus Archive (EVAg, grant no. 653316) for providing funds for this study. M.A. Rodríguez-Pérez thanks Instituto Politécnico Nacional for providing logistical support during field work. J.A. Garza-Hernández thanks SEP-PRODEP (Reference No. UACJ-PTC-399; S11-6/18-9518) for funding that permitted the completion of the present manuscript. L.M. Hernández-Triana thanks Sean Prosser for reviewing the manuscript. Special thanks to the Sequencing Facility at the Animal and Plant Health Agency, especially Jeremy Hawthorn and Steve Shankster, for technical support. In addition, L.M. Hernández-Triana thanks Clara Berrón, Universidad Nacional del Litoral, for the provision of COI sequences from Argentinean mosquitoes, and Pramual Pairot, Mahasarakham University, Thailand, for his help with the Bayesian Inference analysis.