Seasonal dynamics and molecular differentiation of three natural Anopheles species (Diptera: Culicidae) of the Maculatus group (Neocellia series) in malaria hotspot villages of Thailand

Background Anopheles sawadwongporni Rattanarithikul & Green, Anopheles maculatus Theobald and Anopheles pseudowillmori (Theobald) of the Anopheles maculatus group (Diptera: Culicidae) are recognized as potential malaria vectors in many countries from the Indian subcontinent through Southeast Asia to Taiwan. A number of malaria vectors in malaria hotspot areas along the Thai-Myanmar border belong to this complex. However, the species distribution and dynamic trends remain understudied in this malaria endemic region. Methods Mosquitoes of the Maculatus group were collected using CDC light traps every other week from four villages in Tha Song Yang District, Tak Province, Thailand from January to December 2015. Adult female mosquitoes were morphologically identified on site using taxonomic keys. Molecular species identification was performed by multiplex PCR based on the internal transcribed spacer 2 (ITS2) region of ribosomal DNA (rDNA) and sequencing of the cox1 gene at a DNA barcoding region in a subset of 29 specimens. Results A total of 1328 An. maculatus (sensu lato) female mosquitoes were captured with An. maculatus, An. sawadwongporni and An. pseudowilmori accounting for 75.2, 22.1 and 2.7% respectively. The field captured mosquitoes of the Maculatus group were most abundant in the wet season and had a preferred distribution in villages at higher elevations. The phylogenetic relationships of 29 cox1 sequences showed a clear-cut separation of the three member species of the Maculatus group, with the An. pseudowillmori cluster being separated from An. sawadwongporni and An. maculatus. Conclusions This study provides updated information for the species composition, seasonal dynamics and microgeographical distribution of the Maculatus group in malaria-endemic areas of western Thailand. This information can be used to guide the planning and implementation of mosquito control measures in the pursuance of malaria transmission.


Background
Malaria is a serious parasitic disease caused by the protozoan parasites of the genus Plasmodium that occurs mainly throughout tropical and subtropical zones [1]. In 2018, the World Health Organization estimated that there were 228 million cases of malaria, down from 231 million cases in 2017, and 405,000 deaths from malaria worldwide [2]. Anopheles mosquitoes (Diptera: Culicidae) are the sole vectors of malaria. Globally, a total of 494 species of the subfamily Anophelinae are currently recognized, which are divided into the following eight subgenera: Anopheles (190 species), Baimaia (1 species), Cellia (225 species), Christya (2 species), Kerteszia (12 species), Lophopodomyia (6 species), Nyssorhynchus (40 species), and Stethomyia (5 species) [3]. Approximately 70 species of formally recognized Anopheles are human malaria vectors and 40 species are dominant malaria vectors [4,5].
Different species of Anopheles mosquitoes play different roles in malaria transmission in each area. Each Anopheles species has particular characteristics in biology, ecology, and behavior such as host preference, biting indoors or outdoors (endophagic/exophagic), resting behavior (exophilic/endophilic), longevity, and larval habitat preference [6]. Therefore, the first key for the control of malaria vectors is the knowledge of the correct vector species that can transmit malaria. However, species identification of Anopheles using standard morphological methods is difficult in those with similar morphologies including the various species complexes in Southeast Asia [7]. Additionally, the collection of adult mosquitoes in the field using traps often causes damage to the specimens, resulting misidentification [8].
The Maculatus group as a large species assemblage consists of nine member species including Anopheles dispar Rattanarithikul [3,9]. The member species in this group are distributed throughout many countries from the Indian subcontinent through Southeast Asia to Taiwan and are recognized as being associated with malaria [7,10,11]. However, microgeographical distribution and seasonal dynamics of this species complex are unclear, which are an obstacle for their control. Since each member species has similar morphological characteristics, species identification using morphological methods is difficult [7].
In Thailand, seven member species of the Maculatus group have been reported, of which three (An. sawadwongporni, An. maculatus and An. pseudowillmori) are the potential malaria vectors [12][13][14]. In malaria transmission hotspot areas of Thailand such as Tak Province, these species are found in sympatry [12][13][14][15]. Although all these three Anopheles species have similar morphological characteristics, An. maculatus and An. sawadwongporni can be distinguished by a tuft of black scales on the wing vein bifurcation of radius 2 and radius 3, while An. pseudowillmori can be identified by the pattern of scales on the abdomen and wings [9]. Although species identification during routine entomological surveillance can be performed down to the subgenus level, resolution of the species complex requires in-depth morphological analysis by an expert entomologist. Damage of the wing scales, especially the characters needed for species identification, often occurs during sample collection in the rainy season. Thus, molecular techniques are a powerful tool for the identification of morphologically similar species or morphologically indistinguishable species, and have become increasingly popular [7,16,17]. Such methods include PCR, which is relatively quick, straightforward, accurate, and reliable [18]. Sequence markers that have been used in previous research for the identification of mosquito species include the internal transcribed spacer 2 (ITS2) of ribosomal DNA (rDNA) [16,18,19], the second and third domains (D2 and D3) of the rDNA 28S gene [16], and cytochrome c oxidase subunits 1 and 2 (cox1 and cox2) of the mitochondrial DNA [7,20]. The cox1 gene is considered one of the most effective molecular markers because of frequent base substitutions in the third codon position, a lack of introns, limited exposure to recombination, and a haploid mode of inheritance. It is therefore commonly used in barcoding studies for the species identification of mosquito vectors [21,22]. In Thailand, cox1 was used for DNA barcoding identification of sibling species of Anopheles mosquitoes including the Barbirostris group [22], the Hyrcanus group [23], and the Leucosphyrus group [24]. However, the identification of Anopheles species using the cox1 gene alone is not sufficient to make precise conclusions, and thus should be combined with alternative markers to increase the accuracy. Walton et al. [18] developed an effective PCR-based identification method to distinguish five member species of the Maculatus group in Thailand.
The present study aims to provide updated information on the natural distribution of the An. maculatus group on the village scale and assess the taxonomic status of naturally captured species of the Maculatus group in a malaria endemic area in Tak Province, Thailand using two molecular approaches: multiplex PCR of ITS2 and sequencing of the DNA barcode region in the cox1 gene.
The sample sources were also used to define the vector seasonal dynamics and the composition of the Maculatus complex. The results are consistent with the distribution of the Maculatus group in high-transmission area implicating their contribution to malaria transmission. Bottle programmable rotator CDC light traps (Bio-Quip) were set to collect the mosquitoes every 2 h for 6 replacements/week/village. They were placed outdoors from additional 24 selected houses in each village. In the morning, mosquito samples were transferred from the CDC light traps into 60 ml specimen containers, kept in dry ice boxes, and sent to the Department of Medical Entomology, Faculty of Tropical Medicine, Mahidol University, Bangkok for species identification. Sources of samples were recorded with details of the trap types, location, gravidity, and geographical information at each village. Adult females of the maculatus group were morphologically identified on site using taxonomic keys [9,25]. Maculatus group samples were focused and reconfirmed in the laboratory by morphology. Samples were then preserved dry at -20 °C for molecular identification. The meteorological parameters of temperature, humidity, and rainfall were obtained from the local climatology division (code station 376202), Meteorological Department, Ministry of Information and Communication Technology, Bangkok, Thailand. The station is localized in Mae Sod District, Tak Province, about 60 km away from the study site.

Malaria parasite detection in the mosquitoes
Collected Anopheles mosquitoes were kept at -20 °C until the detection of malaria parasite sporozoites by using an enzyme-linked immunosorbent assay (ELISA).

DNA extraction and polymerase chain reaction (PCR)
For a subset of 29 Anopheles specimens that were confirmed by expert entomologists using morphological characteristics, including nine An. sawadwongporni individuals (SAW1-9), nine An. maculatus individuals (MAC1-9), and 11 An. pseudowillmori individuals (PSE1-11), genomic DNA was extracted using the Pure-Link ® Genomic DNA Mini Kit (Invitrogen, Carlsbad, USA) according to the manufacturer's protocol. DNA was stored in separate 1.5 ml microcentrifuge tubes at -20 °C until further analysis. Information on the 29 mosquito specimens is provided in Additional file 1: Table S1.

Sequence alignment and phylogenetic analyses
Each cox1 sequence was compared with previously published sequences in GenBank using the standard nucleotide Basic Logical Alignment Search Tool, available at https ://blast .ncbi.nlm.nih.gov/Blast .cgi and in The Barcode of Life Database (BOLD) available at https ://www.barco dingl ife.org. The chromatograms (traces) of cox1 sequences were visualized and manually edited using Chromas software version 2.6.6 (https ://techn elysi um.com.au/wp/). Multiple cox1 sequence alignment was performed using Clustal X version 2.1 [28]. Intraspecific and interspecific pairwise sequence divergences of all individuals were calculated using the distance model Kimura 2-parameter (K2P), available within the Molecular Evolutionary Genetics Analysis (MEGA) software version 7 based on the number of nucleotide substitutions per site between two DNA sequences [29]. The detailed sequences of 29 individuals of three members of the Maculatus group were submitted and uploaded to GenBank (Additional file 2: Table S2). A phylogenetic tree was built to provide a graphical representation of the clustering pattern among member species of the Maculatus group. To identify the genetic relationships among the Maculatus group by phylogenetic analysis, we used the reference genome for other member species in the Maculatus group from GenBank (Additional file 2: Table S2). The sequences of Aedes aegypti, Ae. albopictus, Culex quinquefasciatus, and Mansonia bonneae were used as outgroup species (Additional file 2: Table S2) . Phylogenetic analysis was performed by using the maximum likelihood method with 1000 bootstrap replications implemented in the MEGA 7.0 software. The general time reversible plus gamma distributed with invariant sites (GTR + G + I) was used as the nucleotide substitution model.

Morphological identification of An. maculatus group mosquitoes
In 2015, a total of 7519 Anopheles mosquitoes were collected using the two different light traps. Twenty-four species were identified morphologically, among which An. minimus, An. culicifacies, An. maculatus, An. annularis were the most abundant, accounting for 49.5, 14.0, 13.3 and 11.3% of the total collection, respectively (Table 1). In total, 1328 An. maculatus (s.l.) were captured from the four villages. Two species of the Maculatus group, An. maculatus and An. sawadwongporni, were initially assigned on site by the field sample collection team using morphological characters. However, morphological reconfirmation by expert entomologists in the medical entomology laboratory revealed that 75.2, 22.1 and 2.7% of these specimens were An. maculatus, An. sawadwongporni and An. pseudowilmori, respectively ( Table 2). The monthly dynamics of the Maculatus group captured from the four villages showed two peaks during the rainy season with the main peak occurring in August to September (Fig. 2a). While most samples of this species group were collected from the high-transmission sites (KN and SO), An. maculatus was the major species captured from both high-and low-transmission villages. KN village had the highest abundance of the Maculatus group with the highest proportion of An. maculatus 58.6% (778/1328), which peaked in August. In other villages, the captured An. maculatus ranged from 4.4% to 6.7%. Anopheles sawadwongporni had the same seasonal trend with high prevalence during June-August. The majority (13.10%, 174/1328) of the An. sawadwongporni samples was collected in KN, whereas 1.7-13.1% were collected in other villages. Anopheles pseudowilmori was found in three of the four villages except SO, with the highest numbers of 17 in KN followed by 12 in NB (Fig. 2b). The rotating CDC light traps were used to monitor the times of mosquito activities. The results showed that the first peak of the Maculatus group occurred at 20.00-22.00 h and the second peak at 00.00-2.00 h (Fig. 2c). Of all of the Maculatus species captured, 48.9% An. maculatus, 15.3% An. sawadwongporni, and 2.1% An. pseudowilmori, respectively, were collected by the outdoor traps (Fig. 2d, Table 2). Of the 1328 Maculatus species collected, 25 (1.9%) were blood-fed. For An. maculatus, 2.3% (20/999) were blood-fed, and the ratio of blood-fed An. maculatus (8.6%) was the highest in NB (Table 2).

Plasmodium infectivity in the Maculatus group mosquitoes
ELISA analysis of all mosquitoes of the Maculatus group revealed that only one An. maculatus mosquito captured from an indoor light trap in KN village in August was positive for P. vivax (PV210).

Molecular identification based on the species-specific ITS2 region
The 29 samples that were identified to the species level using morphological characteristics were investigated by multiplex PCR based on the ITS2 region. All Anopheles species in the Maculatus group were confirmed to match their morphological confirmation and cox1 sequence identification. The 9 An. sawadwongporni displayed the largest PCR band of 242 bp, followed by 11 An. pseudowillmori (203 bp), and 9 An. maculatus (180 bp) (Fig. 3).

Mitochondrial cox1 gene sequence analysis
The cox1 sequences of the three species in the Maculatus group (658 bp) were analyzed between and within species. No noise or double peaks were found in the chromatograms of all obtained sequences. The adenosine and thymine (AT) content was similar for the three species: 69% in An. sawadwongporni; 68.5% in An. maculatus; and 67.8% in An. pseudowillmori. These sequences had the highest average A + T content at the first codon position, including 93.6% in An. sawadwongporni, 91.1% in An. maculatus, and 90.5% in An. pseudowillmori. The average intraspecific divergence of An. sawadwongporni, An. maculatus, and An. pseudowillmori based on K2P distances was 0.002, 0.004, and 0.008, respectively (Table 3), and the overall divergence among the three species was 0.069. Average interspecific divergences were 0.069 between An. sawadwongporni and An. maculatus, 0.114 between An. sawadwongporni and An. pseudowillmori, and 0.104 between An. maculatus and An. pseudowillmori (Table 3).
For species identification, sequences from the three members in the Maculatus group were compared with available sequences in the GenBank and BOLD systems. The 29 cox1 sequences generated from the three species of morphologically identified specimens showed more than 99% homology to An. maculatus (GenBank: JQ728164), An. sawadwongporni (GenBank: Q728408), and An. pseudowillmori (GenBank: JQ728241), respectively.

Phylogenetic analysis
The relationships among the 29 cox1 sequences of An. sawadwongporni, An. maculatus and An. pseudowillmori of the Maculatus group in Thailand and 10 cox1 sequences retrieved from the GenBank database were analyzed by phylogeny (Fig. 4). The maximum likelihood tree showed a clear-cut separation between species members of the Maculatus group clusters. All specimens of the same species were grouped in individual clusters. The addition of a cox1 sequence of An. dravidicus from GenBank showed that it also belongs to the Maculatus group. The An. pseudowillmori cluster was separated from the An. dravidicus cluster and An. sawadwongporni and An. maculatus sub-cluster. Additionally, outgroup species including Ae. aegypti, Ae. albopictus, Cx. quinquefasciatus, and Ma. bonneae were clearly separated from clusters of Anopheles species in the Maculatus group. Altogether, the results of cox1 sequence analysis including the genetic divergence of K2P distances and phylogenetic tree clearly showed consistency with species identification derived from the species-specific ITS2 and morphological approach from the re-examination by entomological experts.

Discussion
Anopheles maculatus, An. sawadwongporni and An. pseudowillmori are recognized as potential malaria vectors in many regions of the world including Thailand [13,30], Motuo County in China's Tibet [31], and southwestern China [32]. This study reported the dynamics of the Maculatus group in western Thailand, which is consistent with their malaria vector status with peak abundance in the wet season [13]. Herein, the study confirmed that there are at least three members of the Maculatus group  with An. maculatus being the predominant potential malaria vector in this region. Although the same members of the Maculatus group were found in the neighboring Kayin State, Myanmar, An. swadwongporni was the main species there, suggesting species abundance within a small geographical region may vary [33]. This study further revealed that even within a small geographical area such as the four villages, vector abundance can vary drastically. More than 58% of the mosquito collection of the Maculatus group was from the village KN, which is near the forest and surrounded by mountains. The identification of a sporozoite-positive An. maculatus from KN also supports active malaria transmission in this village. Since mosquito control efforts consisting of indoor sprays of insecticides twice a year were conducted similarly in these villages by the local Vector-Borne Disease Control Unit, the differences in vector abundance may be due to the variations in local environment. With this information, more intensive mosquito control efforts should be undertaken in the village KN. This study identified that the Maculatus group are generally more exophilic, though there were slight variations in the endophilic patterns of the different species within this group (Fig. 2d) as confirmed in another site of the Myanmar-Thai border [33]. Although the endoor exophilicity is referred from the sites where they are collected (indoor vs outdoor) and may not reflect the true preference, this information nonetheless should be useful to guide the implementation of outdoor control efforts to prevent outdoor malaria transmission. Another important finding is the relatively persistent activity of the Maculatus group throughout the collection time in the night, though two peak collection times were noticed at 20.00-22.00 h and 00.00-2.00 h. The first peak time coincides with the family time before bed, suggesting that LLINs may offer ineffective protection for these earlybiting mosquitoes. Alternative measures preventing mosquito biting before bedtime are needed [34].
Discriminating Anopheles species within species complexes can be difficult because of morphological similarity. Molecular techniques have been used to overcome this difficulty. This includes amplification of the ITS2 region of rDNA [18], D3 of rDNA, and cox2 of mitochondrial DNA [16]. cox1 could be used for DNA barcoding in mosquito identification [17], but has not been evaluated properly among member species of the Maculatus group. The present study determined the efficiency of using cox1 barcodes for the identification of An. sawadwongporni, An. maculatus, and An. pseudowillmori in Thailand, and established a group of sequences associated with the identified species. Species-specific multiplex PCR based on the ITS2 region was also used to confirm the findings based on cox1 analysis.
Phylogenetic relationships based on cox1 sequences clearly revealed grouping between species in the Maculatus group. The phylogeny based on this gene indicated that An. sawadwongporni is more similar to An. maculatus, while An. pseudowillmori is more distantly related to these two species [17]. These results were consistent with the analysis of other molecular markers in previous studies including ITS2 of rDNA [11,18] and combined ITS2 and D3 [16], which separated the An. pseudowillmori cluster from the An. sawadwongporni and An. maculatus clusters. Recently, these three species were identified by analyzing wing geometry which revealed consistent findings with DNA analysis, while modern morphometric findings were consistent with cox1 phylogenetic analysis and wing morphology [15]. Again, this finding based on the cox1 results is consistent with the current mosquito taxonomy, which placed An. pseudowillmori in the Maculatus group, and An. sawadwongporni and An. maculatus in the Sawadwongporni and Maculatus subgroups, respectively [9].
Our multiplex PCR results based on ITS2 were consistent with the cox1 gene findings which identified species according to their morphological identification. This confirms the accuracy of species identification of the Maculatus group by cox1 DNA barcodes. The cox1 sequences of the three potential malaria vectors in Thailand submitted to the GenBank database could be used as reference cox1 sequences in the understudied Maculatus group and for future taxonomic studies of Anopheles vectors.