Anopheles neivai Howard, Dyar & Knab, 1913, is distributed from the Yucatán Peninsula (México) to the Guayas province (Ecuador) (1). In Colombia, it has been recognized as a secondary malaria vector in the Pacific Coast with local importance in small towns such as Charambirá, Santa Bárbara Iscuandé and Buenaventura (2-4). There are also reports of naturally infected females with Plasmodium falciparum in Buenaventura (5). As regards its biology, immature stages develop in bromeliads at mangrove (1). Activity in adults occurs at dawn and dusk (2) when younger females prefer to feed (6) and can be found in homes and in their surroundings (3,7). There are also biting reports in mangrove environments and in boats during fishing activities, which increase the epidemiological risk for local habitants (2).
Anophelines are recognized as a complex group, where it is common to find complexes of cryptic species (8). In Neotropical Anopheles, there are recognized cryptic species complexes that have importance as malaria vectors such as A. triannulatus (9) and A. albitarsis (10). For A. neivai, previous research suggests the existence of a species complex (2,11) based on morphological polymorphism and exploratory information on molecular variability (12).
The use of mitochondrial DNA (mtDNA) has played a vital role in understanding the phylogenetic relationships in Anopheles, including diversification patterns, divergence time estimation (13,14) and the recognition of cryptic species complexes (15-17). There are several advantages to using mtDNA, including variable mutation rates in eukaryotes, lack of introns, and high abundance of copies, which makes mtDNA relatively easy to obtain even from degraded samples (18). In contrast, phylogenetic signals of mtDNA can be altered by mitochondrial introgression (19) and by indirect selection and linkage disequilibrium caused by maternally inherited symbionts in arthropods (20). For phylogeny studies, fragments of cytochrome oxidase I gene (COI), have been extensively used (21-23). However, other genes in the mtDNA genome can be useful for reconstructing evolutionary relationships, including those involved in the respiratory chain complex I and III such as NADH dehydrogenase subunit I (ND1) and cytochrome b (Cytb) (24,25).
In addition, the transfer RNA genes (tRNA), also present in the mitochondrial genome, have a pivotal role in protein biosynthesis (26). The critical function of these genes limits the mutations they can accommodate, which is a desirable tool for studying evolution among closely related organisms. This is particularly important where selective pressure, mutation bias and genetic drift become relevant factors in modeling the codon use for each species (27,28). Furthermore, comparative analysis of the secondary structure of these genes improves the quality of alignments and, thereby, leads to more robust phylogenetic inferences (29,30).
In this work, we explored the phylogenetic utility of a mitochondrial region with fragments of Cytb, ND1 and the secondary structure models of serine transfer RNA (SertRNA) variants to differentiate specimens identified morphologically as A. neivai.
Materials and methods
Specimens from Kerteszia were collected from six localities in Chocó, Colombia: Bahía Solano (Playa Potes and Nabugá), Litoral de San Juan (Charambirá), Acandí, Nuquí and Jurubidá, and the Caribbean region near the Sierra Nevada de Santa Marta. In Central America, the type locality for A. neivai at Portobelo (Panamá) was sampled, as well as two localities in Guatemala: Puerto Barrios and Chiquimula (figure 1, table 1). All the immature stages were collected in bromeliads using manual pipetting and reared under laboratory conditions to obtain adults (31), while field adults were collected using Shannon traps. Specimen identification was performed with dichotomous identification keys (32,33).
For DNA extraction, two legs were removed from each adult, whereas partial abdomens were used for the larval stages. The DNA extraction procedure for all specimens was based on a grind buffer protocol (34). The subsequent PCR conditions for the Cytb-tRNASer-IG1-ND1 region were as follows: 17.1 µL (dH2O); 6 µL (buffer 5X); 2 µL (dntp [2.5 mM]); 1 µL (MgCl2 [25 mM]); 0.3 µL Gotaq® Flexi (Promega); and 1 µL [10 mM] for primers CB3CF (CAYATTCAACCWGAATGATA) and NINFRGGT AYWTTGCCTCGAWTTCGWTATGA) (35). Primer sequences were modified from previously used primers CB-J-11338 and NI-N-11841 (36). The PCR reactions used the following conditions: 94°C/30 s, 47°C/30 s, 72°C/60 s. All the PCR experiments were performed on a PTC-100 BioRad ® thermal cycler. Subsequent fragments were examined by agarose gel electrophoresis (1%). Successive positive PCR fragments were sequenced using an ABI 3500XL ® (Applied Biosystems) automated capillary electrophoresis sequencer at the Centers for Disease Control and Prevention (CDC), Atlanta.
Consensus sequences were assembled using Geneious 6.0 (37). Sequences were checked against available records in the National Center for Biotechnology Information (NCBI) (http://www.ncbi.nlm.nih.gov/) using the BLAST program (38). In addition, the sequences were verified for possible mitochondrial copies at the nucleus (NUMT), following the procedure suggested by Hlaing, et al. (39). The sequences of each region were aligned using Clustal X 2 (40) and included reference sequences of A. albitarsis (NC_020662) from NCBI.
The utility of this region was first assessed by identifying polymorphic sites supporting variants between specimens of A. neivai with other specimens from the subgenus Kerteszia and A. albitarsis (subgenus Nyssorhynchus) using DnaSP 5 (41). Then, tRNASer secondary structures for each variant were modeled using the RNA Fold package available in Vienna RNA 2.0 (42). Parameters were based on Mathews, et al. (43) using minimum free energy and partition, and no GU pairs at the end of helices, allowing dangling energies on both sides of the helix at 37 °C. The phylogenetic signal of this region was estimated using maximum parsimony (MP) with bootstrap resampling using 100,000 replicates through MEGA 7 (44), and a maximum likelihood (ML) with bootstrap support of 1,000 replicates. A Generalized Time Reversible modelGTR (45) was implemented using PhyML (46), and Bayesian inference (BI) under GTR model, assuming relaxed clock rates and Yule speciation process for BI (47) using BEAST 2.4 (48).
Results
Adults from A. neivai, including topotypic specimens from Portobelo (Panamá), were identified based on the following set of characters from (33): mesanepimeron with an upper row seta (figure 2A), abdominal terguites I to IV without dark decumbent scales (figure 2B), acrostical area without white setae (figure 2C), and scale patterns on hind tarsomers (figure 2D). Additionally, specimens from immature stages were identified using the VIth abdominal segment and C setae from VIIth abdominal segment described in (32). No variability for the two immature specimens of A. pholidotus was observed.
![](/img/revistas/bio/v37s2//0120-4157-bio-37-s2-00143-gf2.jpg)
Figure 2 Diagnostic characters used for the identification of A. neivai as described by Harrison, et al. (33). A. Detail of mesanepimeron (Mesu) with upper setae. B. Detail of abdominal segments I to IV (1-IV), without dark decumbent scales. C. Detail of scale patterns in hind tarsomer (Ta-III5). D. Detail of acrostical area (AcS), without white scales
The entire set of specimens produced a total of 36 sequences (table 2): 34 for A. neivai and two for A. pholidotus. Consequently, the amplified fragment of 411 bp consisted of 304 bp for Cytb, 68 bp for tRNASer, 25 bp for IG1, and 11 bp for ND1. Among them, 51 polymorphisms were identified for Cytb, six for tRNASer (including 2 indels), eight for IG1 (including 5 indels) and four for ND1.
Table 2 Specimen collection data for A. neivai and A. pholidotus
![](/img/revistas/bio/v37s2//0120-4157-bio-37-s2-00143-gt2.jpg)
Collection method: Manual pipetting (Mp.), Shannon trap (Sh.) 3 Stage used for specimen identification/DNA sequencing: L (larvae), F (female), M (male) 4 Specimen voucher deposited at Museo Entomológico Francisco Luis Gallego 5 Accession records at the National Center for Biotechnology Information
Three variant secondary structures were associated with tRNASer polymorphisms (table 3). The first two structures (A and B) were found with A. neivai specimens, and the other (C) was found in A. pholidotus. In the A. neivai variants, the only difference is located at the dihydrouridine loop (DHU) as a result of a base substitution (U → C) (figure 3). When comparing A. neivai and A. pholidotus further differences were evident, including a length difference in the DHU loop as a consequence of an adenine insertion, followed by a transversion in A. pholidotus. All the remaining differences between A. neivai and A. pholidotus consisted of transversions in the pseudouridine loop (TΨC) as synonymous substitutions detected in A. pholidotus (A → U, A → U and U → A)
Table 3 Secondary structure (tRNASer) and intergenic spacer (IG1) variants for A. neivai and A. pholidotus
![](/img/revistas/bio/v37s2//0120-4157-bio-37-s2-00143-gt3.jpg)
![](/img/revistas/bio/v37s2//0120-4157-bio-37-s2-00143-gf3.jpg)
Figure 3 tRNASer secondary structure models. A. A. neivai (variant A). B. A. neivai (variant B). C. A. pholidotus (variant C)
The intergenic spacer (IG1) between tRNASer and ND1 exhibited differences in length. A larger intergenic spacer was present in most of the specimens of A. neivai (24-25 bp), in contrast to A. pholidotus (20 bp) and A. albitarsis (18). The nucleotide variability in this region for A. neivai consisted of four variants, base substitutions found in variants 1 to 3, and an insertion on position 373 for variant 4 (Chiquimula, Guatemala). The same insertion was found in A. pholidotus (figure 4).
The comparison of three different approaches (MP, ML and IB) revealed similar patterns. One group with unresolved specimens for the majority of localities and another group including specimens from Panamá and Guatemala (figure 5). The overall pattern of nucleotide substitution estimated by ML across the Cytb-tRNASer-IG1-ND1 regions recovered from several trees revealed phylogenetic signal for the Cytb fragment and tRNASer gene (figure 5b). Equivalent groupings based on trees that included both fragments (figure 6a and 6b), were evident in specimens of A. pholidotus from Colombia and in A. neivai from Central America (Panamá and Guatemala localities). In contrast, when these regions were excluded, leaving only the IG1ND1 fragment (figure 6c), the relationships between A. neivai specimens in the Colombian localities and those from Central America were impossible to resolve due to the presence of polytomies.
![](/img/revistas/bio/v37s2//0120-4157-bio-37-s2-00143-gf5.jpg)
Figure 5 Comparison of tree topologies for A. neivai based on several different analyses for Cytb-tRNASer-IG1-ND1. A. Maximum parsimony. B. Maximum likelihood (GTR model). C. Bayesian inference (GTR model). Bootstrap branch support for each group is provided for maximum likelihood and parsimony. Posterior probability for each node is provided for Bayesian inference
Discussion
Phylogenetic signal is the tendency of closely related species to resemble one another (49). It is also a property of stochastic evolution along a hierarchical tree, depending on sample size, mutational rate and branch length (50). Therefore, testing the phylogenetic signal of a molecular marker is essential for phylogeny reconstruction of closely related and cryptic species, before proposing evolutionary patterns for biological entities or delimiting species based on trees (50).
One method for estimating the phylogenetic signal is tree topology comparison using several methods of analysis. For DNA-based matrices, ML and BI are often suggested over MP because they use an explicit model of evolution, have a consistent approach to parameter estimation problems, have lower variance problems, do not violate many of the assumptions of the evolutionary model, and better account for branch length (51). In addition, ML can be more accurate than BI due to possible long branch attraction in trees (52,53). Our results suggest a consistent topology from both approaches (ML and BI) but a weak phylogenetic signal to differentiate Central America A. neivai specimens from those in Colombia.
For anopheline phylogeny, the use of a 460 bp fragment of Cytb usually lacks meaningful phylogenetic signal due to extreme conservation at the protein level and rapid saturation of synonymous positions for reconstructing the phylogeny among eight subgenera in Anopheles (Anopheles, Cellia, Nyssorhynchus, Kerteszia, Stethomyia, Bironellaand Chagasia) (13). However, the use of a fragment of this gene, as with additional loci (RNASer-IG1ND1) in this work, improved the estimation of tree topology in A. neivai.
In addition, transfer RNA is one of the most central and ancient molecules of the cell (54). Theevolution of this molecule supports the hypothesis of accretion of its structural parts where acceptor arm and T?C loop are known to be more ancient than the anticodon, the DHU and the variable loop (55). Because of its importance in cell protein biosynthesis, the tRNA secondary structure is highly conserved at DHU and T?C loops (56). However, polymorphisms at stems caused by a switch from a Watson-Crick base pairing (A-U, or G-C) to nonWatson-Crick pairing (A-C or G-U), could affect the fitness of individuals (57) . Our results from tRNASer secondary structure revealed there is a switch from G-C (variant B) to G-U (variant A) at the DHU stem. In evolving lineages, the fixation at population level of intermediates from this switch phenomenon is rare (56). However, the switch from Watson-Crick to non-Watson-Crick base pairing can compromise the development of bristle and decrease fecundity in other insects as reported in Drosophila (58).
The use of this region in anophelines or other culicids is scarce, despite the efficacy of tRNASer for resolving cryptic species in insects such as sandflies Lutzomyia pia and L. tihuilensis (Diptera: Psychodidae) (59) and differentiation among satyr butterflies (29). From our results, the polymorphism in A. neivai from a single transition occurring at the DHU stem cannot be considered different biological entities but instead a putative ancestral variant (B) and a derived intermediate form (variant A), as variant B shares the same Watson-Crick pairing at DHU stem with variant C from A. pholidotus.
In insect mitogenomes, a small conserved region between tRNASer and ND1 has been reported previously (35,60,61). For anophelines and other insects the size of this region can vary as result of insertions and deletions. The phylogenetic utility of IG1for resolving species complexes was demonstrated when comparing mutations (insertions) between A. gambiae and A. quadrimaculatus (62). However, the current results suggest that substitution patterns and indels in IG1 make it unsuitable for phylogenetic reconstruction and molecular taxonomy as consequence of length (25 bp) and low polymorphism. In addition, the motif (ATACTAA) for A. neivai and A. pholidotus corresponds to the binding site of the transcription termination peptide (MtTERM), signifying the end of the major strand coding region (63,64) shared among insect mtDNA (63).
Currently, the use of integrative taxonomy is multifaceted, using morphological and genetic evidence for the delimitation of the biological entities and integrating different perspectives, such as molecular systematics, phylogeography, comparative biology, and population genetics (65). From our perspective, the utility of the Cytb-SertRNA-IG1-ND1 region revealed interesting polymorphisms in A. neivai from a wide geographic range from Guatemala, the type locality in Panamá, the Pacific coast and the Urabá gulf in Colombia. However, the variations were not considered adequate to uncover cryptic species. Several variants may be adaptive and appear to be distributed in relation to altitude (Chiquimula, Guatemala) and geographic distribution (comparing Panamá, Guatemala and Colombia). The significance of these genetic variants requires further verification by additional studies on biting behavior, human affinity, mitogenomics and population dynamics.